]>
Commit | Line | Data |
---|---|---|
3e87ef69 | 1 | // @(#) $Id$ |
e8541476 | 2 | |
b661165c | 3 | // Author: Anders Vestbo <mailto:vestbo@fi.uib.no> |
3e87ef69 | 4 | //*-- Copyright © ALICE HLT Group |
108615fc | 5 | |
118c26c3 | 6 | #include "AliL3StandardIncludes.h" |
108615fc | 7 | #include <TFile.h> |
8 | #include <TH1.h> | |
9 | #include <TParticle.h> | |
10 | #include <TTree.h> | |
11 | #include <TClonesArray.h> | |
12 | ||
4499ed26 | 13 | #include <AliRun.h> |
14 | #include <AliSimDigits.h> | |
15 | #include <AliTPC.h> | |
16 | #include <AliTPCcluster.h> | |
17 | #include <AliTPCClustersArray.h> | |
18 | #include <AliTPCClustersRow.h> | |
19 | #include <AliTPCParam.h> | |
20 | #include <AliComplexCluster.h> | |
b2a02bce | 21 | #include <AliStack.h> |
108615fc | 22 | |
0bd0c1ef | 23 | #if __GNUC__ == 3 |
118c26c3 | 24 | #include <fstream> |
25 | #include <iosfwd> | |
26 | #else | |
27 | #include <fstream.h> | |
28 | #endif | |
29 | ||
30 | #include "AliL3Logging.h" | |
108615fc | 31 | #include "AliL3Transform.h" |
32 | #include "AliL3SpacePointData.h" | |
33 | #include "AliL3Track.h" | |
34 | #include "AliL3FileHandler.h" | |
35 | #include "AliL3TrackArray.h" | |
36 | #include "AliL3Evaluate.h" | |
118c26c3 | 37 | |
0bd0c1ef | 38 | #if __GNUC__ == 3 |
118c26c3 | 39 | using namespace std; |
40 | #endif | |
108615fc | 41 | |
3e87ef69 | 42 | /** \class AliL3Evaluate |
43 | <pre> | |
b661165c | 44 | //_____________________________________________________________ |
45 | // AliL3Evaluate | |
46 | // | |
47 | // Evaluation class for tracking. Plots efficiencies etc.. | |
84122e15 | 48 | // |
3e87ef69 | 49 | </pre> |
50 | */ | |
108615fc | 51 | |
52 | ClassImp(AliL3Evaluate) | |
53 | ||
54 | AliL3Evaluate::AliL3Evaluate() | |
55 | { | |
e0f4d6b2 | 56 | fTracks = 0; |
c0217812 | 57 | fNFastPoints = 0; |
58 | fMcIndex = 0; | |
59 | fMcId = 0; | |
1f79afc0 | 60 | fMinSlice=0; |
61 | fMaxSlice=0; | |
e0f4d6b2 | 62 | fGoodTracks=0; |
b2a02bce | 63 | fNtupleRes=0; |
64 | fDigitsTree=0; | |
65 | fPtRes=0; | |
66 | fNGoodTracksPt=0; | |
67 | fNFoundTracksPt=0; | |
68 | fNFakeTracksPt=0; | |
69 | fTrackEffPt=0; | |
70 | fFakeTrackEffPt=0; | |
71 | fNGoodTracksEta=0; | |
72 | fNFoundTracksEta=0; | |
73 | fNFakeTracksEta=0; | |
74 | fTrackEffEta=0; | |
75 | fFakeTrackEffEta=0; | |
76 | fMcIndex=0; | |
77 | fMcId=0; | |
78 | fNtuppel=0; | |
79 | fStandardComparison=kTRUE; | |
108615fc | 80 | } |
81 | ||
b2a02bce | 82 | AliL3Evaluate::AliL3Evaluate(Char_t *datapath,Int_t min_clusters,Int_t minhits,Double_t minpt,Double_t maxpt,Int_t *slice) |
108615fc | 83 | { |
3e87ef69 | 84 | |
e0f4d6b2 | 85 | if(slice) |
86 | { | |
87 | fMinSlice=slice[0]; | |
88 | fMaxSlice=slice[1]; | |
89 | } | |
90 | else | |
91 | { | |
92 | fMinSlice=0; | |
93 | fMaxSlice=35; | |
94 | } | |
3e87ef69 | 95 | sprintf(fPath,"%s",datapath); |
e0f4d6b2 | 96 | fMaxFalseClusters = 0.1; |
23908b9b | 97 | fGoodFound = 0; |
98 | fGoodGen = 0; | |
e0f4d6b2 | 99 | fMinPointsOnTrack = min_clusters; |
100 | fMinHitsFromParticle = minhits; | |
101 | fMinGoodPt = minpt; | |
b2a02bce | 102 | fMaxGoodPt = maxpt; |
3e87ef69 | 103 | memset(fClusters,0,36*6*sizeof(AliL3SpacePointData*)); |
104 | fTracks=0; | |
105 | fGoodTracks=0; | |
b2a02bce | 106 | fNtupleRes=0; |
107 | fDigitsTree=0; | |
108 | fPtRes=0; | |
109 | fNGoodTracksPt=0; | |
110 | fNFoundTracksPt=0; | |
111 | fNFakeTracksPt=0; | |
112 | fTrackEffPt=0; | |
113 | fFakeTrackEffPt=0; | |
114 | fNGoodTracksEta=0; | |
115 | fNFoundTracksEta=0; | |
116 | fNFakeTracksEta=0; | |
117 | fTrackEffEta=0; | |
118 | fFakeTrackEffEta=0; | |
119 | fMcIndex=0; | |
120 | fMcId=0; | |
121 | fNtuppel=0; | |
122 | fStandardComparison=kTRUE; | |
3e87ef69 | 123 | } |
124 | ||
125 | void AliL3Evaluate::LoadData(Int_t event,Bool_t sp) | |
126 | { | |
e0f4d6b2 | 127 | Char_t fname[1024]; |
40896357 | 128 | AliL3FileHandler *clusterfile[36][6]; |
108615fc | 129 | for(Int_t s=fMinSlice; s<=fMaxSlice; s++) |
130 | { | |
e0f4d6b2 | 131 | for(Int_t p=0; p<AliL3Transform::GetNPatches(); p++) |
108615fc | 132 | { |
3e87ef69 | 133 | Int_t patch; |
134 | if(sp==kTRUE) | |
135 | patch=-1; | |
136 | else | |
137 | patch=p; | |
138 | ||
139 | if(fClusters[s][p]) | |
140 | delete fClusters[s][p]; | |
1f79afc0 | 141 | fClusters[s][p] = 0; |
108615fc | 142 | clusterfile[s][p] = new AliL3FileHandler(); |
3e87ef69 | 143 | if(event<0) |
144 | sprintf(fname,"%s/points_%d_%d.raw",fPath,s,patch); | |
145 | else | |
146 | sprintf(fname,"%s/points_%d_%d_%d.raw",fPath,event,s,patch); | |
108615fc | 147 | if(!clusterfile[s][p]->SetBinaryInput(fname)) |
148 | { | |
149 | LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open") | |
150 | <<"Inputfile "<<fname<<" does not exist"<<ENDLOG; | |
bc2f4f0e | 151 | delete clusterfile[s][p]; |
152 | clusterfile[s][p] = 0; | |
153 | continue; | |
108615fc | 154 | } |
155 | fClusters[s][p] = (AliL3SpacePointData*)clusterfile[s][p]->Allocate(); | |
156 | clusterfile[s][p]->Binary2Memory(fNcl[s][p],fClusters[s][p]); | |
157 | clusterfile[s][p]->CloseBinaryInput(); | |
3e87ef69 | 158 | if(sp==kTRUE) |
159 | break; | |
108615fc | 160 | } |
161 | } | |
b2a02bce | 162 | |
3e87ef69 | 163 | sprintf(fname,"%s/tracks_%d.raw",fPath,event); |
108615fc | 164 | AliL3FileHandler *tfile = new AliL3FileHandler(); |
e0f4d6b2 | 165 | if(!tfile->SetBinaryInput(fname)){ |
bc2f4f0e | 166 | LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open") |
e0f4d6b2 | 167 | <<"Inputfile "<<fname<<" does not exist"<<ENDLOG; |
bc2f4f0e | 168 | return; |
169 | } | |
3e87ef69 | 170 | if(fTracks) |
171 | delete fTracks; | |
108615fc | 172 | fTracks = new AliL3TrackArray(); |
173 | tfile->Binary2TrackArray(fTracks); | |
174 | tfile->CloseBinaryInput(); | |
175 | delete tfile; | |
b2a02bce | 176 | fTracks->QSort(); |
bc2f4f0e | 177 | } |
178 | ||
108615fc | 179 | |
e0f4d6b2 | 180 | AliL3Evaluate::~AliL3Evaluate() |
108615fc | 181 | { |
e0f4d6b2 | 182 | if(fGoodTracks) delete fGoodTracks; |
183 | if(fDigitsTree) fDigitsTree->Delete(); | |
184 | if(fTracks) delete fTracks; | |
185 | if(fPtRes) delete fPtRes; | |
186 | if(fNGoodTracksPt) delete fNGoodTracksPt; | |
187 | if(fNFoundTracksPt) delete fNFoundTracksPt; | |
188 | if(fNFakeTracksPt) delete fNFakeTracksPt; | |
189 | if(fTrackEffPt) delete fTrackEffPt; | |
190 | if(fFakeTrackEffPt) delete fFakeTrackEffPt; | |
191 | if(fNGoodTracksEta) delete fNGoodTracksEta; | |
192 | if(fNFoundTracksEta) delete fNFoundTracksEta; | |
193 | if(fNFakeTracksEta) delete fNFakeTracksEta; | |
194 | if(fTrackEffEta) delete fTrackEffEta; | |
195 | if(fFakeTrackEffEta) delete fFakeTrackEffEta; | |
196 | if(fMcIndex) delete [] fMcIndex; | |
197 | if(fMcId) delete [] fMcId; | |
198 | if(fNtuppel) delete fNtuppel; | |
b2a02bce | 199 | if(fNtupleRes) delete fNtupleRes; |
108615fc | 200 | } |
201 | ||
8af8a456 | 202 | |
108615fc | 203 | void AliL3Evaluate::AssignIDs() |
204 | { | |
205 | //Assign MC id to the tracks. | |
e0f4d6b2 | 206 | #ifndef do_mc |
207 | cerr<<"AliL3Evaluate::AssignIDs() : You need to compile with the do_mc flag!"<<endl; | |
208 | return; | |
209 | #endif | |
3e87ef69 | 210 | fGoodFound=0; |
108615fc | 211 | fTracks->QSort(); |
212 | LOG(AliL3Log::kDebug,"AliL3Evaluate::AssignIDs","Track Loop") | |
213 | <<"Assigning MC id to the found tracks...."<<ENDLOG; | |
214 | for(Int_t i=0; i<fTracks->GetNTracks(); i++) | |
215 | { | |
216 | AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i); | |
c0217812 | 217 | if(!track) continue; |
108615fc | 218 | if(track->GetNumberOfPoints() < fMinPointsOnTrack) break; |
23908b9b | 219 | |
23908b9b | 220 | fGoodFound++; |
c0217812 | 221 | Int_t tID = GetMCTrackLabel(track); |
108615fc | 222 | track->SetMCid(tID); |
108615fc | 223 | } |
3e87ef69 | 224 | //cout<<"Found "<<fGoodFound<<" good tracks "<<endl; |
108615fc | 225 | } |
226 | ||
227 | ||
228 | struct S {Int_t lab; Int_t max;}; | |
c0217812 | 229 | Int_t AliL3Evaluate::GetMCTrackLabel(AliL3Track *track){ |
108615fc | 230 | //Returns the MCtrackID of the belonging clusters. |
231 | //If MCLabel < 0, means that track is fake. | |
3e87ef69 | 232 | //Definitions are identical to offline. |
233 | //Fake track means: | |
234 | // - more than 10 percent of clusters are assigned incorrectly | |
235 | // - more than half of the innermost 10% of clusters were assigned incorrectly. | |
236 | ||
494fad94 | 237 | |
238 | #ifdef do_mc | |
108615fc | 239 | Int_t num_of_clusters = track->GetNumberOfPoints(); |
240 | S *s=new S[num_of_clusters]; | |
241 | Int_t i; | |
242 | for (i=0; i<num_of_clusters; i++) s[i].lab=s[i].max=0; | |
243 | UInt_t *hitnum = track->GetHitNumbers(); | |
244 | UInt_t id; | |
108615fc | 245 | |
108615fc | 246 | Int_t lab=123456789; |
247 | for (i=0; i<num_of_clusters; i++) | |
248 | { | |
249 | //Tricks to get the clusters belonging to this track: | |
250 | id = hitnum[i]; | |
251 | Int_t slice = (id>>25) & 0x7f; | |
252 | Int_t patch = (id>>22) & 0x7; | |
253 | UInt_t pos = id&0x3fffff; | |
254 | ||
255 | AliL3SpacePointData *points = fClusters[slice][patch]; | |
c0217812 | 256 | if(!points) continue; |
108615fc | 257 | if(pos>=fNcl[slice][patch]) |
258 | { | |
259 | LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray") | |
260 | <<AliL3Log::kDec<<"ERROR"<<ENDLOG; | |
261 | continue; | |
262 | } | |
263 | ||
264 | //Get the label of the cluster: | |
3e87ef69 | 265 | lab=points[pos].fTrackID[0]; |
266 | ||
108615fc | 267 | Int_t j; |
268 | for (j=0; j<num_of_clusters; j++) | |
269 | if (s[j].lab==lab || s[j].max==0) break; | |
270 | s[j].lab=lab; | |
271 | s[j].max++; | |
272 | } | |
273 | ||
274 | Int_t max=0; | |
275 | for (i=0; i<num_of_clusters; i++) | |
276 | if (s[i].max>max) {max=s[i].max; lab=s[i].lab;} | |
277 | ||
3e87ef69 | 278 | if(lab == -1) |
279 | return -1; //If most clusters is -1, this is a noise track. | |
280 | if(lab < 0) | |
281 | cerr<<"AliL3Evaluate::GetMCTrackLabel : Track label negative :"<<lab<<endl; | |
282 | ||
108615fc | 283 | delete[] s; |
284 | ||
285 | for (i=0; i<num_of_clusters; i++) | |
286 | { | |
287 | id = hitnum[i]; | |
288 | Int_t slice = (id>>25) & 0x7f; | |
289 | Int_t patch = (id>>22) & 0x7; | |
290 | UInt_t pos = id&0x3fffff; | |
291 | ||
292 | AliL3SpacePointData *points = fClusters[slice][patch]; | |
c0217812 | 293 | if(!points) continue; |
108615fc | 294 | if(pos>=fNcl[slice][patch]) |
295 | { | |
296 | LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray") | |
297 | <<AliL3Log::kDec<<"ERROR"<<ENDLOG; | |
298 | continue; | |
299 | } | |
300 | ||
355debe1 | 301 | if (abs(points[pos].fTrackID[1]) == lab || |
302 | abs(points[pos].fTrackID[2]) == lab ) max++; | |
108615fc | 303 | } |
304 | ||
3e87ef69 | 305 | |
306 | //Check if more than 10% of the clusters were assigned incorrectly: | |
e0f4d6b2 | 307 | if (1.-Float_t(max)/num_of_clusters > fMaxFalseClusters) |
108615fc | 308 | { |
309 | return -lab; | |
310 | } | |
3e87ef69 | 311 | else //Check if at least half of the 10% innermost clusters are assigned correctly. |
312 | { | |
313 | Int_t tail=Int_t(0.10*num_of_clusters); | |
314 | max=0; | |
315 | for (i=1; i<=tail; i++) | |
316 | { | |
317 | id = hitnum[num_of_clusters - i]; | |
318 | Int_t slice = (id>>25) & 0x7f; | |
319 | Int_t patch = (id>>22) & 0x7; | |
320 | UInt_t pos = id&0x3fffff; | |
321 | ||
322 | AliL3SpacePointData *points = fClusters[slice][patch]; | |
323 | if(lab == abs(points[pos].fTrackID[0]) || | |
324 | lab == abs(points[pos].fTrackID[1]) || | |
325 | lab == abs(points[pos].fTrackID[2])) max++; | |
326 | } | |
327 | if (max < Int_t(0.5*tail)) return -lab; | |
328 | } | |
329 | ||
108615fc | 330 | return lab; |
494fad94 | 331 | #else //If we are running with mc_ids or not |
332 | return 0; | |
333 | #endif | |
334 | ||
108615fc | 335 | } |
336 | ||
c0217812 | 337 | void AliL3Evaluate::GetFastClusterIDs(Char_t *path) |
108615fc | 338 | { |
339 | //Get the MC id of space points in case of using the fast simulator. | |
c0217812 | 340 | char fname[256]; |
341 | sprintf(fname,"%s/point_mc.dat",path); | |
342 | FILE *infile = fopen(fname,"r"); | |
343 | if(!infile) return; | |
108615fc | 344 | Int_t hitid,hitmc,i; |
345 | ||
346 | for(i=0; ; i++) | |
347 | if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break; | |
108615fc | 348 | rewind(infile); |
c0217812 | 349 | fNFastPoints = i; |
350 | fMcId = new Int_t[fNFastPoints]; | |
351 | fMcIndex = new UInt_t[fNFastPoints]; | |
108615fc | 352 | |
c0217812 | 353 | for(i=0; i<fNFastPoints; i++) |
108615fc | 354 | { |
355 | if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break; | |
c0217812 | 356 | fMcId[i] = hitmc; |
357 | fMcIndex[i] = hitid; | |
108615fc | 358 | } |
359 | fclose(infile); | |
c0217812 | 360 | } |
361 | ||
e0f4d6b2 | 362 | void AliL3Evaluate::CreateHistos(Int_t nbin,Float_t xlow,Float_t xup) |
c0217812 | 363 | { |
364 | //Create the histograms | |
365 | ||
40896357 | 366 | LOG(AliL3Log::kInformational,"AliL3Evaluate::CreateHistos","Allocating") |
367 | <<"Creating histograms..."<<ENDLOG; | |
368 | ||
1f79afc0 | 369 | fNtuppel = new TNtuple("fNtuppel","Pt resolution","pt_gen:pt_found:nHits"); |
a8a3f5c2 | 370 | fNtuppel->SetDirectory(0); |
c0217812 | 371 | fPtRes = new TH1F("fPtRes","Relative Pt resolution",30,-10.,10.); |
372 | fNGoodTracksPt = new TH1F("fNGoodTracksPt","Good tracks vs pt",nbin,xlow,xup); | |
373 | fNFoundTracksPt = new TH1F("fNFoundTracksPt","Found tracks vs pt",nbin,xlow,xup); | |
374 | fNFakeTracksPt = new TH1F("fNFakeTracksPt","Fake tracks vs pt",nbin,xlow,xup); | |
375 | fTrackEffPt = new TH1F("fTrackEffPt","Tracking efficiency vs pt",nbin,xlow,xup); | |
376 | fFakeTrackEffPt = new TH1F("fFakeTrackEffPt","Efficiency for fake tracks vs pt",nbin,xlow,xup); | |
377 | ||
378 | fNGoodTracksEta = new TH1F("fNGoodTracksEta","Good tracks vs eta",20,-50,50); | |
379 | fNFoundTracksEta = new TH1F("fNFoundTracksEta","Found tracks vs eta",20,-50,50); | |
380 | fNFakeTracksEta = new TH1F("fNFakeTracksEta","Fake tracks vs eta",20,-50,50); | |
381 | fTrackEffEta = new TH1F("fTrackEffEta","Tracking efficienct vs eta",20,-50,50); | |
382 | fFakeTrackEffEta = new TH1F("fFakeTrackEffEta","Efficiency for fake tracks vs eta",20,-50,50); | |
e0f4d6b2 | 383 | |
384 | } | |
385 | ||
3e87ef69 | 386 | void AliL3Evaluate::GetGoodParticles(Char_t *path,Int_t event,Int_t *padrowrange) |
e0f4d6b2 | 387 | { |
388 | //Read the good particles from file. This file should already have been | |
389 | //generated by macro AliTPCComparison.C. | |
390 | ||
391 | Char_t filename[1024]; | |
3e87ef69 | 392 | if(event<0 && !padrowrange) |
a8a3f5c2 | 393 | sprintf(filename,"%s/good_tracks_tpc",path); |
3e87ef69 | 394 | else if(event>=0 && !padrowrange) |
395 | sprintf(filename,"%s/good_tracks_tpc_%d",path,event); | |
a8a3f5c2 | 396 | else |
3e87ef69 | 397 | sprintf(filename,"%s/good_tracks_tpc_%d_%d_%d",path,event,padrowrange[0],padrowrange[1]); |
e0f4d6b2 | 398 | ifstream in(filename); |
399 | if(!in) | |
400 | { | |
401 | cerr<<"AliL3Evaluate::GetGoodParticles : Problems opening file :"<<filename<<endl; | |
402 | return; | |
403 | } | |
404 | Int_t MaxTracks=20000; | |
3e87ef69 | 405 | if(fGoodTracks) |
406 | delete [] fGoodTracks; | |
407 | fGoodGen=0; | |
e0f4d6b2 | 408 | fGoodTracks = new GoodTrack[MaxTracks]; |
409 | ||
b2a02bce | 410 | if(fStandardComparison){ |
411 | while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>> | |
412 | fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>> | |
413 | fGoodTracks[fGoodGen].x>>fGoodTracks[fGoodGen].y>>fGoodTracks[fGoodGen].z) | |
414 | { | |
415 | fGoodTracks[fGoodGen].nhits=-1; | |
416 | fGoodTracks[fGoodGen].sector=-1; | |
417 | fGoodGen++; | |
418 | if (fGoodGen==MaxTracks) | |
419 | { | |
420 | cerr<<"AliL3Evaluate::GetGoodParticles : Too many good tracks !\n"; | |
421 | break; | |
422 | } | |
423 | } | |
424 | } else { | |
425 | while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>> | |
426 | fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>> | |
427 | fGoodTracks[fGoodGen].x>>fGoodTracks[fGoodGen].y >>fGoodTracks[fGoodGen].z>>fGoodTracks[fGoodGen].nhits>>fGoodTracks[fGoodGen].sector) | |
428 | { | |
429 | fGoodGen++; | |
430 | if (fGoodGen==MaxTracks) | |
431 | { | |
432 | cerr<<"AliL3Evaluate::GetGoodParticles : Too many good tracks !\n"; | |
433 | break; | |
434 | } | |
435 | } | |
436 | } | |
c0217812 | 437 | } |
eeddc64d | 438 | |
b2a02bce | 439 | //has to be modified for fakes. |
a8a3f5c2 | 440 | |
e1a32fa5 | 441 | void AliL3Evaluate::FillEffHistos() |
c0217812 | 442 | { |
e0f4d6b2 | 443 | if(!fGoodTracks) |
444 | { | |
445 | cerr<<"AliL3Evaluate::FillEffHistos : No good tracks"<<endl; | |
446 | return; | |
447 | } | |
3e87ef69 | 448 | //cout<<"Comparing "<<fGoodGen<<" good tracks ..."<<endl; |
e1a32fa5 | 449 | for(Int_t i=0; i<fGoodGen; i++) |
c0217812 | 450 | { |
e0f4d6b2 | 451 | //cout<<"Checking particle "<<i<<endl; |
b2a02bce | 452 | if(!fStandardComparison) |
453 | if(fGoodTracks[i].nhits < fMinHitsFromParticle) continue; | |
e0f4d6b2 | 454 | Float_t ptg = TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py); |
b2a02bce | 455 | if(ptg < fMinGoodPt || ptg > fMaxGoodPt) continue; |
e0f4d6b2 | 456 | Float_t pzg=fGoodTracks[i].pz; |
c0217812 | 457 | Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi(); |
b2a02bce | 458 | |
459 | //If we are only considering tracks on one side of the TPC: | |
460 | if(fMaxSlice <= 17) | |
461 | if(dipangle < 0) | |
462 | continue; | |
e0f4d6b2 | 463 | |
c0217812 | 464 | fNGoodTracksPt->Fill(ptg); |
465 | fNGoodTracksEta->Fill(dipangle); | |
466 | Int_t found = 0; | |
e0f4d6b2 | 467 | |
c0217812 | 468 | for(Int_t k=0; k<fTracks->GetNTracks(); k++) |
469 | { | |
470 | AliL3Track *track = fTracks->GetCheckedTrack(k); | |
471 | if(!track) continue; | |
472 | Int_t nHits = track->GetNumberOfPoints(); | |
473 | if(nHits < fMinPointsOnTrack) break; | |
c0217812 | 474 | Int_t tracklabel; |
475 | tracklabel = track->GetMCid(); | |
e0f4d6b2 | 476 | |
e1a32fa5 | 477 | if(TMath::Abs(tracklabel) != fGoodTracks[i].label) continue; |
c0217812 | 478 | found=1; |
c0217812 | 479 | Float_t pt=track->GetPt(); |
3e87ef69 | 480 | if(tracklabel == fGoodTracks[i].label) |
481 | { | |
482 | fNFoundTracksPt->Fill(ptg); | |
483 | fNFoundTracksEta->Fill(dipangle); | |
484 | fNtuppel->Fill(ptg,pt,nHits); | |
485 | fPtRes->Fill((pt-ptg)/ptg*100.); | |
486 | } | |
487 | else | |
488 | { | |
489 | fNFakeTracksPt->Fill(ptg); | |
490 | fNFakeTracksEta->Fill(dipangle); | |
491 | } | |
492 | //fPtRes->Fill((pt-ptg)/ptg*100.); | |
493 | //fNtuppel->Fill(ptg,pt,nHits); | |
c0217812 | 494 | break; |
495 | ||
496 | } | |
3e87ef69 | 497 | //if(!found) |
b2a02bce | 498 | //cout<<"Track "<<fGoodTracks[i].label<<" was not found"<<endl; |
c0217812 | 499 | } |
500 | } | |
eeddc64d | 501 | |
e0f4d6b2 | 502 | void AliL3Evaluate::FillEffHistosNAIVE() |
355debe1 | 503 | { |
504 | //Fill the efficiency histograms. | |
505 | ||
506 | cout<<endl<<"Note: Doing NAIVE evaluation "<<endl; | |
507 | for(Int_t i=0; i<fGoodGen; i++) | |
508 | { | |
b2a02bce | 509 | if(!fStandardComparison) |
510 | if(fGoodTracks[i].nhits < fMinHitsFromParticle) continue; | |
e0f4d6b2 | 511 | Double_t ptg=TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py); |
b2a02bce | 512 | if(ptg < fMinGoodPt || ptg > fMaxGoodPt) continue; |
e0f4d6b2 | 513 | Double_t pzg=fGoodTracks[i].pz; |
355debe1 | 514 | Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi(); |
515 | //printf("filling particle with pt %f and dipangle %f\n",ptg,dipangle); | |
516 | fNGoodTracksPt->Fill(ptg); | |
517 | fNGoodTracksEta->Fill(dipangle); | |
518 | ||
519 | } | |
520 | ||
521 | for(Int_t k=0; k<fTracks->GetNTracks(); k++) | |
522 | { | |
523 | AliL3Track *track = fTracks->GetCheckedTrack(k); | |
524 | if(!track) continue; | |
525 | Int_t nHits = track->GetNumberOfPoints(); | |
526 | if(nHits < fMinPointsOnTrack) break; | |
b2a02bce | 527 | if(track->GetPt()<fMinGoodPt || track->GetPt() > fMaxGoodPt) continue; |
355debe1 | 528 | if(fabs(track->GetPseudoRapidity())>0.9) continue; |
529 | ||
530 | fNFoundTracksPt->Fill(track->GetPt()); fNFoundTracksEta->Fill(track->GetPseudoRapidity()); | |
531 | //Float_t pt=track->GetPt(); | |
532 | //fPtRes->Fill((pt-ptg)/ptg*100.); | |
533 | //fNtuppel->Fill(ptg,pt,nHits); | |
534 | ||
535 | } | |
536 | } | |
c0217812 | 537 | |
3e87ef69 | 538 | void AliL3Evaluate::CalcEffHistos() |
539 | { | |
540 | ||
c0217812 | 541 | Stat_t ngood=fNGoodTracksPt->GetEntries(); |
542 | Stat_t nfound=fNFoundTracksPt->GetEntries(); | |
543 | Stat_t nfake=fNFakeTracksPt->GetEntries(); | |
3e87ef69 | 544 | |
c0217812 | 545 | LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency") |
546 | <<AliL3Log::kDec<<"There was "<<ngood<<" generated good tracks"<<ENDLOG; | |
547 | LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency") | |
548 | <<AliL3Log::kDec<<"Found "<<nfound<<" tracks"<<ENDLOG; | |
549 | LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency") | |
550 | <<AliL3Log::kDec<<"Integral efficiency is about "<<nfound/ngood*100<<ENDLOG; | |
551 | LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency") | |
552 | <<AliL3Log::kDec<<"Fake tracks relative is about "<<nfake/ngood*100<<ENDLOG; | |
3e87ef69 | 553 | //LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency") |
554 | //<<AliL3Log::kDec<<"Naive efficiency "<<(Double_t)fGoodFound/(Double_t)fGoodGen<<ENDLOG; | |
355debe1 | 555 | |
c0217812 | 556 | fNFoundTracksPt->Sumw2(); fNGoodTracksPt->Sumw2(); |
557 | fTrackEffPt->Divide(fNFoundTracksPt,fNGoodTracksPt,1,1.,"b"); | |
558 | fFakeTrackEffPt->Divide(fNFakeTracksPt,fNGoodTracksPt,1,1.,"b"); | |
559 | fTrackEffPt->SetMaximum(1.4); | |
560 | fTrackEffPt->SetXTitle("P_{T} [GeV]"); | |
561 | fTrackEffPt->SetLineWidth(2); | |
562 | fFakeTrackEffPt->SetFillStyle(3013); | |
563 | fTrackEffPt->SetLineColor(4); | |
564 | fFakeTrackEffPt->SetFillColor(2); | |
3e87ef69 | 565 | |
c0217812 | 566 | fNFoundTracksEta->Sumw2(); fNGoodTracksEta->Sumw2(); |
567 | fTrackEffEta->Divide(fNFoundTracksEta,fNGoodTracksEta,1,1.,"b"); | |
568 | fFakeTrackEffEta->Divide(fNFakeTracksEta,fNGoodTracksEta,1,1.,"b"); | |
569 | fTrackEffEta->SetMaximum(1.4); | |
570 | fTrackEffEta->SetXTitle("#lambda [degrees]"); | |
571 | fTrackEffEta->SetLineWidth(2); | |
572 | fFakeTrackEffEta->SetFillStyle(3013); | |
573 | fTrackEffEta->SetLineColor(4); | |
574 | fFakeTrackEffEta->SetFillColor(2); | |
575 | ||
108615fc | 576 | } |
577 | ||
578 | void AliL3Evaluate::Write2File(Char_t *outputfile) | |
579 | { | |
580 | //Write histograms to file: | |
581 | ||
e0f4d6b2 | 582 | TFile *of = TFile::Open(outputfile,"RECREATE"); |
108615fc | 583 | if(!of->IsOpen()) |
584 | { | |
585 | LOG(AliL3Log::kError,"AliL3Evaluate::Write2File","File Open") | |
586 | <<"Problems opening rootfile"<<ENDLOG; | |
587 | return; | |
588 | } | |
589 | ||
590 | of->cd(); | |
1f79afc0 | 591 | fNtuppel->Write(); |
108615fc | 592 | fPtRes->Write(); |
593 | fNGoodTracksPt->Write(); | |
594 | fNFoundTracksPt->Write(); | |
595 | fNFakeTracksPt->Write(); | |
596 | fTrackEffPt->Write(); | |
597 | fFakeTrackEffPt->Write(); | |
598 | fNGoodTracksEta->Write(); | |
599 | fNFoundTracksEta->Write(); | |
600 | fNFakeTracksEta->Write(); | |
601 | fTrackEffEta->Write(); | |
602 | fFakeTrackEffEta->Write(); | |
603 | ||
604 | of->Close(); | |
108615fc | 605 | |
606 | } | |
607 | ||
b2a02bce | 608 | TNtuple *AliL3Evaluate::GetNtuple() |
1f79afc0 | 609 | { |
b2a02bce | 610 | if(!fNtupleRes) |
611 | { | |
612 | fNtupleRes = new TNtuple("ntuppel","Residuals","residual_trans:residual_long:zHit:pt:dipangle:beta:padrow:nHits"); | |
613 | fNtupleRes->SetDirectory(0);//Bug in older version of root. | |
614 | } | |
615 | return fNtupleRes; | |
616 | } | |
1f79afc0 | 617 | |
b2a02bce | 618 | void AliL3Evaluate::CalculateResiduals() |
619 | { | |
620 | ||
621 | TNtuple *ntuppel = GetNtuple(); | |
3e87ef69 | 622 | |
623 | for(Int_t i=0; i<fTracks->GetNTracks(); i++) | |
1f79afc0 | 624 | { |
3e87ef69 | 625 | |
626 | AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i); | |
627 | if(!track) continue; | |
b2a02bce | 628 | if(track->GetNHits() < fMinPointsOnTrack) break; |
3e87ef69 | 629 | |
630 | track->CalculateHelix(); | |
631 | UInt_t *hitnum = track->GetHitNumbers(); | |
632 | UInt_t id; | |
633 | ||
634 | Float_t xyz[3]; | |
635 | Int_t padrow; | |
636 | for(Int_t j=0; j<track->GetNumberOfPoints()-1; j++) | |
1f79afc0 | 637 | { |
3e87ef69 | 638 | id = hitnum[j]; |
639 | Int_t slice = (id>>25) & 0x7f; | |
640 | Int_t patch = (id>>22) & 0x7; | |
641 | UInt_t pos = id&0x3fffff; | |
b2a02bce | 642 | |
643 | //if(slice<18) continue; | |
1f79afc0 | 644 | |
3e87ef69 | 645 | AliL3SpacePointData *points = fClusters[slice][patch]; |
1f79afc0 | 646 | |
3e87ef69 | 647 | if(!points) |
1f79afc0 | 648 | { |
3e87ef69 | 649 | LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray") |
650 | <<"No points at slice "<<slice<<" patch "<<patch<<" pos "<<pos<<ENDLOG; | |
651 | continue; | |
1f79afc0 | 652 | } |
3e87ef69 | 653 | if(pos>=fNcl[slice][patch]) |
654 | { | |
655 | LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray") | |
656 | <<AliL3Log::kDec<<"ERROR"<<ENDLOG; | |
657 | continue; | |
658 | } | |
659 | ||
660 | xyz[0] = points[pos].fX; | |
661 | xyz[1] = points[pos].fY; | |
662 | xyz[2] = points[pos].fZ; | |
663 | padrow = points[pos].fPadRow; | |
5a31e9df | 664 | //AliL3Transform::Global2Local(xyz,slice,kTRUE); |
665 | AliL3Transform::Global2LocHLT(xyz,slice); | |
3e87ef69 | 666 | |
667 | Float_t angle = 0; | |
668 | AliL3Transform::Local2GlobalAngle(&angle,slice); | |
b2a02bce | 669 | if(!track->CalculateReferencePoint(angle,AliL3Transform::Row2X(padrow))) |
670 | { | |
671 | LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Crossing point") | |
672 | <<"Track does not crossing padrow "<<padrow<<" in slice "<<slice<<ENDLOG; | |
673 | continue; | |
674 | } | |
675 | ||
3e87ef69 | 676 | Float_t xyz_cross[3] = {track->GetPointX(),track->GetPointY(),track->GetPointZ()}; |
5a31e9df | 677 | //AliL3Transform::Global2Local(xyz_cross,slice,kTRUE); |
678 | AliL3Transform::Global2LocHLT(xyz_cross,slice); | |
3e87ef69 | 679 | |
680 | Double_t beta = track->GetCrossingAngle(padrow,slice); | |
681 | ||
682 | Double_t yres = xyz_cross[1] - xyz[1]; | |
683 | Double_t zres = xyz_cross[2] - xyz[2]; | |
3e87ef69 | 684 | Double_t dipangle = atan(track->GetTgl()); |
685 | ntuppel->Fill(yres,zres,xyz_cross[2],track->GetPt(),dipangle,beta,padrow,track->GetNumberOfPoints()); | |
686 | ||
1f79afc0 | 687 | } |
1f79afc0 | 688 | } |
1f79afc0 | 689 | } |
690 | ||
3e87ef69 | 691 | enum tagprimary {kPrimaryCharged = 0x4000}; |
692 | void AliL3Evaluate::EvaluatePoints(Char_t *rootfile,Char_t *exactfile,Char_t *tofile,Int_t nevent,Bool_t offline,Bool_t sp) | |
78127e35 | 693 | { |
e1a32fa5 | 694 | //Compare points to the exact crossing points of track and padrows. |
1b76fd18 | 695 | //The input file to this function, contains the exact clusters calculated |
696 | //in AliTPC::Hits2ExactClusters. | |
8af8a456 | 697 | |
5bfdc6ad | 698 | #ifndef do_mc |
699 | cerr<<"AliL3Evaluate::EvaluatePoints : Compile with do_mc flag!"<<endl; | |
3e87ef69 | 700 | return; |
5bfdc6ad | 701 | #else |
355debe1 | 702 | cout<<"Evaluating points"<<endl; |
3e87ef69 | 703 | TNtuple *ntuppel = new TNtuple("ntuppel_res","Cluster properties", |
704 | "slice:padrow:charge:resy:resz:zHit:pt:beta:sigmaY2:sigmaZ2:psigmaY2:psigmaZ2"); | |
8af8a456 | 705 | ntuppel->SetDirectory(0); |
78127e35 | 706 | |
3e87ef69 | 707 | TNtuple *ntuppel2 = new TNtuple("ntuppel_eff","Efficiency","slice:padrow:nfound:ngen"); |
708 | ntuppel2->SetDirectory(0); | |
709 | ||
710 | TFile *exfile = TFile::Open(rootfile); | |
355debe1 | 711 | if(!exfile) |
712 | { | |
3e87ef69 | 713 | cerr<<"Error opening rootfile "<<rootfile<<endl; |
714 | return; | |
355debe1 | 715 | } |
8af8a456 | 716 | gAlice = (AliRun*)exfile->Get("gAlice"); |
717 | if (!gAlice) | |
718 | { | |
719 | LOG(AliL3Log::kError,"AliL3Evaluate::InitMC","gAlice") | |
720 | <<"AliRun object non existing on file"<<ENDLOG; | |
3e87ef69 | 721 | return; |
8af8a456 | 722 | } |
3e87ef69 | 723 | |
5e0f9911 | 724 | AliTPCParam *param = (AliTPCParam*)exfile->Get(AliL3Transform::GetParamName()); |
e1a32fa5 | 725 | |
3e87ef69 | 726 | TFile *exact = TFile::Open(exactfile); |
727 | if(!exact) | |
78127e35 | 728 | { |
3e87ef69 | 729 | cerr<<"AliL3Evaluate::EvaluatePoints : Problems opening file :"<<exactfile<<endl; |
730 | return; | |
731 | } | |
732 | ||
b2a02bce | 733 | AliStack *astack=gAlice->Stack(); |
734 | ||
3e87ef69 | 735 | AliTPCClustersArray *arr=0; |
736 | for(Int_t event=0; event<nevent; event++) | |
737 | { | |
738 | LoadData(event,sp); | |
739 | exfile->cd(); | |
740 | if(arr) | |
741 | delete arr; | |
742 | Int_t nparticles = gAlice->GetEvent(event); | |
743 | Int_t nprimaries = 0;//FindPrimaries(nparticles); | |
744 | cout<<"Event "<<event<<" had "<<nparticles<<" particles and "<<nprimaries<<" primaries"<<endl; | |
745 | exact->cd(); | |
78127e35 | 746 | |
3e87ef69 | 747 | //Get the exact clusters from file: |
748 | AliTPCClustersArray *arr = new AliTPCClustersArray; | |
749 | arr->Setup(param); | |
750 | arr->SetClusterType("AliComplexCluster"); | |
751 | char treeName[500]; | |
752 | sprintf(treeName,"TreeCExact_%s_%d",param->GetTitle(),event); | |
753 | Bool_t clusterok = arr->ConnectTree(treeName);//Segment Tree (for offline clusters) | |
754 | if(!clusterok) {printf("AliL3Evaluate::EvaluatePoints : Error in clusterloading\n"); return;} | |
78127e35 | 755 | |
3e87ef69 | 756 | //cout<<"Entering loop with "<<(Int_t)arr->GetTree()->GetEntries()<<endl; |
757 | for(Int_t i=0; i<arr->GetTree()->GetEntries(); i++) | |
8af8a456 | 758 | { |
3e87ef69 | 759 | //Get the exact clusters for this row: |
760 | Int_t cursec,currow; | |
761 | AliSegmentID *s = arr->LoadEntry(i); | |
762 | param->AdjustSectorRow(s->GetID(),cursec,currow); | |
763 | ||
764 | AliTPCClustersRow *ro = (AliTPCClustersRow *)arr->GetRow(cursec,currow); | |
765 | TClonesArray *clusters = ro->GetArray(); | |
766 | int num_of_offline=clusters->GetEntriesFast(); | |
767 | ||
768 | //Get the found clusters: | |
769 | Int_t slice,padrow; | |
770 | AliL3Transform::Sector2Slice(slice,padrow,cursec,currow); | |
771 | if(slice < fMinSlice) continue; | |
772 | if(slice > fMaxSlice) break; | |
773 | ||
774 | Int_t patch = AliL3Transform::GetPatch(padrow); | |
775 | if(sp) | |
776 | patch=0; | |
777 | AliL3SpacePointData *points = fClusters[slice][patch]; | |
778 | if(!points) | |
779 | continue; | |
780 | ||
781 | //cout<<"Slice "<<slice<<" padrow "<<padrow<<" has "<<num_of_offline<<" clusters "<<endl; | |
782 | Int_t clustercount=0; | |
783 | Int_t crosscount=0; | |
e1a32fa5 | 784 | for(Int_t m=0; m<num_of_offline; m++) |
78127e35 | 785 | { |
e1a32fa5 | 786 | AliComplexCluster *cluster = (AliComplexCluster *)clusters->UncheckedAt(m); |
b2a02bce | 787 | #ifdef use_newio |
788 | Int_t mcId = cluster->GetTrack(0); | |
789 | #else | |
e1a32fa5 | 790 | Int_t mcId = cluster->fTracks[0]; |
b2a02bce | 791 | #endif |
294ee877 | 792 | if(mcId <0) continue; |
b2a02bce | 793 | |
794 | #ifdef use_newio | |
795 | if(cluster->GetY() < 1 || cluster->GetY() > AliL3Transform::GetNPads(padrow) - 2 || | |
796 | cluster->GetX() < 1 || cluster->GetX() > AliL3Transform::GetNTimeBins() - 2) | |
797 | continue; | |
798 | #else | |
3e87ef69 | 799 | if(cluster->fY < 1 || cluster->fY > AliL3Transform::GetNPads(padrow) - 2 || |
800 | cluster->fX < 1 || cluster->fX > AliL3Transform::GetNTimeBins() - 2) | |
e1a32fa5 | 801 | continue; |
b2a02bce | 802 | #endif |
3e87ef69 | 803 | Float_t xyz_ex[3]; |
804 | ||
b2a02bce | 805 | #ifdef use_newio |
806 | AliL3Transform::Raw2Local(xyz_ex,cursec,currow,cluster->GetY(),cluster->GetX()); | |
807 | #else | |
494fad94 | 808 | AliL3Transform::Raw2Local(xyz_ex,cursec,currow,cluster->fY,cluster->fX); |
b2a02bce | 809 | #endif |
8af8a456 | 810 | //In function AliTPC::Hits2ExactClusters the time offset is not included, |
811 | //so we have to substract it again here. | |
3e87ef69 | 812 | if(slice<18) |
813 | xyz_ex[2]-=AliL3Transform::GetZOffset(); | |
814 | else | |
815 | xyz_ex[2]+=AliL3Transform::GetZOffset(); | |
8af8a456 | 816 | |
3e87ef69 | 817 | //Outside our cone: |
818 | if(param->GetPadRowRadii(cursec,currow)<230./250.*fabs(xyz_ex[2])) | |
819 | continue; | |
820 | ||
b2a02bce | 821 | TParticle *part = astack->Particle(mcId); |
3e87ef69 | 822 | crosscount++; |
823 | ||
824 | if(part->Pt() < fMinGoodPt) continue; | |
825 | ||
826 | //Dont take secondaries, because in width calculation we assume primaries: | |
827 | //if(!(part->TestBit(kPrimaryCharged))) continue; | |
828 | if(part->GetFirstMother()>=0) continue; | |
829 | ||
830 | Int_t tempcount=0; | |
831 | for(UInt_t c=0; c<fNcl[slice][patch]; c++) | |
832 | { | |
833 | if((Int_t)points[c].fPadRow!=padrow) continue; | |
834 | Float_t xyz_cl[3] = {points[c].fX,points[c].fY,points[c].fZ}; | |
835 | ||
836 | if(!offline) | |
837 | AliL3Transform::Global2Local(xyz_cl,cursec); | |
838 | tempcount++; | |
839 | ||
840 | if(points[c].fTrackID[0] != mcId && | |
841 | points[c].fTrackID[1] != mcId && | |
842 | points[c].fTrackID[2] != mcId) | |
843 | continue; | |
844 | ||
845 | //Residuals: | |
846 | Float_t resy = xyz_cl[1] - xyz_ex[1]; | |
847 | Float_t resz = xyz_cl[2] - xyz_ex[2]; | |
848 | ||
849 | //Cluster shape | |
850 | Int_t charge = (Int_t)points[c].fCharge; | |
851 | Float_t beta = GetCrossingAngle(part,slice,padrow,xyz_ex); | |
852 | Double_t tanl = xyz_ex[2]/sqrt(xyz_ex[0]*xyz_ex[0]+xyz_ex[1]*xyz_ex[1]); | |
853 | Float_t psigmaY2 = AliL3Transform::GetParSigmaY2(padrow,xyz_ex[2],beta); | |
854 | Float_t psigmaZ2 = AliL3Transform::GetParSigmaZ2(padrow,xyz_ex[2],tanl); | |
855 | Float_t sigmaY2 = points[c].fSigmaY2; | |
856 | Float_t sigmaZ2 = points[c].fSigmaZ2; | |
857 | ntuppel->Fill(slice,padrow,charge,resy,resz,xyz_ex[2],part->Pt(),beta,sigmaY2,sigmaZ2,psigmaY2,psigmaZ2); | |
858 | } | |
859 | clustercount=tempcount; | |
78127e35 | 860 | } |
3e87ef69 | 861 | ntuppel2->Fill(slice,padrow,clustercount,crosscount); |
862 | arr->ClearRow(cursec,currow); | |
863 | } | |
78127e35 | 864 | } |
3e87ef69 | 865 | exfile->Close(); |
866 | exact->Close(); | |
e1a32fa5 | 867 | |
3e87ef69 | 868 | TFile *ofile = TFile::Open(tofile,"RECREATE"); |
869 | ntuppel->Write(); | |
870 | ntuppel2->Write(); | |
871 | ofile->Close(); | |
872 | ||
5bfdc6ad | 873 | #endif |
78127e35 | 874 | } |
355debe1 | 875 | |
3e87ef69 | 876 | void AliL3Evaluate::GetCFeff(Char_t *path,Char_t *outfile,Int_t nevent,Bool_t sp) |
78127e35 | 877 | { |
3e87ef69 | 878 | //Evaluate the cluster finder efficiency. |
78127e35 | 879 | |
3e87ef69 | 880 | #ifndef do_mc |
881 | cerr<<"AliL3Evaluate::GetCFeff : Compile with do_mc flag"<<endl; | |
882 | return; | |
883 | #else | |
884 | TNtuple *ntuppel = new TNtuple("ntuppel","Cluster finder efficiency","slice:row:ncrossings:nclusters"); | |
885 | ntuppel->SetDirectory(0); | |
355debe1 | 886 | |
3e87ef69 | 887 | Char_t filename[1024]; |
888 | sprintf(filename,"%s/alirunfile.root",path); | |
889 | TFile *rfile = TFile::Open(filename); | |
890 | gAlice = (AliRun*)rfile->Get("gAlice"); | |
b2a02bce | 891 | |
892 | AliStack *astack=gAlice->Stack(); | |
355debe1 | 893 | |
3e87ef69 | 894 | AliTPCParam *param = (AliTPCParam*)rfile->Get(AliL3Transform::GetParamName()); |
895 | ||
896 | Int_t zero=param->GetZeroSup(); | |
897 | ||
898 | sprintf(filename,"%s/digitfile.root",path); | |
899 | TFile *dfile = TFile::Open(filename); | |
355debe1 | 900 | |
3e87ef69 | 901 | for(Int_t event=0; event<nevent; event++) |
78127e35 | 902 | { |
3e87ef69 | 903 | LoadData(event,sp); |
904 | rfile->cd(); | |
905 | gAlice->GetEvent(event); | |
b2a02bce | 906 | Int_t np = astack->GetNtrack(); |
3e87ef69 | 907 | cout<<"Processing event "<<event<<" with "<<np<<" particles "<<endl; |
908 | dfile->cd(); | |
909 | sprintf(filename,"TreeD_75x40_100x60_150x60_%d",event); | |
910 | TTree *TD=(TTree*)gDirectory->Get(filename); | |
911 | AliSimDigits da, *digits=&da; | |
912 | TD->GetBranch("Segment")->SetAddress(&digits); | |
913 | ||
914 | Int_t crossed=0,recs=0; | |
915 | Int_t *count = new Int_t[np]; //np number of particles. | |
916 | Int_t i; | |
917 | Float_t xyz[3]; | |
918 | for (i=0; i<np; i++) count[i]=0; | |
919 | ||
920 | ||
921 | Int_t sec,row,sl,sr; | |
922 | for(Int_t i=0; i<(Int_t)TD->GetEntries(); i++) | |
78127e35 | 923 | { |
3e87ef69 | 924 | crossed=recs=0; |
925 | if (!TD->GetEvent(i)) continue; | |
926 | param->AdjustSectorRow(digits->GetID(),sec,row); | |
927 | AliL3Transform::Sector2Slice(sl,sr,sec,row); | |
928 | if(sl < fMinSlice) continue; | |
929 | if(sl > fMaxSlice) break; | |
930 | cout<<"Processing slice "<<sl<<" row "<<sr<<endl; | |
931 | digits->First(); | |
355debe1 | 932 | do { |
3e87ef69 | 933 | Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn(); |
934 | Short_t dig = digits->GetDigit(it,ip); | |
355debe1 | 935 | |
3e87ef69 | 936 | if(dig<=param->GetZeroSup()) continue; |
494fad94 | 937 | AliL3Transform::Raw2Local(xyz,sec,row,ip,it); |
3e87ef69 | 938 | if(param->GetPadRowRadii(sec,row)<230./250.*fabs(xyz[2])) |
355debe1 | 939 | continue; |
940 | ||
3e87ef69 | 941 | Int_t idx0=digits->GetTrackID(it,ip,0); |
942 | Int_t idx1=digits->GetTrackID(it,ip,1); | |
943 | Int_t idx2=digits->GetTrackID(it,ip,2); | |
355debe1 | 944 | |
945 | if (idx0>=0 && dig>=zero) count[idx0]+=1; | |
946 | if (idx1>=0 && dig>=zero) count[idx1]+=1; | |
947 | if (idx2>=0 && dig>=zero) count[idx2]+=1; | |
3e87ef69 | 948 | } while (digits->Next()); |
355debe1 | 949 | for (Int_t j=0; j<np; j++) |
78127e35 | 950 | { |
b2a02bce | 951 | TParticle *part = astack->Particle(j); |
3e87ef69 | 952 | if(part->Pt() < fMinGoodPt) continue; |
953 | if(part->GetFirstMother() >= 0) continue; | |
355debe1 | 954 | if (count[j]>1) //at least two digits at this padrow |
494fad94 | 955 | { |
956 | crossed++; | |
957 | count[j]=0; | |
958 | } | |
355debe1 | 959 | } |
3e87ef69 | 960 | |
961 | Int_t patch = AliL3Transform::GetPatch(sr); | |
962 | if(sp==kTRUE) | |
963 | patch=0; | |
964 | AliL3SpacePointData *points = fClusters[sl][patch]; | |
965 | if(!points) | |
966 | continue; | |
967 | for(UInt_t k=0; k<fNcl[sl][patch]; k++) | |
355debe1 | 968 | { |
3e87ef69 | 969 | if(points[k].fPadRow!=sr) continue; |
355debe1 | 970 | recs++; |
78127e35 | 971 | } |
3e87ef69 | 972 | ntuppel->Fill(sl,sr,crossed,recs); |
78127e35 | 973 | } |
494fad94 | 974 | |
3e87ef69 | 975 | TD->Delete(); |
976 | delete[] count; | |
78127e35 | 977 | } |
355debe1 | 978 | TFile *file = TFile::Open(outfile,"RECREATE"); |
979 | ntuppel->Write(); | |
980 | file->Close(); | |
981 | ||
3e87ef69 | 982 | rfile->Close(); |
983 | dfile->Close(); | |
984 | #endif | |
985 | } | |
986 | ||
dd7d3870 | 987 | Float_t AliL3Evaluate::GetCrossingAngle(TParticle *part,Int_t slice,Int_t /*padrow*/,Float_t *xyz) |
3e87ef69 | 988 | { |
989 | //Calculate the padrow crossing angle of the particle | |
990 | ||
991 | Double_t kappa = AliL3Transform::GetBField()*AliL3Transform::GetBFact()/part->Pt(); | |
992 | ||
993 | Double_t radius = 1/fabs(kappa); | |
994 | if(part->GetPdgCode() > 0) kappa = -kappa; | |
995 | ||
996 | Float_t angl[1] = {part->Phi()}; | |
997 | ||
998 | AliL3Transform::Global2LocalAngle(angl,slice); | |
999 | ||
1000 | Double_t charge = -1.*kappa; | |
1001 | ||
1002 | Double_t trackPhi0 = angl[0] + charge*0.5*AliL3Transform::Pi()/fabs(charge); | |
1003 | ||
1004 | Double_t x0=0; | |
1005 | Double_t y0=0; | |
1006 | Double_t xc = x0 - radius * cos(trackPhi0); | |
1007 | Double_t yc = y0 - radius * sin(trackPhi0); | |
1008 | ||
1009 | Double_t tangent[2]; | |
1010 | tangent[0] = -1.*(xyz[1] - yc)/radius; | |
1011 | tangent[1] = (xyz[0] - xc)/radius; | |
1012 | ||
1013 | Double_t perp_padrow[2] = {1,0}; //locally in slice | |
1014 | ||
1015 | Double_t cos_beta = fabs(tangent[0]*perp_padrow[0] + tangent[1]*perp_padrow[1]); | |
1016 | if(cos_beta > 1) cos_beta=1; | |
1017 | return acos(cos_beta); | |
78127e35 | 1018 | } |
355debe1 | 1019 | |
3e87ef69 | 1020 | Int_t AliL3Evaluate::FindPrimaries(Int_t nparticles) |
1021 | { | |
1022 | // cuts: | |
1023 | Double_t vertcut = 0.001; | |
1024 | Double_t decacut = 3.; | |
1025 | Double_t timecut = 0.; | |
1026 | Int_t nprch1=0; | |
b2a02bce | 1027 | AliStack *astack=gAlice->Stack(); |
1028 | TParticle * part = astack->Particle(0); | |
3e87ef69 | 1029 | Double_t xori = part->Vx(); |
1030 | Double_t yori = part->Vy(); | |
1031 | Double_t zori = part->Vz(); | |
1032 | for(Int_t iprim = 0; iprim<nparticles; iprim++){ //loop on tracks | |
1033 | ||
b2a02bce | 1034 | part = astack->Particle(iprim); |
62bb4b3d | 1035 | const char * xxx=strstr(part->GetName(),"XXX"); |
3e87ef69 | 1036 | if(xxx)continue; |
1037 | ||
1038 | TParticlePDG *ppdg = part->GetPDG(); | |
1039 | if(TMath::Abs(ppdg->Charge())!=3)continue; // only charged (no quarks) | |
1040 | ||
1041 | Double_t dist=TMath::Sqrt((part->Vx()-xori)*(part->Vx()-xori)+(part->Vy()-yori)*(part->Vy()-yori)+(part->Vz()-zori)*(part->Vz()-zori)); | |
1042 | if(dist>vertcut)continue; // cut on the vertex | |
1043 | ||
1044 | if(part->T()>timecut)continue; | |
1045 | ||
1046 | Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz()); | |
1047 | if(ptot==(TMath::Abs(part->Pz())))continue; // no beam particles | |
1048 | ||
1049 | Bool_t prmch = kTRUE; // candidate primary track | |
1050 | Int_t fidau=part->GetFirstDaughter(); // cut on daughters | |
1051 | Int_t lasdau=0; | |
1052 | Int_t ndau=0; | |
1053 | if(fidau>=0){ | |
1054 | lasdau=part->GetLastDaughter(); | |
1055 | ndau=lasdau-fidau+1; | |
1056 | } | |
1057 | if(ndau>0){ | |
1058 | for(Int_t j=fidau;j<=lasdau;j++){ | |
b2a02bce | 1059 | TParticle *dau=astack->Particle(j); |
3e87ef69 | 1060 | Double_t distd=TMath::Sqrt((dau->Vx()-xori)*(dau->Vx()-xori)+(dau->Vy()-yori)*(dau->Vy()-yori)+(dau->Vz()-zori)*(dau->Vz()-zori)); |
1061 | if(distd<decacut)prmch=kFALSE; // eliminate if the decay is near the vertex | |
1062 | } | |
1063 | } | |
1064 | ||
1065 | if(prmch){ | |
1066 | nprch1++; | |
1067 | part->SetBit(kPrimaryCharged); | |
1068 | } | |
1069 | } | |
1070 | ||
1071 | return nprch1; | |
1072 | } |