L3 becomes HLT
[u/mrichter/AliRoot.git] / HLT / exa / trigger_pp.C
1 // $Id$
2
3 /**
4    trigger_pp -> Calculate efficiency of trigger and fake tracks in pp pileup events
5    plot_pp -> Plot vertex resolution
6    plot_pp_vert -> Plot vertex resolution
7    If you are looking for the old trigger macro search for 
8    OLD_TRIGGER_MACRO.
9
10    Author: C. Loizides
11 */
12
13
14 /* ---trigger_pp---
15    Calculate the efficiency of found tracks in pileup event 
16    and found fake tracks after applying different DCA cuts to 
17    the assumed vertex of (0,0,0). Tracks are stored in the ntuple
18    root file (see fill_pp.C). Good tracks can either be used from
19    the single pp event or the good_tracks_tpc event (the macro
20    will detect which ntuple file you specified).
21
22    Script: Track single pp events and store in ntuple:
23 --------------------------------------------------------
24 #!/bin/bash
25
26 dir=/mnt/data/loizides/alidata/head
27 wdir=$dir/pp/getrackted-0-375/
28 mdir=~/work/l3/level3code/exa/
29 sfile=~/pp-tracks.root
30 rfile=$dir/PythiaPP_1000Events_TPConly_ev0-375_digits.root
31
32 cd $dir
33 ln -s -f $rfile digitfile.root
34 cd $mdir;
35
36 for i in `seq 0 365`; do
37  echo $i
38  trigev=$i
39  cd $mdir
40  aliroot -b -q runtracker_pp.C\($trigev,0,\"$rfile\",\"$wdir\"\)
41  aliroot -b -q fill_pp.C\($trigev,\"$wdir\",\"$sfile\"\)
42 done
43 --------------------------------------------------------
44
45    Script: Track pileup events and store in ntuple:
46 -----------------------------------------------------
47 #!/bin/bash
48
49 dir=/mnt/data/loizides/alidata/head
50 wdir=/mnt/data/loizides/alidata/head/pp/pileup-25
51 mdir=~/l3/level3code/exa
52 sfile=~/pileup25-tracks.root
53
54 cd $mdir;
55
56 for i in `seq 0 100`; do
57  echo $i
58  pdir=$wdir/pileup-$i/
59  cd $pdir
60  trigev=`ls digits_*_0_-1.raw | cut -d_ -f 2`
61  cd $mdir
62  aliroot -b -q runtracker_pp.C\($trigev,\"$pdir\",0,\"$pdir\"\)
63  aliroot -b -q fill_pp.C\($trigev,\"$pdir\",\"$sfile\",$i\)
64 done
65 -----------------------------------------------------
66
67   Parameter to trigger_pp are:
68   Pileupfile is the ntuple root file containing the information of the pileup event
69   Ppfile gives the ntuple root file containing the reference information of the single pp event (can either be good_particles or single tracked pp events)
70   Evi and Eve give start and end event number to process.
71 */
72
73 void trigger_pp(Char_t *pileupfile, Char_t *ppfile, Int_t evi=0, Int_t eve=100)
74 {
75
76   //-----------------------------------------------------------
77   const Char_t *LTEXT_="for 25 piles, min. 2 trigger tracks";  
78   const Int_t MINTRACKS_=2; //min. tracks of the triggered 
79   const Int_t N_=100;       //number of data points
80   const Float_t WIDTH_=1; //the width of the xbin
81   //-----------------------------------------------------------
82
83
84   //"pt:phi:eta:xvert:yvert:zvert:imp:nhits:px:py:pz:event:mc")
85   Float_t pt1,phi1,eta1,xvert1,yvert1,zvert1,imp1,nhits1,px1,py1,pz1,event1,mc1; //pileup
86   Float_t pt2,phi2,eta2,xvert2,yvert2,zvert2,imp2,nhits2,px2,py2,pz2,event2,mc2; //pp
87
88   Int_t countgoodevents=0;
89   Int_t nent1,nent2;
90   Char_t dummy[1024];
91
92   Int_t tnorm;
93   Int_t fnorm;
94   Int_t tcounts[N_];
95   Int_t fcounts[N_];
96
97   Float_t tallcounts[N_];
98   Float_t fallcounts[N_];
99   Int_t totnorm=0;
100
101   Float_t thcounts[N_];
102   Float_t fhcounts[N_];
103   Float_t thmeancounts[N_];
104   Float_t fhmeancounts[N_];
105   Float_t thsigmacounts[N_];
106   Float_t fhsigmacounts[N_];
107   Float_t xcut[N_];
108   for(Int_t k=0;k<N_;k++){
109     tallcounts[k]=0;
110     fallcounts[k]=0;
111     thmeancounts[k]=0;
112     fhmeancounts[k]=0;
113     thsigmacounts[k]=0;
114     fhsigmacounts[k]=0;
115     xcut[k]=(k+1)*WIDTH_; 
116   }
117
118   TFile file1=TFile(pileupfile,"READ");
119   TFile file2=TFile(ppfile,"READ");
120   if(!file1.IsOpen() || !file2.IsOpen()) return;
121
122   for(Int_t i=evi;i<eve;i++){ //loop over pileup-events
123
124     sprintf(dummy,"good-pptracks-%d",i);
125     TNtuple *ntuppel1 = (TNtuple*)file1.Get(dummy);
126     Float_t nhits1=0,pdg1=0,sector1=0;
127     if(ntuppel1){
128       ntuppel1->SetBranchAddress("mc",&mc1);
129       ntuppel1->SetBranchAddress("pdg",&pdg1);
130       ntuppel1->SetBranchAddress("px",&px1);
131       ntuppel1->SetBranchAddress("py",&py1);
132       ntuppel1->SetBranchAddress("pz",&pz1);
133       ntuppel1->SetBranchAddress("xvert",&xvert1);
134       ntuppel1->SetBranchAddress("yvert",&yvert1);
135       ntuppel1->SetBranchAddress("zvert",&zvert1);
136       ntuppel1->SetBranchAddress("nhits",&nhits1);
137       ntuppel1->SetBranchAddress("sector",&sector1);
138       event1=i;
139    } else {
140       sprintf(dummy,"pptracks-%d",i);
141       TNtuple *ntuppel1 = (TNtuple*)file1.Get(dummy);
142       if(!ntuppel1) continue;
143       ntuppel1->SetBranchAddress("pt",&pt1);
144       ntuppel1->SetBranchAddress("phi",&phi1);
145       ntuppel1->SetBranchAddress("eta",&eta1);
146       ntuppel1->SetBranchAddress("xvert",&xvert1);
147       ntuppel1->SetBranchAddress("yvert",&yvert1);
148       ntuppel1->SetBranchAddress("zvert",&zvert1);
149       ntuppel1->SetBranchAddress("imp",&imp1);
150       ntuppel1->SetBranchAddress("px",&px1);
151       ntuppel1->SetBranchAddress("py",&py1);
152       ntuppel1->SetBranchAddress("pz",&pz1);
153       ntuppel1->SetBranchAddress("event",&event1);
154       ntuppel1->SetBranchAddress("mc",&mc1);
155     }
156
157     tnorm=0;
158     fnorm=0;
159     for(Int_t j=0;j<N_;j++){
160       tcounts[j]=0;
161       fcounts[j]=0;
162       thcounts[j]=0;
163       fhcounts[j]=0;
164     }
165
166     nent1=ntuppel1->GetEntries();
167     for(Int_t j=0;j<nent1;j++){
168       ntuppel1->GetEvent(j);
169       if(nhits1!=0){ // file1 is good-particle file
170         pt1=TMath::Sqrt(px1*px1+py1*py1);
171         Float_t p=TMath::Sqrt(px1*px1+py1*py1+pz1*pz1);
172         eta1=0.5*TMath::log((p+pz1)/(p-pz1));
173         phi1=TMath::atan(py1/px1);
174         imp1=TMath::Sqrt(xvert1*xvert1+yvert1*yvert1+zvert1*zvert1);
175       }
176       if(j==0){ //get info of triggered event from file2
177         sprintf(dummy,"good-pptracks-%d",event1);
178         TNtuple *ntuppel2 = (TNtuple*)file2.Get(dummy);
179         if(ntuppel2){
180           Float_t nhits2,pdg2,sector2;
181           ntuppel2->SetBranchAddress("mc",&mc2);
182           ntuppel2->SetBranchAddress("pdg",&pdg2);
183           ntuppel2->SetBranchAddress("px",&px2);
184           ntuppel2->SetBranchAddress("py",&py2);
185           ntuppel2->SetBranchAddress("pz",&pz2);
186           ntuppel2->SetBranchAddress("xvert",&xvert2);
187           ntuppel2->SetBranchAddress("yvert",&yvert2);
188           ntuppel2->SetBranchAddress("zvert",&zvert2);
189           ntuppel2->SetBranchAddress("nhits",&nhits2);
190           ntuppel2->SetBranchAddress("sector",&sector2);
191           event2=event1;
192           nent2=ntuppel2->GetEntries();
193           for(Int_t k=0;k<nent2;k++){
194             ntuppel2->GetEvent(k);
195             if(mc2<0)continue; //dont count fake tracks
196             //if(TMath::fabs(zvert2)>30) continue;
197             pt2=TMath::Sqrt(px2*px2+py2*py2);
198             if(nhits2<63) continue;
199             if(pt2<0.1) continue;
200             if((px2==0)&&(py2==0)&&(pz2==0)) continue;
201             Float_t p=TMath::Sqrt(px2*px2+py2*py2+pz2*pz2);
202             eta2=0.5*TMath::log((p+pz2)/(p-pz2));
203             if(TMath::fabs(eta2)>0.9) continue;
204             phi2=TMath::atan(py2/px2);
205             imp2=TMath::Sqrt(xvert2*xvert2+yvert2*yvert2+zvert2*zvert2);
206             tnorm++;
207           }
208         } else {
209           sprintf(dummy,"pptracks-%d",event1);
210           TNtuple *ntuppel2 = (TNtuple*)file2.Get(dummy);
211           if(ntuppel2){
212             ntuppel2->SetBranchAddress("pt",&pt2);
213             ntuppel2->SetBranchAddress("phi",&phi2);
214             ntuppel2->SetBranchAddress("eta",&eta2);
215             ntuppel2->SetBranchAddress("xvert",&xvert2);
216             ntuppel2->SetBranchAddress("yvert",&yvert2);
217             ntuppel2->SetBranchAddress("zvert",&zvert2);
218             ntuppel2->SetBranchAddress("imp",&imp2);
219             ntuppel2->SetBranchAddress("px",&px2);
220             ntuppel2->SetBranchAddress("py",&py2);
221             ntuppel2->SetBranchAddress("pz",&pz2);
222             ntuppel2->SetBranchAddress("event",&event2);
223             ntuppel2->SetBranchAddress("mc",&mc2);
224
225             nent2=ntuppel2->GetEntries();
226             for(Int_t k=0;k<nent2;k++){
227               ntuppel2->GetEvent(k);
228               if(mc2<0)continue; //dont count fake tracks
229               if(pt2<0.1) continue;
230               //if(TMath::fabs(zvert2)>30) continue;
231               if(TMath::fabs(eta2)>0.9) continue;
232               tnorm++;
233               //cout << k << ": " << pt2 << " " << mc2 <<endl;
234             }
235           }
236         }
237       }
238       if(tnorm<MINTRACKS_){
239         tnorm=0;
240         break; //triggered event has to have one track at least
241       }
242
243       if(pt1<0.1) continue;
244       if(TMath::fabs(eta1>0.9)) continue;
245
246       if(mc1<-1) continue; //dont count fake tracks
247       else if(mc1==-1) fnorm++;
248
249       for(Int_t k=0;k<N_;k++){ //do the counting
250         Float_t cut=xcut[k];
251         //Float_t impt=TMath::Sqrt(xvert1*xvert1+yvert1*yvert1);
252         Float_t impz=TMath::abs(zvert1);
253         if(impz>cut) continue;
254         if(mc1==-1) fcounts[k]++;
255         else tcounts[k]++;
256       }
257
258     }//pp track loop
259
260     if(tnorm==0) continue; //break from above
261     countgoodevents++; //for normalization
262
263     totnorm+=tnorm;
264
265     for(Int_t k=0;k<N_;k++){ //do the counting
266       if(tcounts[k]> tnorm)tcounts[k]=tnorm; 
267       thcounts[k]=(Float_t)tcounts[k]/tnorm;
268       fhcounts[k]=(Float_t)fcounts[k]/tnorm;
269       tallcounts[k]+=tcounts[k];
270       fallcounts[k]+=fcounts[k];
271       thmeancounts[k]+=thcounts[k];
272       fhmeancounts[k]+=fhcounts[k];
273       thsigmacounts[k]+=thcounts[k]*thcounts[k];
274       fhsigmacounts[k]+=fhcounts[k]*fhcounts[k];
275       if((k==N_-1)&&(thcounts[k]< 0.99)){ 
276         //cout << "Warning " << i << " " << event2 << " " << tcounts[k] << " " << tnorm << endl;
277       }
278     }
279   }//pileup loop
280
281   file1.Close();
282   file2.Close(); 
283
284   cout << "Events used: " << countgoodevents << endl;
285   for(Int_t k=0;k<N_;k++){ //do the counting
286     thmeancounts[k]/=countgoodevents;
287     fhmeancounts[k]/=countgoodevents;
288
289     tallcounts[k]/=totnorm;
290     fallcounts[k]/=totnorm;
291
292     if(countgoodevents>1){
293       Float_t stemp=thsigmacounts[k]/countgoodevents;
294       Float_t s2temp=fhsigmacounts[k]/countgoodevents;
295       Int_t N=countgoodevents-1;
296       thsigmacounts[k]=TMath::sqrt((stemp -thmeancounts[k]*thmeancounts[k])/N);
297       fhsigmacounts[k]=TMath::sqrt((s2temp-fhmeancounts[k]*fhmeancounts[k])/N);
298     }else{
299       thsigmacounts[k]=0;
300       fhsigmacounts[k]=0;
301     }
302     cout << k << ": " << thmeancounts[k] << "+-"<< thsigmacounts[k] << " " << fhmeancounts[k] << "+-" << fhsigmacounts[k] << endl;
303   }
304
305   TCanvas *c1 = new TCanvas("c1","PileUp",1000,500);
306   //c1->SetFillColor(42);
307   //c1->SetGrid();
308   //c1->GetFrame()->SetFillColor(21);
309   //c1->GetFrame()->SetBorderSize(12);
310
311   TGraphErrors *g1=new TGraphErrors(N_,xcut,thmeancounts,0,&thsigmacounts[0]);
312   //TGraphErrors *g1=new TGraphErrors(N_,xcut,tallcounts);
313   sprintf(dummy,"Good tracks %s",LTEXT_);
314   g1->SetTitle(dummy);
315   g1->SetMarkerColor(4);
316   g1->GetHistogram()->SetXTitle("Z_{cut} [cm]");
317   //g1->GetHistogram()->SetYTitle("Efficiency");
318   g1->SetMarkerStyle(21);
319   
320   TGraphErrors *g2=new TGraphErrors(N_,xcut,fhmeancounts,0,fhsigmacounts);
321   //TGraphErrors *g2=new TGraphErrors(N_,xcut,fallcounts);
322   g2->SetTitle("Fake tracks (relative to good triggered tracks)");
323   g2->GetHistogram()->SetXTitle("Z_{cut} [cm]");
324   g2->SetMarkerColor(4);
325   g2->SetMarkerStyle(21);
326
327   c1->Divide(2,1);
328   c1->cd(1);
329   g1->Draw("AP");
330   c1->cd(2);
331   g2->Draw("AP");
332   c1->Update();
333 }
334
335 /* Plot the vertex reconstruction resolution */
336
337 void plot_pp_vert(Int_t evi=0,Int_t evs=100,Char_t *infile="pp-tracks.root")
338 {
339   const Int_t zcut=3;
340   const Int_t zcutplot=5;
341
342   Char_t dummy[1000];
343   Char_t fname[1000];
344   Char_t fname2[1000];
345
346   Float_t pt,phi,eta,xvert,yvert,zvert,imp,nhits,px,py,pz,event,mc;
347   TH1F *vhist = new TH1F("zhist","Vertex reconstruction resolution",50,-zcutplot,zcutplot);
348
349   TFile file = TFile(infile,"READ");
350   if(!file.IsOpen()) return;
351
352   TNtuple *ntuppel;
353
354   for(int ev=evi; ev<evs; ev++){
355     //sprintf(dummy,"good-pptracks-%d",ev); /*dont use this*/
356     sprintf(dummy,"pptracks-%d",ev);
357
358     if(ntuppel) delete ntuppel;
359     ntuppel = (TNtuple*)file.Get(dummy);
360     if(!ntuppel) continue;
361     if(ntuppel.GetEntries()==0) continue;
362
363     sprintf(fname,"(eta >-0.9) && (eta<0.9) && (pt>0.1) && (nhits>63) && (zvert<%d) && (zvert>-%d)",zcut,zcut);
364     ntuppel->Draw("zvert>>histo",fname,"groff");
365     Float_t mean = histo->GetMean();
366     //cout << mean << endl;
367     vhist->Fill(mean);
368   }
369   
370   gStyle->SetStatColor(10);
371   gStyle->SetOptFit(1);
372
373   TF1 *f1 = new TF1("f1","gaus",-(Float_t)zcutplot/2.,(Float_t)zcutplot/2.);
374   vhist->Fit("f1","R");
375
376   vhist->SetXTitle("Z* [cm]");
377   vhist->Draw("E1P");
378 }
379
380 /* Plot the impact parameter resolution */
381  
382 void plot_pp(Int_t evi=0,Int_t evs=100,Char_t *infile="pp-tracks.root")
383 {
384   const Int_t zcut=3;
385   const Int_t zcutplot=5;
386
387   Char_t dummy[1000];
388   Char_t fname[1000];
389   Char_t fname2[1000];
390
391   TFile file = TFile(infile,"READ");
392   if(!file.IsOpen()) return;
393   
394   TH1F *vhist = new TH1F("zhist","Impact parameter resolution",50,-zcut,zcut);
395   
396   TNtuple *ntuppel;
397   Int_t a=0;
398   for(int ev=evi; ev<evs; ev++){
399     //sprintf(dummy,"good-pptracks-%d",ev); /*dont use this*/
400     sprintf(dummy,"pptracks-%d",ev);
401
402     if(ntuppel) delete ntuppel;
403     ntuppel = (TNtuple*)file.Get(dummy);
404     if(!ntuppel) continue;
405     if(ntuppel.GetEntries()==0) continue;
406
407     sprintf(fname,"(eta>-0.9) && (eta<0.9) && (pt>0.1) && (nhits>63) && (zvert<%d) && (zvert>-%d)",zcut,zcut);
408     ntuppel->Draw("zvert>>histo",fname,"groff");
409     Float_t mean = histo->GetMean();
410     a++;
411     //cout << mean << endl;
412
413     sprintf(fname2,"(zvert-%f)>>+zhist",mean);
414     ntuppel->Draw(fname2,fname);
415   }
416   
417   gStyle->SetStatColor(10);
418   gStyle->SetOptFit(1);
419   cout << a <<endl;
420   TF1 *f1 = new TF1("f1","gaus",-(Float_t)zcutplot/2.,(Float_t)zcutplot/2.);
421   vhist->Fit("f1","R");
422
423   vhist->SetXTitle("Z_{DCA}-Z* [cm]");
424   vhist->DrawCopy("E1P");
425   file->Close();
426 }
427
428
429 #ifdef OLD_TRIGGER_MACRO
430 void trigger_pp(char *outfile="results.root")
431 {
432   
433   TNtuple *ntuppel = new TNtuple("ntuppel","","pt:eta:xvert:yvert:zvert:nhits:px:py:pz:event");
434   Float_t meas[10];  
435   
436   for(int event=0; event<1; event++)
437     {
438       char fname[256];
439       sprintf(fname,"/data1/AliRoot/pp/pileup/tracks_0.raw");
440       //sprintf(fname,"aliruntfile.root");
441
442       //Get the tracks:
443       /*AliHLTTrackArray *tracks = new AliHLTTrackArray();
444       AliHLTFileHandler *file = new AliHLTFileHandler();
445       file->SetBinaryInput(fname);
446       file->Binary2TrackArray(tracks);
447       file->CloseBinaryInput();
448       delete file;*/
449
450       sprintf(fname,"/data1/AliRoot/pp/pileup/");
451         
452       AliHLTEvaluate *eval=new AliHLTEvaluate(fname,63,63);
453       eval->LoadData(event,-1);
454       eval->AssignIDs();
455       
456       AliHLTTrackArray *tracks=eval->GetTracks();
457
458       sprintf(fname,"/data1/AliRoot/pp/pileup/");
459       //sprintf(fname,"/prog/alice/data/Rawdata/1_patch/pp/recon_%d/",event);
460
461       Int_t ntracks=0;
462       Double_t xc,yc,zc;
463       Double_t impact;
464       AliHLTVertex vertex;
465
466       Int_t mcid = 0;
467
468       AliHLTTrackArray *ftracks = new AliHLTTrackArray();
469       
470       for(int i=0; i<tracks->GetNTracks(); i++)
471         {
472           track = (AliHLTTrack*)tracks->GetCheckedTrack(i);
473           if(!track) continue;
474           
475           //Assign MCid
476           mcid = track->GetMCid();
477           //if (mcid != -1) 
478           //cout << "MCid " << mcid << endl; 
479           
480           track->CalculateHelix();
481           //cout<<"Pt "<<track->GetPt()<<" eta "<<track->GetPseudoRapidity()<<" Nhits "<<track->GetNHits()<<endl;
482           //cout<<"Before second fit; pt "<<track->GetPt()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
483
484           track->CalculateHelix();
485           track->GetClosestPoint(&vertex,xc,yc,zc);
486           meas[0]=track->GetPt();
487           meas[1]=track->GetPseudoRapidity();
488           meas[2]=xc;
489           meas[3]=yc;
490           meas[4]=zc;
491           meas[5]=track->GetNHits();
492           meas[6]=track->GetPx();
493           meas[7]=track->GetPy();
494           meas[8]=track->GetPz();
495           meas[9]=event;
496           ntuppel->Fill(meas);
497           //continue;
498           if(fabs(track->GetPseudoRapidity())>0.9) continue;
499           if(track->GetNHits() < 100) continue;
500           if(track->GetPt()<0.2) continue;
501
502           impact = sqrt(xc*xc+yc*yc+zc*zc);
503           if(fabs(zc)>3) continue;
504           ftracks->AddLast(track);
505           cout<<"-------------------------------------"<<endl;
506           cout<<"Number of hits "<<track->GetNHits()<<endl;
507           cout<<"Transversal impact "<<sqrt(xc*xc+yc*yc)<<endl;
508           cout<<"Longitudinal impact "<<zc<<endl;
509           cout<<"Total impact "<<sqrt(xc*xc+yc*yc+zc*zc)<<endl;
510           cout<<"xc "<<xc<<" yc "<<yc<<" zc "<<zc<<endl;
511           
512           //cout<<"After second fit; pt "<<track->GetPt()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
513           cout<<"pt "<<track->GetPt()<<" eta "<<track->GetPseudoRapidity()<<endl;
514           
515           ntracks++;
516         }
517       cout<<endl<<"There was "<<ntracks<<" accepted tracks, out of total "<<tracks->GetNTracks()<<" found tracks"<<endl;
518       
519       //display(ftracks,fname);
520       delete tracks;
521       delete ftracks;
522     }
523   
524   TFile *rfile = TFile::Open(outfile,"RECREATE");
525   ntuppel->Write();
526   rfile->Close();
527   delete ntuppel;
528 }
529
530 void display(AliHLTTrackArray *tracks,char *path)
531 {
532   int slice[2]={0,35};
533   d = new AliHLTDisplay(slice);
534   d->Setup("tracks_0.raw",path);
535   d->SetTracks(tracks);
536   //d->DisplayClusters();
537   d->DisplayAll();
538   //d->DisplayTracks();
539   
540 }
541
542 void ploteff()
543 {
544   //double x[6]={1,2,4,6,8,10};
545   //double y[6]={0.873684,1.01379,1.17751,1.28614,1.31638,1.32022};
546   
547   double x[6]={0.8,1.6,3.2,4.8,6.4,8};
548   double y[6]={0.497006,0.88024,1.19162,1.23952,1.39521,1.40719};//->sigmas
549
550   
551   TGraph *gr = new TGraph(6,x,y);
552   c1 = new TCanvas("c1","",2);
553   SetCanvasOptions(c1);
554   gr->SetMarkerStyle(20);
555   gr->SetTitle("");
556   gr->Draw("APL");
557   gr->GetHistogram()->SetXTitle("Z_{CUT} [x #sigma_{Z}]");
558   gr->GetHistogram()->SetYTitle("Efficiency");
559   SetTGraphOptions(gr);
560 }
561
562 double geteff(char *infile,char *singlefile,double cut)
563 {
564   gStyle->SetStatColor(10);
565   gStyle->SetOptFit(1);
566   Int_t pileups[25],single[25];
567   file = TFile::Open(infile,"READ");
568   
569   char name[256];
570   for(int i=0; i<25; i++)
571     {
572       sprintf(name,"zvert < %f && zvert > %f && pt>0.2 && nhits>100 && eta>-0.9 && eta<0.9 && event==%d",cut,-1.*cut,i);
573       ntuppel->Draw(">>eventlist",name);
574       int entries = eventlist->GetN();
575       pileups[i]=entries;
576     }
577   //eventlist->Print("all");
578   file->Close();
579   
580   file = TFile::Open(singlefile,"read");
581   for(int i=0; i<25; i++)
582     {
583       sprintf(name,"zvert < %f && zvert > %f && pt>0.2 && nhits>100 && eta>-0.9 && eta<0.9 && event==%d",3,-3,i);
584       ntuppel->Draw(">>eventlist",name);
585       int entries = eventlist->GetN();
586       single[i]=entries;
587     }
588   file->Close();
589   double totsingle=0,totpileup=0;
590   for(int i=0; i<25; i++)
591     {
592       totsingle += single[i];
593       totpileup += pileups[i];
594     }
595   double toteff = totpileup/totsingle;
596   cout<<"Total eff "<<toteff<<endl;
597   
598   return toteff;
599 }
600
601
602 void plot(char *infile)
603 {
604   gStyle->SetStatColor(10);
605   gStyle->SetOptFit(1);
606   file = TFile::Open(infile,"READ");
607   
608   histo = new TH1F("histo","histo",20,-3,3);
609   
610   vhist = new TH1F("vhist","",20,-3,3);
611   SetTH1Options(vhist);
612   
613   char fname[256];
614   char fname2[256];
615   for(int ev=0; ev<25; ev++)
616     {
617       sprintf(fname,"pt>0.2 && nhits>100 && eta>-0.9 && eta<0.9 && event==%d",ev);
618       ntuppel->Draw("zvert>>histo",fname,"goff");
619       float mean = histo->GetMean();
620       vhist->Fill(mean);
621       continue;
622       sprintf(fname2,"(zvert-(%f))>>+vhist",mean);
623       cout<<fname2<<endl;
624       ntuppel->Draw(fname2,fname,"goff");
625     }
626   
627   c1 = new TCanvas("c1","",2);
628   SetCanvasOptions(c1);
629   vhist->SetXTitle("Z* [cm]");
630   vhist->Draw();
631   return;
632   //ntuppel->Draw("zvert>>histo","pt>0.2");
633   TF1 *f1 = new TF1("f1","gaus",-3,3);
634   histo->Fit("f1","R");
635
636   //histo->Draw();
637   //file->Close();
638 }
639
640 enum tagprimary {kPrimaryCharged=0x4000};
641 void LoadEvent(Int_t event=0)
642 {
643   //Load the generated particles
644   
645   gROOT->LoadMacro("$(ALICE_ROOT)/macros/loadlibs.C");
646   loadlibs();
647   //TFile *rootfile = TFile::Open("/prog/alice/data/pro/25event_pp.root","READ");
648   TFile *rootfile = TFile::Open("/prog/alice/pp/pileup/100pileup/alirunfile.root","READ");
649   gAlice = (AliRun*)rootfile->Get("gAlice");
650   
651   //  TNtuple *ntup = new TNtuple(
652   
653   Int_t nparticles = gAlice->GetEvent(event);
654   Int_t nprim = FindPrimaries(nparticles,0.,0.,0.);
655   cout<<"Number of primaries "<<nprim<<endl;
656   int co=0;
657   for(Int_t i=0; i<nparticles; i++)
658     {
659       TParticle *part = gAlice->Particle(i);
660       if(!part->TestBit(kPrimaryCharged)) continue;
661       if(fabs(part->Eta())>0.9) continue;
662       cout<<part->GetName()<<" pt "<<part->Pt()<<" eta "<<part->Eta()<<" xvert "<<part->Vx()<<" yvert "<<part->Vy()<<" zvert "<<part->Vz()<<endl;
663       co++;
664     }
665   cout<<endl<<"Number of primary tracks in the detector: "<<co<<endl;
666   gAlice=0;
667   rootfile->Close();
668   
669 }
670
671 Int_t FindPrimaries(Int_t nparticles, Double_t xori, Double_t yori, Double_t zori)
672 {
673   //Define primary particles in a pp-event. Code taken from offline.
674
675
676   // cuts:
677   //Double_t vertcut = 0.001;  // cut on the vertex position
678   Double_t vertcut = 2.0;
679   Double_t decacut = 3.;     // cut if the part. decays close to the vert.
680   Double_t timecut = 0.;
681   
682   TList *listprim = new TList();
683   listprim->SetOwner(kFALSE);
684   
685   Int_t nprch1=0;
686   Int_t nprch2=0;
687   for(Int_t iprim = 0; iprim<nparticles; iprim++){   //loop on  tracks
688     
689     TParticle * part = gAlice->Particle(iprim);
690     char *xxx=strstr(part->GetName(),"XXX");
691     if(xxx)continue;
692     
693     TParticlePDG *ppdg = part->GetPDG();
694     //if(TMath::Abs(ppdg->Charge())!=3)continue;  // only charged (no quarks)
695     if(TMath::Abs(ppdg->Charge())<1)continue;  // only charged (no quarks)
696     
697     Double_t dist=TMath::Sqrt((part->Vx()-xori)*(part->Vx()-xori)+(part->Vy()-yori)*(part->Vy()-yori)+(part->Vz()-zori)*(part->Vz()-zori));
698     if(dist>vertcut)continue;  // cut on the vertex
699
700     if(part->T()>timecut)continue;
701
702     Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
703     if(ptot==(TMath::Abs(part->Pz())))continue; // no beam particles
704
705     Bool_t prmch = kTRUE;   // candidate primary track
706     Int_t fidau=part->GetFirstDaughter();  // cut on daughters
707     Int_t lasdau=0;
708     Int_t ndau=0;
709     if(fidau>=0){
710       lasdau=part->GetLastDaughter();
711       ndau=lasdau-fidau+1;
712     }
713     if(ndau>0){
714       for(Int_t j=fidau;j<=lasdau;j++){
715         TParticle *dau=gAlice->Particle(j);
716         Double_t distd=TMath::Sqrt((dau->Vx()-xori)*(dau->Vx()-xori)+(dau->Vy()-yori)*(dau->Vy()-yori)+(dau->Vz()-zori)*(dau->Vz()-zori));
717         Double_t rad = TMath::Sqrt((dau->Vx()-xori)*(dau->Vx()-xori)+(dau->Vy()-yori)*(dau->Vy()-yori));
718         if(distd<decacut)prmch=kFALSE;  // eliminate if the decay is near the vertex
719         //if(rad < 20)prmch=kFALSE;
720       }
721     }
722
723     if(prmch){
724       nprch1++;
725       part->SetBit(kPrimaryCharged);
726       listprim->Add(part);    // list of primary particles (before cleanup)
727     }
728   }
729
730
731   nprch2=0;
732   for(Int_t iprim = 0; iprim<nparticles; iprim++){ // cleanup loop
733     TParticle * part = gAlice->Particle(iprim);
734     if(part->TestBit(kPrimaryCharged)){
735       Int_t mothind=part->GetFirstMother();
736       if(mothind<0)continue;
737       TParticle *moth=gAlice->Particle(mothind);
738       TParticle *mothb=(TParticle*)listprim->FindObject(moth);
739       if(mothb){
740         listprim->Remove(moth);
741         moth->ResetBit(kPrimaryCharged);
742         nprch2++;
743       }
744     }
745   }
746
747   listprim->Clear("nodelete");
748   delete listprim;
749   return nprch1-nprch2;
750   
751 }
752 #endif