Using TMath::Abs instead of fabs
[u/mrichter/AliRoot.git] / HLT / exa / trigger_pp.C
1 // $Id$
2
3 void trigger_pp(char *outfile="results.root")
4 {
5   
6   TNtuple *ntuppel = new TNtuple("ntuppel","","pt:eta:xvert:yvert:zvert:nhits:px:py:pz:event");
7   Float_t meas[10];  
8   
9   for(int event=0; event<1; event++)
10     {
11       char fname[256];
12       sprintf(fname,"/prog/alice/data/Rawdata/1_patch/pp/pileups/recon_%d/tracks.raw",event);
13       //sprintf(fname,"/prog/alice/data/Rawdata/1_patch/pp/recon_%d/tracks.raw",event);
14
15       //Get the tracks:
16       AliL3TrackArray *tracks = new AliL3TrackArray();
17       AliL3FileHandler *file = new AliL3FileHandler();
18       file->SetBinaryInput(fname);
19       file->Binary2TrackArray(tracks);
20       file->CloseBinaryInput();
21       delete file;
22       
23       sprintf(fname,"/prog/alice/data/Rawdata/1_patch/pp/pileups/recon_%d/",event);
24       //sprintf(fname,"/prog/alice/data/Rawdata/1_patch/pp/recon_%d/",event);
25
26       Int_t ntracks=0;
27       Double_t xc,yc,zc;
28       Double_t impact;
29       AliL3Vertex vertex;
30       AliL3Fitter *fitter = new AliL3Fitter(&vertex);
31       fitter->LoadClusters(fname);
32       //fitter->NoVertex();
33       
34       AliL3TrackArray *ftracks = new AliL3TrackArray();
35       
36       for(int i=0; i<tracks->GetNTracks(); i++)
37         {
38           track = (AliL3Track*)tracks->GetCheckedTrack(i);
39           if(!track) continue;
40           track->CalculateHelix();
41           //cout<<"Pt "<<track->GetPt()<<" eta "<<track->GetPseudoRapidity()<<" Nhits "<<track->GetNHits()<<endl;
42           //cout<<"Before second fit; pt "<<track->GetPt()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
43           fitter->FitHelix(track);//refit the tracks
44           track->CalculateHelix();
45           track->GetClosestPoint(&vertex,xc,yc,zc);
46           meas[0]=track->GetPt();
47           meas[1]=track->GetPseudoRapidity();
48           meas[2]=xc;
49           meas[3]=yc;
50           meas[4]=zc;
51           meas[5]=track->GetNHits();
52           meas[6]=track->GetPx();
53           meas[7]=track->GetPy();
54           meas[8]=track->GetPz();
55           meas[9]=event;
56           ntuppel->Fill(meas);
57           //continue;
58           if(fabs(track->GetPseudoRapidity())>0.9) continue;
59           if(track->GetNHits() < 100) continue;
60           if(track->GetPt()<0.2) continue;
61           impact = sqrt(xc*xc+yc*yc+zc*zc);
62           if(fabs(zc)>3) continue;
63           ftracks->AddLast(track);
64           cout<<"-------------------------------------"<<endl;
65           cout<<"Number of hits "<<track->GetNHits()<<endl;
66           cout<<"Transversal impact "<<sqrt(xc*xc+yc*yc)<<endl;
67           cout<<"Longitudinal impact "<<zc<<endl;
68           cout<<"Total impact "<<sqrt(xc*xc+yc*yc+zc*zc)<<endl;
69           cout<<"xc "<<xc<<" yc "<<yc<<" zc "<<zc<<endl;
70           
71           //cout<<"After second fit; pt "<<track->GetPt()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
72           cout<<"pt "<<track->GetPt()<<" eta "<<track->GetPseudoRapidity()<<endl;
73           
74           ntracks++;
75         }
76       cout<<endl<<"There was "<<ntracks<<" accepted tracks, out of total "<<tracks->GetNTracks()<<" found tracks"<<endl;
77       
78       display(ftracks,fname);
79       delete tracks;
80       delete fitter;
81       delete ftracks;
82     }
83   
84   TFile *rfile = TFile::Open(outfile,"RECREATE");
85   ntuppel->Write();
86   rfile->Close();
87   delete ntuppel;
88
89 }
90
91 void display(AliL3TrackArray *tracks,char *path)
92 {
93   int slice[2]={0,35};
94   d = new AliL3Display(slice);
95   d->Setup("tracks.raw",path);
96   d->SetTracks(tracks);
97   //d->DisplayClusters();
98     d->DisplayAll();
99     //d->DisplayTracks();
100   
101 }
102
103 void ploteff()
104 {
105   gROOT->LoadMacro("XFunct.C");
106   //double x[6]={1,2,4,6,8,10};
107   //double y[6]={0.873684,1.01379,1.17751,1.28614,1.31638,1.32022};
108   
109   double x[6]={0.8,1.6,3.2,4.8,6.4,8};
110   double y[6]={0.497006,0.88024,1.19162,1.23952,1.39521,1.40719};//->sigmas
111
112   
113   TGraph *gr = new TGraph(6,x,y);
114   c1 = new TCanvas("c1","",2);
115   SetCanvasOptions(c1);
116   gr->SetMarkerStyle(20);
117   gr->SetTitle("");
118   gr->Draw("APL");
119   gr->GetHistogram()->SetXTitle("Z_{CUT} [x #sigma_{Z}]");
120   gr->GetHistogram()->SetYTitle("Efficiency");
121   SetTGraphOptions(gr);
122 }
123
124 double geteff(char *infile,char *singlefile,double cut)
125 {
126   gROOT->LoadMacro("XFunct.C");
127   gStyle->SetStatColor(10);
128   gStyle->SetOptFit(1);
129   Int_t pileups[25],single[25];
130   file = TFile::Open(infile,"READ");
131   
132   char name[256];
133   for(int i=0; i<25; i++)
134     {
135       sprintf(name,"zvert < %f && zvert > %f && pt>0.2 && nhits>100 && eta>-0.9 && eta<0.9 && event==%d",cut,-1.*cut,i);
136       ntuppel->Draw(">>eventlist",name);
137       int entries = eventlist->GetN();
138       pileups[i]=entries;
139     }
140   //eventlist->Print("all");
141   file->Close();
142   
143   file = TFile::Open(singlefile,"read");
144   for(int i=0; i<25; i++)
145     {
146       sprintf(name,"zvert < %f && zvert > %f && pt>0.2 && nhits>100 && eta>-0.9 && eta<0.9 && event==%d",3,-3,i);
147       ntuppel->Draw(">>eventlist",name);
148       int entries = eventlist->GetN();
149       single[i]=entries;
150     }
151   file->Close();
152   double totsingle=0,totpileup=0;
153   for(int i=0; i<25; i++)
154     {
155       totsingle += single[i];
156       totpileup += pileups[i];
157     }
158   double toteff = totpileup/totsingle;
159   cout<<"Total eff "<<toteff<<endl;
160   
161   return toteff;
162 }
163
164
165 void plot(char *infile)
166 {
167   gROOT->LoadMacro("XFunct.C");
168   gStyle->SetStatColor(10);
169   gStyle->SetOptFit(1);
170   file = TFile::Open(infile,"READ");
171   
172   histo = new TH1F("histo","histo",20,-3,3);
173   
174   vhist = new TH1F("vhist","",20,-3,3);
175   SetTH1Options(vhist);
176   
177   char fname[256];
178   char fname2[256];
179   for(int ev=0; ev<25; ev++)
180     {
181       sprintf(fname,"pt>0.2 && nhits>100 && eta>-0.9 && eta<0.9 && event==%d",ev);
182       ntuppel->Draw("zvert>>histo",fname,"goff");
183       float mean = histo->GetMean();
184       vhist->Fill(mean);
185       continue;
186       sprintf(fname2,"(zvert-(%f))>>+vhist",mean);
187       cout<<fname2<<endl;
188       ntuppel->Draw(fname2,fname,"goff");
189     }
190   
191   c1 = new TCanvas("c1","",2);
192   SetCanvasOptions(c1);
193   vhist->SetXTitle("Z* [cm]");
194   vhist->Draw();
195   return;
196   //ntuppel->Draw("zvert>>histo","pt>0.2");
197   TF1 *f1 = new TF1("f1","gaus",-3,3);
198   histo->Fit("f1","R");
199
200   //histo->Draw();
201   //file->Close();
202 }
203
204 enum tagprimary {kPrimaryCharged=0x4000};
205 void LoadEvent(Int_t event=0)
206 {
207   //Load the generated particles
208   
209   gROOT->LoadMacro("$(ALICE_ROOT)/macros/loadlibs.C");
210   loadlibs();
211   TFile *rootfile = TFile::Open("/prog/alice/data/pro/25event_pp.root","READ");
212   gAlice = (AliRun*)rootfile->Get("gAlice");
213   
214   //  TNtuple *ntup = new TNtuple(
215   
216   Int_t nparticles = gAlice->GetEvent(event);
217   Int_t nprim = FindPrimaries(nparticles,0.,0.,0.);
218   cout<<"Number of primaries "<<nprim<<endl;
219   int co=0;
220   for(Int_t i=0; i<nparticles; i++)
221     {
222       TParticle *part = gAlice->Particle(i);
223       if(!part->TestBit(kPrimaryCharged)) continue;
224       if(fabs(part->Eta())>0.9) continue;
225       cout<<part->GetName()<<" pt "<<part->Pt()<<" eta "<<part->Eta()<<" xvert "<<part->Vx()<<" yvert "<<part->Vy()<<" zvert "<<part->Vz()<<endl;
226       co++;
227     }
228   cout<<endl<<"Number of primary tracks in the detector: "<<co<<endl;
229   gAlice=0;
230   rootfile->Close();
231   
232 }
233
234 Int_t FindPrimaries(Int_t nparticles, Double_t xori, Double_t yori, Double_t zori)
235 {
236   //Define primary particles in a pp-event. Code taken from offline.
237
238
239   // cuts:
240   Double_t vertcut = 0.001;  // cut on the vertex position
241   Double_t decacut = 3.;     // cut if the part. decays close to the vert.
242   Double_t timecut = 0.;
243   
244   TList *listprim = new TList();
245   listprim->SetOwner(kFALSE);
246   
247   Int_t nprch1=0;
248   Int_t nprch2=0;
249   for(Int_t iprim = 0; iprim<nparticles; iprim++){   //loop on  tracks
250     
251     TParticle * part = gAlice->Particle(iprim);
252     char *xxx=strstr(part->GetName(),"XXX");
253     if(xxx)continue;
254     
255     TParticlePDG *ppdg = part->GetPDG();
256     if(TMath::Abs(ppdg->Charge())!=3)continue;  // only charged (no quarks)
257     
258     Double_t dist=TMath::Sqrt((part->Vx()-xori)*(part->Vx()-xori)+(part->Vy()-yori)*(part->Vy()-yori)+(part->Vz()-zori)*(part->Vz()-zori));
259     if(dist>vertcut)continue;  // cut on the vertex
260
261     if(part->T()>timecut)continue;
262
263     Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
264     if(ptot==(TMath::Abs(part->Pz())))continue; // no beam particles
265
266     Bool_t prmch = kTRUE;   // candidate primary track
267     Int_t fidau=part->GetFirstDaughter();  // cut on daughters
268     Int_t lasdau=0;
269     Int_t ndau=0;
270     if(fidau>=0){
271       lasdau=part->GetLastDaughter();
272       ndau=lasdau-fidau+1;
273     }
274     if(ndau>0){
275       for(Int_t j=fidau;j<=lasdau;j++){
276         TParticle *dau=gAlice->Particle(j);
277         Double_t distd=TMath::Sqrt((dau->Vx()-xori)*(dau->Vx()-xori)+(dau->Vy()-yori)*(dau->Vy()-yori)+(dau->Vz()-zori)*(dau->Vz()-zori));
278         Double_t rad = TMath::Sqrt((dau->Vx()-xori)*(dau->Vx()-xori)+(dau->Vy()-yori)*(dau->Vy()-yori));
279         if(distd<decacut)prmch=kFALSE;  // eliminate if the decay is near the vertex
280         //if(rad < 20)prmch=kFALSE;
281       }
282     }
283
284     if(prmch){
285       nprch1++;
286       part->SetBit(kPrimaryCharged);
287       listprim->Add(part);    // list of primary particles (before cleanup)
288     }
289   }
290
291
292   nprch2=0;
293   for(Int_t iprim = 0; iprim<nparticles; iprim++){ // cleanup loop
294     TParticle * part = gAlice->Particle(iprim);
295     if(part->TestBit(kPrimaryCharged)){
296       Int_t mothind=part->GetFirstMother();
297       if(mothind<0)continue;
298       TParticle *moth=gAlice->Particle(mothind);
299       TParticle *mothb=(TParticle*)listprim->FindObject(moth);
300       if(mothb){
301         listprim->Remove(moth);
302         moth->ResetBit(kPrimaryCharged);
303         nprch2++;
304       }
305     }
306   }
307
308   listprim->Clear("nodelete");
309   delete listprim;
310   return nprch1-nprch2;
311   
312 }