Using TMath::Abs instead of fabs
[u/mrichter/AliRoot.git] / HLT / exa / trigger_pp.C
CommitLineData
086f41d8 1// $Id$
2
edd80d97 3void trigger_pp(char *outfile="results.root")
b7556f61 4{
5
edd80d97 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++)
b7556f61 10 {
edd80d97 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();
b7556f61 33
edd80d97 34 AliL3TrackArray *ftracks = new AliL3TrackArray();
b7556f61 35
edd80d97 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;
b7556f61 77
edd80d97 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
91void 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
103void 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
124double 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;
b7556f61 150 }
edd80d97 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
165void 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();
b7556f61 202}
203
204enum tagprimary {kPrimaryCharged=0x4000};
205void LoadEvent(Int_t event=0)
206{
207 //Load the generated particles
edd80d97 208
b7556f61 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");
edd80d97 213
214 // TNtuple *ntup = new TNtuple(
215
b7556f61 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
234Int_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.;
edd80d97 243
b7556f61 244 TList *listprim = new TList();
245 listprim->SetOwner(kFALSE);
edd80d97 246
b7556f61 247 Int_t nprch1=0;
248 Int_t nprch2=0;
249 for(Int_t iprim = 0; iprim<nparticles; iprim++){ //loop on tracks
edd80d97 250
b7556f61 251 TParticle * part = gAlice->Particle(iprim);
252 char *xxx=strstr(part->GetName(),"XXX");
253 if(xxx)continue;
edd80d97 254
b7556f61 255 TParticlePDG *ppdg = part->GetPDG();
256 if(TMath::Abs(ppdg->Charge())!=3)continue; // only charged (no quarks)
edd80d97 257
b7556f61 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));
edd80d97 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;
b7556f61 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}