14 #include <TStopwatch.h>
15 #include <TObjArray.h>
16 #include <TObjString.h>
20 #include <TGraphErrors.h>
25 #include <TPaveText.h>
28 #include <AliRunLoader.h>
29 #include <AliHeader.h>
30 #include <AliGenPythiaEventHeader.h>
31 #include <AliJetParticle.h>
32 #include <AliJetParticlesReader.h>
33 #include <AliJetParticlesReaderKine.h>
34 #include <AliJetParticlesReaderKineGoodTPC.h>
35 #include <AliJetParticlesReaderESD.h>
36 #include <AliJetParticlesReaderHLT.h>
37 #include <AliJetEventParticles.h>
43 void Sort(Int_t na,const Int_t *a,Int_t *index);
44 void anaTracks(Int_t esd1,Int_t esd2,Float_t qcut,Int_t nMaxEvents=-1);
45 void anaTracks(Char_t *filename1,Char_t *filename2,Float_t qcut,Int_t nMaxEvents=-1);
46 void anaTracks(TObjArray *fnames1,TObjArray *fnames2,Float_t qcut,Int_t nMaxEvents=-1);
47 void anaEtTracks(Int_t esd1,Int_t esd2,Float_t qcut,
48 Int_t etfrom=1, Int_t etto=16, Int_t nMaxEvents=-1);
50 void anaTracks(Int_t esd1,Int_t esd2,Float_t qcut,Int_t nMaxEvents)
52 const Int_t nes[5]={50,100,150,200,250};
53 const Char_t *dir="/data/aliroot/jets";
54 TObjArray *f1=new TObjArray(1);
55 TObjArray *f2=new TObjArray(1);
56 Char_t filename1[1024];
57 Char_t filename2[1024];
59 for(Int_t i=0;i<5;i++){
60 sprintf(filename1,"%s/%d/aliev-type%d.root",dir,nes[i],esd1);
61 f1->Add(new TObjString(filename1));
62 sprintf(filename2,"%s/%d/aliev-type%d.root",dir,nes[i],esd2);
63 f2->Add(new TObjString(filename2));
65 anaTracks(f1,f2,qcut,nMaxEvents);
68 void anaEtTracks(Int_t esd1,Int_t esd2,Float_t qcut,Int_t etfrom, Int_t etto, Int_t nMaxEvents)
70 const Char_t *dir="/data/aliroot/scan";
71 TObjArray *f1=new TObjArray(1);
72 TObjArray *f2=new TObjArray(1);
73 Char_t filename1[1024];
74 Char_t filename2[1024];
75 for(Int_t i=etfrom;i<=etto;i++){
76 sprintf(filename1,"%s/et%d/aliev-type%d.root",dir,i,esd1);
77 f1->Add(new TObjString(filename1));
78 sprintf(filename2,"%s/et%d/aliev-type%d.root",dir,i,esd2);
79 f2->Add(new TObjString(filename2));
81 anaTracks(f1,f2,qcut,nMaxEvents);
84 void anaTracks(Char_t *filename1,Char_t *filename2,Float_t qcut,Int_t nMaxEvents)
86 filename1 = monte events
87 filename2 = tracked events
90 TObjArray *f1=new TObjArray(1);
91 f1->Add(new TObjString(filename1));
92 TObjArray *f2=new TObjArray(1);
93 f2->Add(new TObjString(filename2));
95 anaTracks(f1,f2,qcut,nMaxEvents);
98 void anaTracks(TObjArray *fnames1,TObjArray *fnames2,Float_t qcut,Int_t nMaxEvents)
100 filename1 = monte events
101 filename2 = tracked events
107 TChain *theTree1 = new TChain("AJEPtree");
108 for(Int_t i=0;i<fnames1->GetEntries();i++){
109 theTree1->Add(((TObjString*)fnames1->At(i))->String());
111 AliJetEventParticles *ev1=new AliJetEventParticles();
112 theTree1->SetBranchAddress("particles",&ev1);
115 TChain *theTree2 = new TChain("AJEPtree");
116 for(Int_t i=0;i<fnames2->GetEntries();i++){
117 theTree2->Add(((TObjString*)fnames2->At(i))->String());
119 AliJetEventParticles *ev2=new AliJetEventParticles();
120 theTree2->SetBranchAddress("particles",&ev2);
122 Int_t treeentries1=(Int_t)theTree1->GetEntries();
123 Int_t treeentries2=(Int_t)theTree2->GetEntries();
125 if((nMaxEvents<0) || (nMaxEvents>treeentries1))
126 nMaxEvents=treeentries1;
128 cout << "Found " << treeentries1 << ", but use " << nMaxEvents << " in chain1" << endl;
129 cout << "Found " << treeentries2 << ", but use " << nMaxEvents << " in chain2" << endl;
131 const Int_t nmax=25000;
132 Int_t *ind2=new Int_t[nmax];
133 const Float_t ptcut=0.5;
134 const Int_t nbinsx=20;
135 const Float_t binsx[nbinsx]={1,2,4,6,8,10,12.5,15,17.5,20,25,30,40,45,50,60,70,80,90,100};
138 TH1F *hInvQ2=new TH1F("hInvQ2","",1000,0,1.5);
139 TH1F *hInvQ2lab=new TH1F("hInvQ2lab","",1000,0,1.5);
140 TH1F *hEffPt1=new TH1F("hEffPt1","",nbinsx-1,binsx);
141 TH1F *hEffPt2=new TH1F("hEffPt2","",nbinsx-1,binsx);
142 TH1F *hEffEta1=new TH1F("hEffEta1","",40,-1,1);
143 TH1F *hEffEta2=new TH1F("hEffEta2","",40,-1,1);
144 TH1F *hEffPhi1=new TH1F("hEffPhi1","",62,0,6.2);
145 TH1F *hEffPhi2=new TH1F("hEffPhi2","",62,0,6.2);
146 TH1F *hFakePt=new TH1F("hFakePt","",nbinsx-1,binsx);
147 TH1F *hFakeEta=new TH1F("hFakeEta","",40,-1,1);
148 TH1F *hFakePhi=new TH1F("hFakePhi","",62,0,6.2);
149 TH1F *hResPt=new TH1F("hResPt","",60,-15,15);
150 TH1F *hResEta=new TH1F("hResEta","",200,-0.05,0.05);
151 TH1F *hResPhi=new TH1F("hResPhi","",200,-0.05,0.05);
152 TProfile *hpResPt=new TProfile("hpResPt","",nbinsx-1,binsx);
153 TProfile *hpResEta=new TProfile("hpResEta","",nbinsx-1,binsx);
154 TProfile *hpResPhi=new TProfile("hpResPhi","",nbinsx-1,binsx);
155 TH1F **hResPtAll=new TH1F*[nbinsx];
156 TH1F **hResEtaAll=new TH1F*[nbinsx];
157 TH1F **hResPhiAll=new TH1F*[nbinsx];
158 for(Int_t i=0;i<nbinsx;i++){
160 sprintf(fn,"hResPt%d",i);
161 hResPtAll[i]=new TH1F(fn,"",60,-15,15);
162 sprintf(fn,"hResEta%d",i);
163 hResEtaAll[i]=new TH1F(fn,"",200,-0.05,0.05);
164 sprintf(fn,"hResPhi%d",i);
165 hResPhiAll[i]=new TH1F(fn,"",200,-0.05,0.05);
167 //=========================================================================
168 // start the event loop
169 //=========================================================================
175 Int_t nlokfakes2 = 0;
177 while(nEvent1<nMaxEvents && nEvent2<nMaxEvents){
180 if((nEvent1 % 100) == 99) {
181 cout << "Reading event " << nEvent1 << " " << nlabs << endl;
187 theTree1->GetEvent(nEvent1);
188 theTree2->GetEvent(nEvent2);
190 if(ev1->GetEventNr()!=ev2->GetEventNr()){
191 cerr << "Need to skip event: " << nEvent1 << " != " << nEvent2
192 << " (" << ev1->GetEventNr() << " " << ev2->GetEventNr() <<") " << endl;
193 if(ev1->GetEventNr()>ev2->GetEventNr()) nEvent2++;
200 const TClonesArray *p1=ev1->GetParticles();
204 const TClonesArray *p2=ev2->GetParticles();
208 const Int_t n1=p1->GetEntriesFast();
209 const Int_t n2=p2->GetEntriesFast();
213 //clean found particle table
214 for(Int_t j=0;j<n2;j++) {
218 //Compute Q^2 with/without labels
219 for(Int_t i=0;i<n1;i++) {
220 AliJetParticle *p=(AliJetParticle*)p1->At(i);
221 if(p->Pt()<ptcut) continue;
226 Int_t lab=p->GetLabel();
227 if(lab==0) lab=p->GetUID(); //workaround for bug
228 for(Int_t j=0;j<n2;j++) {
229 AliJetParticle *q=(AliJetParticle*)p2->At(j);
230 if(q->Pt()<ptcut) continue;
231 if(q->GetNhts()<MINHITS) continue;
235 Float_t Q2=(px-qx)*(px-qx)+(py-qy)*(py-qy)+(pz-qz)*(pz-qz);
236 Q2=TMath::Sqrt(Q2)/p->P();
238 if(q->GetLabel()!=lab) continue;
243 for(Int_t i=0;i<n1;i++) {
244 AliJetParticle *p=(AliJetParticle*)p1->At(i);
245 if(p->Pt()<ptcut) continue;
249 Int_t lab=p->GetLabel();
250 if(lab==0) lab=p->GetUID(); //workaround for bug
254 for(Int_t j=0;j<n2;j++) {
255 if(ind2[j]!=1) continue; //particle was already found
256 AliJetParticle *q=(AliJetParticle*)p2->At(j);
257 if(q->Pt()<ptcut) continue;
258 if(q->GetNhts()<MINHITS) continue;
262 Float_t Q2=(px-qx)*(px-qx)+(py-qy)*(py-qy)+(pz-qz)*(pz-qz);
263 Q2=TMath::Sqrt(Q2)/p->P();
265 if(TMath::Abs(q->GetLabel())==lab){
270 } else { //take particle with smaller Q2
272 ind2[j2]=2;//mark old found as fake
276 } else ind2[j]=2;//mark new found as fake
280 if(Q2<minQ2){ //take particle with smaller Q2
288 hEffPt1->Fill(p->Pt());
289 hEffEta1->Fill(p->Eta());
290 hEffPhi1->Fill(p->Phi());
292 if(j2>=0){ //take found particle
299 if(minQ2>qcut) continue;
302 ind2[j2]=0; //mark as found
303 AliJetParticle *q=(AliJetParticle*)p2->At(j2);
304 hEffPt2->Fill(q->Pt());
305 hEffEta2->Fill(q->Eta());
306 hEffPhi2->Fill(q->Phi());
307 Float_t ptdiff=p->Pt()-q->Pt();
308 Float_t phidiff=p->Phi()-q->Phi();
309 Float_t etadiff=p->Eta()-q->Eta();
310 hResPt->Fill(ptdiff);
311 hResEta->Fill(etadiff);
312 hResPhi->Fill(phidiff);
313 //cout << ptdiff << etadiff << phidiff << endl;
314 hpResPt->Fill(p->Pt(),TMath::Abs(ptdiff)/p->Pt());
315 hpResEta->Fill(p->Pt(),TMath::Abs(etadiff));
316 hpResPhi->Fill(p->Pt(),TMath::Abs(phidiff));
317 Int_t index=hEffPt1->FindBin(p->Pt());
318 hResPtAll[index]->Fill(ptdiff);
319 hResEtaAll[index]->Fill(etadiff);
320 hResPhiAll[index]->Fill(phidiff);
325 for(Int_t j=0;j<n2;j++) {
326 if(!ind2[j]) continue;
327 AliJetParticle *q=(AliJetParticle*)p2->At(j);
328 if (q->Pt()<ptcut) continue;
329 //cout << ind2[j] << " " << q->GetLabel() << endl;
330 if(q->GetNhts()<MINHITS) continue;
331 hFakePt->Fill(q->Pt());
332 hFakeEta->Fill(q->Eta());
333 hFakePhi->Fill(q->Phi());
337 if(nlokfakes2/n2>0.1) cout << nEvent1 << " has more than 10% fakes: " << nlokfakes2 << endl;
340 Float_t ffakes=nfakes2*100./ntracks1;
341 cout << "Good tracks: " << ntracks1 << endl;
342 cout << "Found tracks: " << ntracks2 << endl;
343 cout << "Fake tracks: " << nfakes2 << endl;
344 cout << "Eff: " << ntracks2*100/ntracks1 << endl;
345 cout << "Fake: " << nfakes2*100/ntracks1 << endl;
346 TF1 *f1 = new TF1("f1","gaus");
348 Float_t m = hResPt->GetMean();
349 Float_t s = hResPt->GetRMS();
350 f1->SetParameter(1,m);
351 f1->SetParameter(2,s);
352 cout << "Res: " << m << " " << s << endl;
354 hResPt->Fit(f1,"IN");
355 m = f1->GetParameter(1);
356 s = f1->GetParameter(2);
357 cout << "Res: " << m << " " << s << endl;
359 TCanvas *c1=new TCanvas("cInvQ2","",600,400);
365 ival=hInvQ2->Integral(hInvQ2->FindBin(0),hInvQ2->FindBin(qq));
366 if(ival>ntracks1) break;
369 cout << "Qcut: " << qq << " " << ival << " " << ntracks1 << endl;
370 TCanvas *c1b=new TCanvas("cInvQ2lab","",600,400);
377 hEffPt2->Divide(hEffPt1);
381 hEffEta2->Divide(hEffEta1);
385 hEffPhi2->Divide(hEffPhi1);
388 hFakePt->Divide(hEffPt1);
390 hFakeEta->Divide(hEffEta1);
392 hFakePhi->Divide(hEffPhi1);
394 TGraphErrors *gpts=new TGraphErrors();
395 TGraphErrors *getas=new TGraphErrors();
396 TGraphErrors *gphis=new TGraphErrors();
399 for(Int_t i=0;i<nbinsx;i++){
400 hResPtAll[i]->Sumw2();
401 m = hResPtAll[i]->GetMean();
402 s = hResPtAll[i]->GetRMS();
404 f1->SetParameter(1,m);
405 f1->SetParameter(2,s);
406 hResPtAll[i]->Fit(f1,"IN");
407 m = f1->GetParameter(1);
408 s = f1->GetParameter(2);
409 //gpts->SetPointError(i,0,f1->GetParError(2));
411 gpts->SetPoint(i,binsx[i],s/binsx[i]*100);
413 hResEtaAll[i]->Sumw2();
414 m = hResEtaAll[i]->GetMean();
415 s = hResEtaAll[i]->GetRMS();
417 f1->SetParameter(1,m);
418 f1->SetParameter(2,s);
419 hResEtaAll[i]->Fit(f1,"QIN");
420 m = f1->GetParameter(1);
421 s = f1->GetParameter(2);
423 getas->SetPoint(i,binsx[i],s);
425 hResPhiAll[i]->Sumw2();
427 m = hResPhiAll[i]->GetMean();
428 s = hResPhiAll[i]->GetRMS();
429 f1->SetParameter(1,m);
430 f1->SetParameter(2,s);
431 hResPhiAll[i]->Fit(f1,"QIN");
432 m = f1->GetParameter(1);
433 s = f1->GetParameter(2);
435 gphis->SetPoint(i,binsx[i],s);
439 TCanvas *c2=new TCanvas("cEffPt","",600,400);
442 TCanvas *c3=new TCanvas("cEffEta","",600,400);
445 TCanvas *c4=new TCanvas("cEffPhi","",600,400);
449 TCanvas *c5=new TCanvas("cResPt","",600,400);
451 gpts->SetMarkerSize(1);
452 gpts->SetMarkerStyle(20);
453 //gpts->Draw("p axis");
455 TCanvas *c6=new TCanvas("cResPhi","",600,400);
457 gphis->SetMarkerSize(1);
458 gphis->SetMarkerStyle(20);
459 //gphis->Draw("p axis");
461 TCanvas *c7=new TCanvas("cResEta","",600,400);
463 getas->SetMarkerSize(1);
464 getas->SetMarkerStyle(20);
465 //getas->Draw("p axis");
468 TCanvas *c8=new TCanvas("call","",1200,800);
472 hFakePt->Draw("same");
475 hFakeEta->Draw("same");
478 hFakePhi->Draw("same");
487 TFile *f=new TFile("tracks.root","recreate");
503 getas->Write("gtps");
504 gphis->Write("gtps");
507 for(Int_t i=0;i<nbinsx;i++){
508 hResPtAll[i]->Write();
509 hResEtaAll[i]->Write();
510 hResPhiAll[i]->Write();
521 //_____________________________________________________________________________
524 void Sort(Int_t na, const Int_t *a, Int_t *index)
527 Int_t *localArr1 = new Int_t[na];
528 Int_t *localArr2 = new Int_t[na];
532 for(iEl = 0; iEl < na; iEl++) {
533 localArr1[iEl] = a[iEl];
534 localArr2[iEl] = iEl;
535 //cout << localArr1[iEl] << " " << localArr2[iEl] << endl;
538 for (iEl = 0; iEl < na; iEl++) {
539 for (iEl2 = na-1; iEl2 > iEl; iEl2--) {
540 //cout << iEl << " " << iEl2 << " " << localArr1[iEl2-1] << " " << localArr1[iEl2] << endl;
542 if (localArr1[iEl2-1] < localArr1[iEl2]) {
543 Int_t tmp = localArr1[iEl2-1];
544 localArr1[iEl2-1] = localArr1[iEl2];
545 localArr1[iEl2] = tmp;
547 Int_t tmp2 = localArr2[iEl2-1];
548 localArr2[iEl2-1] = localArr2[iEl2];
549 localArr2[iEl2] = tmp2;
550 //cout << "Swapped: " << localArr1[iEl2-1] << " " << localArr2[iEl2-1] << endl;
551 //cout << "Swapped: " << localArr1[iEl2] << " " << localArr2[iEl2] << endl;
553 //cout << iEl << " " << iEl2 << " " << localArr1[iEl2-1] << " " << localArr1[iEl2] << endl;
557 for (iEl = 0; iEl < na; iEl++) {
558 index[iEl] = localArr2[iEl];
559 //cout << a[iEl] << " " << index[iEl] << endl;