StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
dAuBackground.cxx
1 #include <stdio.h>
2 #include <Riostream.h>
3 #include <TFile.h>
4 #include <TTree.h>
5 #include <TH2F.h>
6 #include <TPostScript.h>
7 #include <TCanvas.h>
8 #include <TVector3.h>
9 #include <TPad.h>
10 #include <TStyle.h>
11 #include <TString.h>
12 #include <TPaveText.h>
13 #include <TMath.h>
14 
15 #include <StEmcPool/StPhotonCommon/MyEvent.h>
16 #include <StEmcPool/StPhotonCommon/MyPoint.h>
17 #include <StEmcPool/StPhotonCommon/MyMcTrack.h>
18 
19 #include "AnaCuts.h"
20 
21 #include "dAuBackground.h"
22 using namespace std;
23 
24 ClassImp(dAuBackground)
25 
26 Int_t isBad[2400];
27 
28 dAuBackground::dAuBackground()
29 {
30  cuts = 0;
31 
32  gStyle->SetPadGridX(kFALSE);
33  gStyle->SetPadGridY(kFALSE);
34 
35  for(Int_t i=0;i<2400;i++)
36  {
37  isBad[i]=0;
38  }
39  Int_t badIds[]={82,156,240,264,733,756,818,853,931,1289,1292,1355,1359,1527,1759,1943,2011,2025,2090,2152,2369,-1};
40  Int_t j=0;
41  while(badIds[j]>0)
42  {
43  isBad[badIds[j]-1]=1;
44  j++;
45  }
46 
47  cuts = new AnaCuts();
48  cout<<"constructed!!"<<endl;
49 }
50 dAuBackground::~dAuBackground()
51 {
52  if (cuts) delete cuts; cuts = 0;
53  cout<<"destructed!!"<<endl;
54 }
55 
56 void dAuBackground::run(const char *file)
57 {
58  //0-20
59  TH1F *ratio_mb_020=new TH1F("ratio_mb_020","2001&&2003, FTPC-Au:0-20%", 25,0.,1.);
60  TH1F *ratio_ht1_020=new TH1F("ratio_ht1_020","2201, FTPC-Au:0-20%", 25,0.,1.);
61  TH1F *ratio_ht2_020=new TH1F("ratio_ht2_020","2202, FTPC-Au:0-20%", 25,0.,1.);
62  TH1F *ratio_sn_020=new TH1F("ratio_sn_020","2002, FTPC-Au:0-20%", 25,0.,1.);
63 
64  //20-40
65  TH1F *ratio_mb_2040=new TH1F("ratio_mb_2040","2001&&2003, FTPC-Au:20-40%", 25,0.,1.);
66  TH1F *ratio_ht1_2040=new TH1F("ratio_ht1_2040","2201, FTPC-Au:20-40%", 25,0.,1.);
67  TH1F *ratio_ht2_2040=new TH1F("ratio_ht2_2040","2202, FTPC-Au:20-40%", 25,0.,1.);
68  TH1F *ratio_sn_2040=new TH1F("ratio_sn_2040","2002, FTPC-Au:20-40%", 25,0.,1.);
69 
70  //40-100
71  TH1F *ratio_mb_40100=new TH1F("ratio_mb_40100","2001&&2003, FTPC-Au:40-100%", 25,0.,1.);
72  TH1F *ratio_ht1_40100=new TH1F("ratio_ht1_40100","2201, FTPC-Au:40-100%", 25,0.,1.);
73  TH1F *ratio_ht2_40100=new TH1F("ratio_ht2_40100","2202, FTPC-Au:40-100%", 25,0.,1.);
74  TH1F *ratio_sn_40100=new TH1F("ratio_sn_40100","2002, FTPC-Au:40-100%", 25,0.,1.);
75 
76  TH1F *hfrac_mb=new TH1F("frac_mb","2001&&2003, I[0.8,1.0]/I[0.0,1.0], 40-100%",50,0.0,0.05);
77  TH1F *hfrac_ht1=new TH1F("frac_ht1","2201, I[0.8,1.0]/I[0.0,1.0], 40-100%",50,0.0,0.5);
78  TH1F *hfrac_ht2=new TH1F("frac_ht2","2202, I[0.8,1.0]/I[0.0,1.0], 40-100%",50,0.0,0.5);
79  TH1F *hfrac_sn=new TH1F("frac_sn","2000, I[0.8,1.0]/I[0.0,1.0], 40-100%",50,0.0,0.05);
80 
81  //totals
82  //0-20
83  TH1F *Ratio_mb_020=new TH1F("Ratio_mb_020","2001&&2003, FTPC-Au:0-20%", 25,0.,1.);
84  TH1F *Ratio_ht1_020=new TH1F("Ratio_ht1_020","2201, FTPC-Au:0-20%", 25,0.,1.);
85  TH1F *Ratio_ht2_020=new TH1F("Ratio_ht2_020","2202, FTPC-Au:0-20%", 25,0.,1.);
86  TH1F *Ratio_sn_020=new TH1F("Ratio_sn_020","2002, FTPC-Au:0-20%", 25,0.,1.);
87 
88  //20-40
89  TH1F *Ratio_mb_2040=new TH1F("Ratio_mb_2040","2001&&2003, FTPC-Au:20-40%", 25,0.,1.);
90  TH1F *Ratio_ht1_2040=new TH1F("Ratio_ht1_2040","2201, FTPC-Au:20-40%", 25,0.,1.);
91  TH1F *Ratio_ht2_2040=new TH1F("Ratio_ht2_2040","2202, FTPC-Au:20-40%", 25,0.,1.);
92  TH1F *Ratio_sn_2040=new TH1F("Ratio_sn_2040","2002, FTPC-Au:20-40%", 25,0.,1.);
93 
94  //40-100
95  TH1F *Ratio_mb_40100=new TH1F("Ratio_mb_40100","2001&&2003, FTPC-Au:40-100%", 25,0.,1.);
96  TH1F *Ratio_ht1_40100=new TH1F("Ratio_ht1_40100","2201, FTPC-Au:40-100%", 25,0.,1.);
97  TH1F *Ratio_ht2_40100=new TH1F("Ratio_ht2_40100","2202, FTPC-Au:40-100%", 25,0.,1.);
98  TH1F *Ratio_sn_40100=new TH1F("Ratio_sn_40100","2002, FTPC-Au:40-100%", 25,0.,1.);
99 
100  TH2F *EvsE_mb=new TH2F("EvsE_mb","2001&&2003, neutral Energy vs TPC pT",240,0.,80,240,0.,80.);
101  TH2F *EvsE_ht1=new TH2F("EvsE_ht1","2201, neutral Energy vs TPC pT",124,0.,80,240,0.,80.);
102  TH2F *EvsE_ht2=new TH2F("EvsE_ht2","2202, neutral Energy vs TPC pT",124,0.,80,240,0.,80.);
103  TH2F *EvsE_sn=new TH2F("EvsE_sn","2002, neutral Energy vs TPC pT",240,0.,80,240,0.,80.);
104 
105  //spectra
106  //0-20
107  TH1F *ptspec_mb_020=new TH1F("ptspec_mb_020","2001&&2003, pT spectrum neutral points, 0-20%",30,0.,15.);
108  TH1F *ptspec_ht1_020=new TH1F("ptspec_ht1_020","2201, pT spectrum neutral points, 0-20%",30,0.,15.);
109  TH1F *ptspec_ht2_020=new TH1F("ptspec_ht2_020","2202, pT spectrum neutral points, 0-20%",30,0.,15.);
110  TH1F *ptspec_sn_020=new TH1F("ptspec_sn_020","2002, pT spectrum neutral points, 0-20%",30,0.,15.);
111 
112  TH1F *ptspec_mb_cut_020=new TH1F("ptspec_mb_cut_020","2001&&2003, pT spectrum neutral points after cut, 0-20%",30,0.,15.);
113  TH1F *ptspec_ht1_cut_020=new TH1F("ptspec_ht1_cut_020","2201, pT spectrum neutral points after cut, 0-20%",30,0.,15.);
114  TH1F *ptspec_ht2_cut_020=new TH1F("ptspec_ht2_cut_020","2202, pT spectrum neutral points after cut, 0-20%",30,0.,15.);
115  TH1F *ptspec_sn_cut_020=new TH1F("ptspec_sn_cut_020","2002, pT spectrum neutral points after cut, 0-20%",30,0.,15.);
116 
117  TH1F *ptspec_mb_bg_020=new TH1F("ptspec_mb_bg_020","2001&&2003, pT spectrum neutral points, 0-20%",30,0.,15.);
118  TH1F *ptspec_ht1_bg_020=new TH1F("ptspec_ht1_bg_020","2201, pT spectrum neutral points, 0-20%",30,0.,15.);
119  TH1F *ptspec_ht2_bg_020=new TH1F("ptspec_ht2_bg_020","2202, pT spectrum neutral points, 0-20%",30,0.,15.);
120  TH1F *ptspec_sn_bg_020=new TH1F("ptspec_sn_bg_020","2002, pT spectrum neutral points, 0-20%",30,0.,15.);
121 
122  //20-40
123  TH1F *ptspec_mb_2040=new TH1F("ptspec_mb_2040","2001&&2003, pT spectrum neutral points, 20-40%",30,0.,15.);
124  TH1F *ptspec_ht1_2040=new TH1F("ptspec_ht1_2040","2201, pT spectrum neutral points, 20-40%",30,0.,15.);
125  TH1F *ptspec_ht2_2040=new TH1F("ptspec_ht2_2040","2202, pT spectrum neutral points, 20-40%",30,0.,15.);
126  TH1F *ptspec_sn_2040=new TH1F("ptspec_sn_2040","2002, pT spectrum neutral points, 20-40%",30,0.,15.);
127 
128  TH1F *ptspec_mb_cut_2040=new TH1F("ptspec_mb_cut_2040","2001&&2003, pT spectrum neutral points after cut, 20-40%",30,0.,15.);
129  TH1F *ptspec_ht1_cut_2040=new TH1F("ptspec_ht1_cut_2040","2201, pT spectrum neutral points after cut, 20-40%",30,0.,15.);
130  TH1F *ptspec_ht2_cut_2040=new TH1F("ptspec_ht2_cut_2040","2202, pT spectrum neutral points after cut, 20-40%",30,0.,15.);
131  TH1F *ptspec_sn_cut_2040=new TH1F("ptspec_sn_cut_2040","2002, pT spectrum neutral points after cut, 20-40%",30,0.,15.);
132 
133  TH1F *ptspec_mb_bg_2040=new TH1F("ptspec_mb_bg_2040","2001&&2003, pT spectrum neutral points, 20-40%",30,0.,15.);
134  TH1F *ptspec_ht1_bg_2040=new TH1F("ptspec_ht1_bg_2040","2201, pT spectrum neutral points, 20-40%",30,0.,15.);
135  TH1F *ptspec_ht2_bg_2040=new TH1F("ptspec_ht2_bg_2040","2202, pT spectrum neutral points, 20-40%",30,0.,15.);
136  TH1F *ptspec_sn_bg_2040=new TH1F("ptspec_sn_bg_2040","2002, pT spectrum neutral points, 20-40%",30,0.,15.);
137 
138  //40-100
139  TH1F *ptspec_mb_40100=new TH1F("ptspec_mb_40100","2001&&2003, pT spectrum neutral points, 40-100%",30,0.,15.);
140  TH1F *ptspec_ht1_40100=new TH1F("ptspec_ht1_40100","2201, pT spectrum neutral points, 40-100%",30,0.,15.);
141  TH1F *ptspec_ht2_40100=new TH1F("ptspec_ht2_40100","2202, pT spectrum neutral points, 40-100%",30,0.,15.);
142  TH1F *ptspec_sn_40100=new TH1F("ptspec_sn_40100","2002, pT spectrum neutral points, 40-100%",30,0.,15.);
143 
144  TH1F *ptspec_mb_cut_40100=new TH1F("ptspec_mb_cut_40100","2001&&2003, pT spectrum neutral points after cut, 40-100%",30,0.,15.);
145  TH1F *ptspec_ht1_cut_40100=new TH1F("ptspec_ht1_cut_40100","2201, pT spectrum neutral points after cut, 40-100%",30,0.,15.);
146  TH1F *ptspec_ht2_cut_40100=new TH1F("ptspec_ht2_cut_40100","2202, pT spectrum neutral points after cut, 40-100%",30,0.,15.);
147  TH1F *ptspec_sn_cut_40100=new TH1F("ptspec_sn_cut_40100","2002, pT spectrum neutral points after cut, 40-100%",30,0.,15.);
148 
149  TH1F *ptspec_mb_bg_40100=new TH1F("ptspec_mb_bg_40100","2001&&2003, pT spectrum neutral points, 40-100%",30,0.,15.);
150  TH1F *ptspec_ht1_bg_40100=new TH1F("ptspec_ht1_bg_40100","2201, pT spectrum neutral points, 40-100%",30,0.,15.);
151  TH1F *ptspec_ht2_bg_40100=new TH1F("ptspec_ht2_bg_40100","2202, pT spectrum neutral points, 40-100%",30,0.,15.);
152  TH1F *ptspec_sn_bg_40100=new TH1F("ptspec_sn_bg_40100","2002, pT spectrum neutral points, 40-100%",30,0.,15.);
153 
154 
155  TH1F *ftpc_mb=new TH1F("ftpc_mb","2001&&2003, ftpc-Au mult.",60,0.,60.);
156  TH1F *ftpc_ht1=new TH1F("ftpc_ht1","2201, ftpc-Au mult.",60,0.,60.);
157  TH1F *ftpc_ht2=new TH1F("ftpc_ht2","2202, ftpc-Au mult.",60,0.,60.);
158  TH1F *ftpc_sn=new TH1F("ftpc_sn","2002, ftpc-Au mult.",60,0.,60.);
159 
160  TH1F *ftpc_mb_cut=new TH1F("ftpc_mb_cut","2001&&2003, ftpc-Au mult. for r<.8",60,0.,60.);
161  TH1F *ftpc_ht1_cut=new TH1F("ftpc_ht1_cut","2201, ftpc-Au mult. for r<.8",60,0.,60.);
162  TH1F *ftpc_ht2_cut=new TH1F("ftpc_ht2_cut","2202, ftpc-Au mult. for r<.8",60,0.,60.);
163  TH1F *ftpc_sn_cut=new TH1F("ftpc_sn_cut","2002, ftpc-Au mult. for r<.8",60,0.,60.);
164 
165  TH1F *tpc_mb=new TH1F("tpc_mb","2001&&2003, tpc mult.",100,0.,100.);
166  TH1F *tpc_ht1=new TH1F("tpc_ht1","2201, tpc mult.",100,0.,100.);
167  TH1F *tpc_ht2=new TH1F("tpc_ht2","2202, tpc mult.",100,0.,100.);
168  TH1F *tpc_sn=new TH1F("tpc_sn","2002, tpc mult.",100,0.,100.);
169 
170  TH1F *tpc_mb_cut=new TH1F("tpc_mb_cut","2001&&2003, tpc mult. for r<.8",100,0.,100.);
171  TH1F *tpc_ht1_cut=new TH1F("tpc_ht1_cut","2201, tpc mult. for r<.8",100,0.,100.);
172  TH1F *tpc_ht2_cut=new TH1F("tpc_ht2_cut","2202, tpc mult. for r<.8",100,0.,100.);
173  TH1F *tpc_sn_cut=new TH1F("tpc_sn_cut","2002, tpc mult. for r<.8",100,0.,100.);
174 
175  TH1F *bemc_mb=new TH1F("bemc_mb","2001&&2003, bemc mult.",60,0.,60.);
176  TH1F *bemc_ht1=new TH1F("bemc_ht1","2201, bemc mult.",60,0.,60.);
177  TH1F *bemc_ht2=new TH1F("bemc_ht2","2202, bemc mult.",60,0.,60.);
178  TH1F *bemc_sn=new TH1F("bemc_sn","2002, bemc mult.",60,0.,60.);
179 
180  TH1F *bemc_mb_cut=new TH1F("bemc_mb_cut","2001&&2003, bemc mult. for r<.8",60,0.,60.);
181  TH1F *bemc_ht1_cut=new TH1F("bemc_ht1_cut","2201, bemc mult. for r<.8",60,0.,60.);
182  TH1F *bemc_ht2_cut=new TH1F("bemc_ht2_cut","2202, bemc mult. for r<.8",60,0.,60.);
183  TH1F *bemc_sn_cut=new TH1F("bemc_sn_cut","2002, bemc mult. for r<.8",60,0.,60.);
184 
185  TH2F *etaphi_fill3085_sn=new TH2F("etaphi_fill3085_sn","eta/phi for 3085 sn 40-100%",500,0.0,1.0,1000,-3.2,3.2);
186  TH2F *etaphi_fill3085_ht2=new TH2F("etaphi_fill3085_ht2","eta/phi for 3085 ht-2 40-100%",500,0.0,1.0,1000,-3.2,3.2);
187  TH2F *etaphi_fill3088_sn=new TH2F("etaphi_fill3088_sn","eta/phi for 3088 sn 40-100%",500,0.0,1.0,1000,-3.2,3.2);
188  TH2F *etaphi_fill3088_ht2=new TH2F("etaphi_fill3088_ht2","eta/phi for 3088 ht-2 40-100%",500,0.0,1.0,1000,-3.2,3.2);
189  TH2F *etaphi_fill3088b_ht2=new TH2F("etaphi_fill3088b_ht2","eta/phi for 3088 ht-2 40-100% and R<0.7",500,0.0,1.0,1000,-3.2,3.2);
190 
191  TH2F *etaphi_sn=new TH2F("etaphi_sn","eta/phi sn 40-100%",500,0.0,1.0,1000,-3.2,3.2);
192 
193 
194  TH1F *h_trigger=new TH1F("h_trigger","trigger ids after cut",2400,0.5,2400.5);
195 
196  ofstream fout;
197  fout.open("bgfrac.dat");
198  fout<<"pagina\trunid\tmb\tnmb\tht1\tnht1\tht2\tnht2\tsn\tnsn"<<endl<<endl;
199 
200  TFile *inf=new TFile(file,"OPEN");
201  TTree *eventTree=(TTree*)inf->Get("mEventTree");
202 
203  MyEvent *ev=new MyEvent();
204  eventTree->SetBranchAddress("branch",&ev);
205 
206  TPostScript *ps=new TPostScript("dAu_bgperrun.ps",111);
207  TCanvas *c9=new TCanvas("c9","c9",800,600);
208  c9->Divide(4,4);
209 
210  Int_t entry=0;
211  Int_t runprev=0;
212  Int_t pagina=0;
213  while(eventTree->GetEntry(entry))
214  {
215  //if(ev->runId()<4061003)
216  //{
217  //entry++;
218  // continue;
219  //}
220  //if(ev->runId()>4061025)
221  //{
222  // cout<<"end of this fill!"<<endl;
223  // break;
224  //}
225 
226  if(ev->runId()!=runprev && runprev!=0)
227  {
228  if(ev->runId()<runprev)
229  {
230  cout<<"Something wrong!"<<endl;
231  return;
232  }
233  cout<<endl<<"new run: "<<ev->runId()<<endl<<endl;
234  pagina++;
235 
236  Float_t total_mb=ratio_mb_40100->Integral();
237  Float_t bg_mb=ratio_mb_40100->Integral(21,25);
238  Float_t total_ht1=ratio_ht1_40100->Integral();
239  Float_t bg_ht1=ratio_ht1_40100->Integral(21,25);
240  Float_t total_ht2=ratio_ht2_40100->Integral();
241  Float_t bg_ht2=ratio_ht2_40100->Integral(21,25);
242  Float_t total_sn=ratio_sn_40100->Integral();
243  Float_t bg_sn=ratio_sn_40100->Integral(21,25);
244 
245  Float_t frac_mb=total_mb>0 ? bg_mb/total_mb : 0.;
246  Float_t frac_ht1=total_ht1>0 ? bg_ht1/total_ht1 : 0.;
247  Float_t frac_ht2=total_ht2>0 ? bg_ht2/total_ht2 : 0.;
248  Float_t frac_sn=total_sn>0 ? bg_sn/total_sn : 0.;
249  fout<<pagina<<"\t"<<runprev<<"\t"<<frac_mb<<"\t"<<total_mb<<"\t"<<frac_ht1<<"\t"<<total_ht1<<"\t"<<frac_ht2<<"\t"<<total_ht2<<"\t"<<frac_sn<<"\t"<<total_sn<<endl;
250 
251  hfrac_mb->Fill(frac_mb,total_mb);//use weight
252  hfrac_ht1->Fill(frac_ht1,total_ht1);
253  hfrac_ht2->Fill(frac_ht2,total_ht2);
254  hfrac_sn->Fill(frac_sn,total_sn);
255 
256  //create new ps page and plot previous run
257  ps->NewPage();
258  c9->cd(1);
259  ratio_mb_020->Draw();
260  c9->cd(2);
261  ratio_ht1_020->Draw();
262  c9->cd(3);
263  ratio_ht2_020->Draw();
264  c9->cd(4);
265  ratio_sn_020->Draw();
266  c9->cd(5);
267  ratio_mb_2040->Draw();
268  c9->cd(6);
269  ratio_ht1_2040->Draw();
270  c9->cd(7);
271  ratio_ht2_2040->Draw();
272  c9->cd(8);
273  ratio_sn_2040->Draw();
274  c9->cd(9);
275  ratio_mb_40100->Draw();
276  c9->cd(10);
277  ratio_ht1_40100->Draw();
278  c9->cd(11);
279  ratio_ht2_40100->Draw();
280  c9->cd(12);
281  ratio_sn_40100->Draw();
282  c9->cd(13);
283  TPaveText *pave=new TPaveText(.1,.1,.9,.9);
284  char label[20];
285  sprintf(label,"run %d",runprev);
286  pave->AddText(label);
287  pave->Draw();
288  c9->cd(0);
289  c9->Update();
290 
291  //reset hists
292  ratio_mb_020->Reset();
293  ratio_ht1_020->Reset();
294  ratio_ht2_020->Reset();
295  ratio_sn_020->Reset();
296  ratio_mb_2040->Reset();
297  ratio_ht1_2040->Reset();
298  ratio_ht2_2040->Reset();
299  ratio_sn_2040->Reset();
300  ratio_mb_40100->Reset();
301  ratio_ht1_40100->Reset();
302  ratio_ht2_40100->Reset();
303  ratio_sn_40100->Reset();
304 
305  //set runprev
306  runprev=ev->runId();
307 
308  }
309  else if(ev->runId()!=runprev)
310  {
311  cout<<endl<<"first run: "<<ev->runId()<<endl<<endl;
312  runprev=ev->runId();
313  }
314 
315  Int_t nTrig=0;
316  if(ev->runId()==0&&ev->trigger()==0)
317  {
318  nTrig=1;
319  if(ev->highTowerAdc()>cuts->ht1AdcCUT) nTrig+=2;
320  if(ev->highTowerAdc()>cuts->ht2AdcCUT) nTrig+=4;
321  ev->setTrigger(nTrig);
322  }
323 
324 
325  Bool_t eventOK=kTRUE;
326  TVector3 vPos=ev->vertex();
327  //apply EVENT cuts
328  if(TMath::Abs(vPos.Z())>cuts->vertexzCUT) eventOK=kFALSE;
329  if(TMath::Abs(vPos.X())>cuts->vertexxCUT) eventOK=kFALSE;
330  if(TMath::Abs(vPos.Y())>cuts->vertexyCUT) eventOK=kFALSE;
331  if(vPos.Mag()<.000001) eventOK=kFALSE;
332  if(ev->trigger()<1 && ev->getMcTrack()->id()==0) eventOK=kFALSE;
333 
334  if(ev->trigger()&2||ev->trigger()&4)
335  {
336  if(ev->highTowerEnergy()<0.0001) eventOK=kFALSE;
337  if(isBad[ev->highTowerId()-1]==1) eventOK=kFALSE;
338  }
339 
340 
341  if(!eventOK)
342  {
343  entry++;
344  continue;
345  }
346 
347  if(ev->trigger()&2||ev->trigger()&4) h_trigger->Fill(ev->highTowerId());
348 
349  Int_t t=ev->trigger();
350  Float_t en=ev->energyInBarrel();
351  Float_t pt=ev->momentumInTpc();
352  Float_t ratio=TMath::Abs(en+pt)>0 ? en/(en+pt) : -.1;
353 
354  if(ev->refMult()>17)
355  {
356  if(t&1) {ratio_mb_020->Fill(ratio); Ratio_mb_020->Fill(ratio);}
357  if(t&2) {ratio_ht1_020->Fill(ratio); Ratio_ht1_020->Fill(ratio);}
358  if(t&4) {ratio_ht2_020->Fill(ratio); Ratio_ht2_020->Fill(ratio);}
359  if(t&8) {ratio_sn_020->Fill(ratio); Ratio_sn_020->Fill(ratio);}
360  }
361  else if(ev->refMult()>11)
362  {
363  if(t&1) {ratio_mb_2040->Fill(ratio); Ratio_mb_2040->Fill(ratio);}
364  if(t&2) {ratio_ht1_2040->Fill(ratio); Ratio_ht1_2040->Fill(ratio);}
365  if(t&4) {ratio_ht2_2040->Fill(ratio); Ratio_ht2_2040->Fill(ratio);}
366  if(t&8) {ratio_sn_2040->Fill(ratio); Ratio_sn_2040->Fill(ratio);}
367  }
368  else if(ev->refMult()>=0)
369  {
370 
371  if(t&1) {ratio_mb_40100->Fill(ratio); Ratio_mb_40100->Fill(ratio);}
372  if(t&2) {ratio_ht1_40100->Fill(ratio); Ratio_ht1_40100->Fill(ratio);}
373  if(t&4) {ratio_ht2_40100->Fill(ratio); Ratio_ht2_40100->Fill(ratio);}
374  if(t&8) {ratio_sn_40100->Fill(ratio); Ratio_sn_40100->Fill(ratio);}
375  }
376 
377  if(t&1){
378  tpc_mb->Fill(ev->goodPrimaries());
379  ftpc_mb->Fill(ev->refMult());
380  bemc_mb->Fill(ev->numberOfPoints());
381  EvsE_mb->Fill(ev->momentumInTpc(),ev->energyInBarrel());
382  if(ratio<cuts->ratioCUT){
383  tpc_mb_cut->Fill(ev->goodPrimaries());
384  ftpc_mb_cut->Fill(ev->refMult());
385  bemc_mb_cut->Fill(ev->numberOfPoints());
386  }
387  }
388  if(t&2){
389  tpc_ht1->Fill(ev->goodPrimaries());
390  ftpc_ht1->Fill(ev->refMult());
391  bemc_ht1->Fill(ev->numberOfPoints());
392  EvsE_ht1->Fill(ev->momentumInTpc(),ev->energyInBarrel());
393  if(ratio<cuts->ratioCUT){
394  tpc_ht1_cut->Fill(ev->goodPrimaries());
395  ftpc_ht1_cut->Fill(ev->refMult());
396  bemc_ht1_cut->Fill(ev->numberOfPoints());
397  }
398  }
399  if(t&4){
400  tpc_ht2->Fill(ev->goodPrimaries());
401  ftpc_ht2->Fill(ev->refMult());
402  bemc_ht2->Fill(ev->numberOfPoints());
403  EvsE_ht2->Fill(ev->momentumInTpc(),ev->energyInBarrel());
404  if(ratio<cuts->ratioCUT){
405  tpc_ht2_cut->Fill(ev->goodPrimaries());
406  ftpc_ht2_cut->Fill(ev->refMult());
407  bemc_ht2_cut->Fill(ev->numberOfPoints());
408  }
409  }
410  if(t&8){
411  tpc_sn->Fill(ev->goodPrimaries());
412  ftpc_sn->Fill(ev->refMult());
413  bemc_sn->Fill(ev->numberOfPoints());
414  EvsE_sn->Fill(ev->momentumInTpc(),ev->energyInBarrel());
415  if(ratio<cuts->ratioCUT){
416  tpc_sn_cut->Fill(ev->goodPrimaries());
417  ftpc_sn_cut->Fill(ev->refMult());
418  bemc_sn_cut->Fill(ev->numberOfPoints());
419  }
420  }
421 
422  //get neutral points
423  MyPoint *p;
424  TClonesArray *clA=(TClonesArray*)ev->getPointArray();
425  for(Int_t j=0;j<clA->GetEntries();j++)
426  {
427  p=(MyPoint*)clA->At(j);
428  TVector3 pPos=p->position();
429  TVector3 pMom=pPos-vPos;
430  pMom.SetMag(p->energy());
431 
432  if(p->distanceToTrack()<cuts->neutralCUT) continue;
433  if(p->nHitsEta()<cuts->etaHitsCUT) continue;
434  if(p->nHitsPhi()<cuts->phiHitsCUT) continue;
435  if(pPos.PseudoRapidity()<cuts->etaMinCUT) continue;
436  if(pPos.PseudoRapidity()>cuts->etaMaxCUT) continue;
437 
438  Float_t pT=pMom.Pt();
439  if(pMom.PseudoRapidity()>cuts->rapidityMinCUT && pMom.PseudoRapidity()<cuts->rapidityMaxCUT)
440  {
441  if(ev->refMult()>17)
442  {
443  if(t&1)
444  {
445  ptspec_mb_020->Fill(pT);
446  if(ratio>=cuts->ratioCUT)
447  ptspec_mb_bg_020->Fill(pT);
448  else
449  ptspec_mb_cut_020->Fill(pT);
450  }
451  if(t&2)
452  {
453  ptspec_ht1_020->Fill(pT);
454  if(ratio>=cuts->ratioCUT)
455  ptspec_ht1_bg_020->Fill(pT);
456  else
457  ptspec_ht1_cut_020->Fill(pT);
458  }
459  if(t&4)
460  {
461  ptspec_ht2_020->Fill(pT);
462  if(ratio>=cuts->ratioCUT)
463  ptspec_ht2_bg_020->Fill(pT);
464  else
465  ptspec_ht2_cut_020->Fill(pT);
466  }
467  if(t&8)
468  {
469  ptspec_sn_020->Fill(pT);
470  if(ratio>=cuts->ratioCUT)
471  ptspec_sn_bg_020->Fill(pT);
472  else
473  ptspec_sn_cut_020->Fill(pT);
474  }
475  }
476  else if(ev->refMult()>11)
477  {
478  if(t&1)
479  {
480  ptspec_mb_2040->Fill(pT);
481  if(ratio>=cuts->ratioCUT)
482  ptspec_mb_bg_2040->Fill(pT);
483  else
484  ptspec_mb_cut_2040->Fill(pT);
485  }
486  if(t&2)
487  {
488  ptspec_ht1_2040->Fill(pT);
489  if(ratio>=cuts->ratioCUT)
490  ptspec_ht1_bg_2040->Fill(pT);
491  else
492  ptspec_ht1_cut_2040->Fill(pT);
493  }
494  if(t&4)
495  {
496  ptspec_ht2_2040->Fill(pT);
497  if(ratio>=cuts->ratioCUT)
498  ptspec_ht2_bg_2040->Fill(pT);
499  else
500  ptspec_ht2_cut_2040->Fill(pT);
501  }
502  if(t&8)
503  {
504  ptspec_sn_2040->Fill(pT);
505  if(ratio>=cuts->ratioCUT)
506  ptspec_sn_bg_2040->Fill(pT);
507  else
508  ptspec_sn_cut_2040->Fill(pT);
509  }
510  }
511  else if(ev->refMult()>=0)
512  {
513  if(t&1)
514  {
515  ptspec_mb_40100->Fill(pT);
516  if(ratio>=cuts->ratioCUT)
517  ptspec_mb_bg_40100->Fill(pT);
518  else
519  ptspec_mb_cut_40100->Fill(pT);
520  }
521  if(t&2)
522  {
523  ptspec_ht1_40100->Fill(pT);
524  if(ratio>=cuts->ratioCUT)
525  ptspec_ht1_bg_40100->Fill(pT);
526  else
527  ptspec_ht1_cut_40100->Fill(pT);
528  }
529  if(t&4)
530  {
531  ptspec_ht2_40100->Fill(pT);
532  if(ratio<0.7) etaphi_fill3088b_ht2->Fill(pPos.PseudoRapidity(),pPos.Phi());
533  if(ratio>=cuts->ratioCUT)
534  {
535  ptspec_ht2_bg_40100->Fill(pT);
536  if(ev->runId()>=4060055&&ev->runId()<=4060088) etaphi_fill3085_ht2->Fill(pPos.PseudoRapidity(),pPos.Phi());
537  if(ev->runId()>=4061003&&ev->runId()<=4061025) etaphi_fill3088_ht2->Fill(pPos.PseudoRapidity(),pPos.Phi());
538  }
539  else
540  ptspec_ht2_cut_40100->Fill(pT);
541  }
542  if(t&8)
543  {
544  ptspec_sn_40100->Fill(pT);
545  if(ratio>=cuts->ratioCUT)
546  {
547  ptspec_sn_bg_40100->Fill(pT);
548  etaphi_sn->Fill(pPos.PseudoRapidity(),pPos.Phi());
549  if(ev->runId()>=4060055&&ev->runId()<=4060088) etaphi_fill3085_sn->Fill(pPos.PseudoRapidity(),pPos.Phi());
550  if(ev->runId()>=4061003&&ev->runId()<=4061025) etaphi_fill3088_sn->Fill(pPos.PseudoRapidity(),pPos.Phi());
551  }
552  else
553  ptspec_sn_cut_40100->Fill(pT);
554  }
555  }
556 
557  }
558 
559  }
560 
561  entry++;
562  }
563 
564  fout.close();
565 
566  ps->Close();
567 
568  TString oname(file);
569  oname.Prepend("outputBg113005_");
570  TFile o(oname.Data(),"RECREATE");
571 
572  etaphi_fill3085_ht2->Write();
573  etaphi_fill3088_ht2->Write();
574  etaphi_fill3085_sn->Write();
575  etaphi_fill3088_sn->Write();
576  etaphi_fill3088b_ht2->Write();
577  etaphi_sn->Write();
578 
579  hfrac_mb->Write();
580  hfrac_ht1->Write();
581  hfrac_ht2->Write();
582  hfrac_sn->Write();
583  h_trigger->Write();
584 
585  Ratio_mb_020->Write();
586  Ratio_ht1_020->Write();
587  Ratio_ht2_020->Write();
588  Ratio_sn_020->Write();
589  Ratio_mb_2040->Write();
590  Ratio_ht1_2040->Write();
591  Ratio_ht2_2040->Write();
592  Ratio_sn_2040->Write();
593  Ratio_mb_40100->Write();
594  Ratio_ht1_40100->Write();
595  Ratio_ht2_40100->Write();
596  Ratio_sn_40100->Write();
597 
598  EvsE_mb->Write();
599  EvsE_ht1->Write();
600  EvsE_ht2->Write();
601  EvsE_sn->Write();
602 
603  ptspec_mb_020->Write();
604  ptspec_mb_cut_020->Write();
605  ptspec_mb_bg_020->Write();
606  ptspec_ht1_020->Write();
607  ptspec_ht1_cut_020->Write();
608  ptspec_ht1_bg_020->Write();
609  ptspec_ht2_020->Write();
610  ptspec_ht2_cut_020->Write();
611  ptspec_ht2_bg_020->Write();
612  ptspec_sn_020->Write();
613  ptspec_sn_cut_020->Write();
614  ptspec_sn_bg_020->Write();
615 
616  ptspec_mb_2040->Write();
617  ptspec_mb_cut_2040->Write();
618  ptspec_mb_bg_2040->Write();
619  ptspec_ht1_2040->Write();
620  ptspec_ht1_cut_2040->Write();
621  ptspec_ht1_bg_2040->Write();
622  ptspec_ht2_2040->Write();
623  ptspec_ht2_cut_2040->Write();
624  ptspec_ht2_bg_2040->Write();
625  ptspec_sn_2040->Write();
626  ptspec_sn_cut_2040->Write();
627  ptspec_sn_bg_2040->Write();
628 
629  ptspec_mb_40100->Write();
630  ptspec_mb_cut_40100->Write();
631  ptspec_mb_bg_40100->Write();
632  ptspec_ht1_40100->Write();
633  ptspec_ht1_cut_40100->Write();
634  ptspec_ht1_bg_40100->Write();
635  ptspec_ht2_40100->Write();
636  ptspec_ht2_cut_40100->Write();
637  ptspec_ht2_bg_40100->Write();
638  ptspec_sn_40100->Write();
639  ptspec_sn_cut_40100->Write();
640  ptspec_sn_bg_40100->Write();
641 
642  tpc_mb->Write();
643  tpc_mb_cut->Write();
644  tpc_ht1->Write();
645  tpc_ht1_cut->Write();
646  tpc_ht2->Write();
647  tpc_ht2_cut->Write();
648  tpc_sn->Write();
649  tpc_sn_cut->Write();
650  ftpc_mb->Write();
651  ftpc_mb_cut->Write();
652  ftpc_ht1->Write();
653  ftpc_ht1_cut->Write();
654  ftpc_ht2->Write();
655  ftpc_ht2_cut->Write();
656  ftpc_sn->Write();
657  ftpc_sn_cut->Write();
658  bemc_mb->Write();
659  bemc_mb_cut->Write();
660  bemc_ht1->Write();
661  bemc_ht1_cut->Write();
662  bemc_ht2->Write();
663  bemc_ht2_cut->Write();
664  bemc_sn->Write();
665  bemc_sn_cut->Write();
666 
667  o.Close();
668 
669 }
670 
671 
672 void dAuBackground::runSim(const char *file)
673 {
674  TH1F *ratio_pythia=new TH1F("ratio_pythia","pythia, E over pT", 25,0.,1.);
675  TH1F *energy_pythia=new TH1F("energy_pythia","energy spectrum",200,0.,100.);
676  TH1F *tpcpt_pythia=new TH1F("tpcpt_pythia","TPC pT spectrum",200,0.,100.);
677  TH1F *spec_pythia=new TH1F("spec_pythia","pT spectrum neutral points",30,0.,15.);
678  TH1F *spec_pythia_bg=new TH1F("spec_pythia_bg","pT spectrum neutral points",30,0.,15.);
679 
680  TH1F *ratio_pythia_ht1=new TH1F("ratio_pythia_ht1","pythia, E over pT ht1", 25,0.,1.);
681  TH1F *energy_pythia_ht1=new TH1F("energy_pythia_ht1","energy spectrum ht1",200,0.,100.);
682  TH1F *tpcpt_pythia_ht1=new TH1F("tpcpt_pythia_ht1","TPC pT spectrum ht1",200,0.,100.);
683  TH1F *spec_pythia_ht1=new TH1F("spec_pythia_ht1","pT spectrum neutral points ht1",30,0.,15.);
684  TH1F *spec_pythia_bg_ht1=new TH1F("spec_pythia_bg_ht1","pT spectrum neutral points ht1",30,0.,15.);
685  TH1F *spec_pythia_ht1_cut=new TH1F("spec_pythia_ht1_cut","pT spectrum neutral points ht1 after cut",30,0.,15.);
686 
687  TH1F *pythia_tpc_ht1=new TH1F("pythia_tpc_ht1","tpc mult. ht-1 sim",100,0.,100.);
688  TH1F *pythia_tpc_ht1_cut=new TH1F("pythia_tpc_ht1_cut","tpc mult. ht-1 sim after cut",100,0.,100.);
689  TH1F *pythia_bemc_ht1=new TH1F("pythia_bemc_ht1","bemc mult. ht-1 sim",60,0.,60.);
690  TH1F *pythia_bemc_ht1_cut=new TH1F("pythia_bemc_ht1_cut","bemc mult. ht-1 sim after cut",60,0.,60.);
691 
692  TFile *inf=new TFile(file,"OPEN");
693  TTree *eventTree=(TTree*)inf->Get("mEventTree");
694  MyEvent *ev=new MyEvent();
695  eventTree->SetBranchAddress("branch",&ev);
696 
697  TPostScript *ps=new TPostScript("dAu_bgPythia.ps",111);
698  TCanvas *c9=new TCanvas("c9","c9",800,1000);
699  c9->Divide(2,4);
700 
701  Int_t entry=0;
702  while(eventTree->GetEntry(entry))
703  {
704  if(entry%10000==0) cout<<entry<<endl;
705 
706  Bool_t eventOK=kTRUE;
707  TVector3 vPos=ev->vertex();
708 
709  //apply EVENT cuts
710  if(TMath::Abs(vPos.Z())>cuts->vertexzCUT) eventOK=kFALSE;
711  if(TMath::Abs(vPos.X())>cuts->vertexxCUT) eventOK=kFALSE;
712  if(TMath::Abs(vPos.Y())>cuts->vertexyCUT) eventOK=kFALSE;
713  if(vPos.Mag()<.000001) eventOK=kFALSE;
714 
715  if(!eventOK)
716  {
717  entry++;
718  continue;
719  }
720 
721  Float_t en=ev->energyInBarrel();
722  Float_t pt=ev->momentumInTpc();
723  Float_t ratio=TMath::Abs(en+pt)>0 ? en/(en+pt) : 2.;
724 
725  ratio_pythia->Fill(ratio);
726  energy_pythia->Fill(en);
727  tpcpt_pythia->Fill(pt);
728  if(ev->highTowerAdc()>cuts->ht1AdcCUT)
729  {
730  pythia_tpc_ht1->Fill(ev->goodPrimaries());
731  pythia_bemc_ht1->Fill(ev->numberOfPoints());
732  ratio_pythia_ht1->Fill(ratio);
733  energy_pythia_ht1->Fill(en);
734  tpcpt_pythia_ht1->Fill(pt);
735  if(ratio<cuts->ratioCUT)
736  {
737  pythia_tpc_ht1_cut->Fill(ev->goodPrimaries());
738  pythia_bemc_ht1_cut->Fill(ev->numberOfPoints());
739  }
740  }
741 
742 
743  //get neutral points
744  MyPoint *p;
745  TClonesArray *clA=(TClonesArray*)ev->getPointArray();
746  for(Int_t j=0;j<clA->GetEntries();j++)
747  {
748  p=(MyPoint*)clA->At(j);
749  TVector3 pPos=p->position();
750  TVector3 pMom=pPos-vPos;
751  pMom.SetMag(p->energy());
752 
753  if(p->distanceToTrack()<cuts->neutralCUT) continue;
754  if(p->nHitsEta()<cuts->etaHitsCUT) continue;
755  if(p->nHitsPhi()<cuts->phiHitsCUT) continue;
756  if(pPos.PseudoRapidity()<cuts->etaMinCUT) continue;
757  if(pPos.PseudoRapidity()>cuts->etaMaxCUT) continue;
758 
759  Float_t pT=pMom.Pt();
760  if(pMom.PseudoRapidity()>cuts->rapidityMinCUT && pMom.PseudoRapidity()<cuts->rapidityMaxCUT)
761  {
762  spec_pythia->Fill(pT);
763  if(ratio>=cuts->ratioCUT) spec_pythia_bg->Fill(pT);
764  if(ev->highTowerAdc()>cuts->ht1AdcCUT)
765  {
766  spec_pythia_ht1->Fill(pT);
767  if(ratio>=cuts->ratioCUT) spec_pythia_bg_ht1->Fill(pT);
768  else spec_pythia_ht1_cut->Fill(pT);
769  }
770  }
771 
772 
773 
774  }
775 
776  entry++;
777  }
778 
779  TFile f("hijing_tryout.root","RECREATE");
780  spec_pythia_ht1->Write();
781  spec_pythia_ht1_cut->Write();
782  spec_pythia_bg_ht1->Write();
783  ratio_pythia_ht1->Write();
784 
785  pythia_tpc_ht1->Write();
786  pythia_tpc_ht1_cut->Write();
787  pythia_bemc_ht1->Write();
788  pythia_bemc_ht1_cut->Write();
789  f.Close();
790 
791  ps->NewPage();
792  c9->cd(1);
793  ratio_pythia->Draw();
794  c9->cd(2);
795  ratio_pythia_ht1->Draw();
796  c9->cd(3);
797  energy_pythia->Draw();
798  gPad->SetLogy();
799  c9->cd(4);
800  energy_pythia_ht1->Draw();
801  gPad->SetLogy();
802  c9->cd(5);
803  tpcpt_pythia->Draw();
804  gPad->SetLogy();
805  c9->cd(6);
806  tpcpt_pythia_ht1->Draw();
807  gPad->SetLogy();
808  c9->cd(7);
809  spec_pythia->Draw();
810  spec_pythia_bg->SetLineColor(6);
811  spec_pythia_bg->Draw("same");
812  gPad->SetLogy();
813  c9->cd(8);
814  spec_pythia_ht1->Draw();
815  spec_pythia_bg_ht1->SetLineColor(6);
816  spec_pythia_bg_ht1->Draw("same");
817  gPad->SetLogy();
818  c9->cd(0);
819  c9->Update();
820  ps->Close();
821 }
MyPoint
Definition: MyPoint.h:7
Generated by   doxygen 1.8.5