]>
Commit | Line | Data |
---|---|---|
b661165c | 1 | // Author: Anders Vestbo <mailto:vestbo@fi.uib.no> |
2 | //*-- Copyright © ASV | |
108615fc | 3 | |
4 | #include <stdio.h> | |
e170e0f9 | 5 | #include <math.h> |
40896357 | 6 | #include <fstream.h> |
108615fc | 7 | #include <TFile.h> |
8 | #include <TH1.h> | |
9 | #include <TParticle.h> | |
10 | #include <TTree.h> | |
11 | #include <TClonesArray.h> | |
12 | ||
13 | #include "AliRun.h" | |
14 | #include "AliSimDigits.h" | |
15 | #include "AliTPC.h" | |
16 | #include "AliTPCClustersArray.h" | |
17 | #include "AliTPCClustersRow.h" | |
18 | #include "AliTPCcluster.h" | |
19 | #include "AliTPCParam.h" | |
e1a32fa5 | 20 | #include "AliComplexCluster.h" |
108615fc | 21 | |
22 | #include "AliL3Transform.h" | |
23 | #include "AliL3SpacePointData.h" | |
24 | #include "AliL3Track.h" | |
25 | #include "AliL3FileHandler.h" | |
26 | #include "AliL3TrackArray.h" | |
27 | #include "AliL3Evaluate.h" | |
28 | #include "AliL3Logging.h" | |
29 | ||
b661165c | 30 | //_____________________________________________________________ |
31 | // AliL3Evaluate | |
32 | // | |
33 | // Evaluation class for tracking. Plots efficiencies etc.. | |
84122e15 | 34 | // |
35 | ||
108615fc | 36 | |
37 | ClassImp(AliL3Evaluate) | |
38 | ||
39 | AliL3Evaluate::AliL3Evaluate() | |
40 | { | |
e1a32fa5 | 41 | fDigitsFile = NULL; |
108615fc | 42 | fTracks = NULL; |
43 | fMCclusterfile = NULL; | |
c0217812 | 44 | fNFastPoints = 0; |
45 | fMcIndex = 0; | |
46 | fMcId = 0; | |
1f79afc0 | 47 | fMinSlice=0; |
48 | fMaxSlice=0; | |
108615fc | 49 | } |
50 | ||
40896357 | 51 | AliL3Evaluate::AliL3Evaluate(Char_t *mcfile,Int_t *slice) |
108615fc | 52 | { |
53 | //Normal constructor. Input are the rootfile containing the | |
54 | //original MC information of the simulated event. | |
55 | ||
e1a32fa5 | 56 | fEventFile = new TFile(mcfile,"READ"); |
40896357 | 57 | |
e1a32fa5 | 58 | fParam = (AliTPCParam*)fEventFile->Get("75x40_100x60"); |
108615fc | 59 | |
60 | fMinSlice = slice[0]; | |
61 | fMaxSlice = slice[1]; | |
e170e0f9 | 62 | fMinGoodPt = 0.1; |
23908b9b | 63 | fNoOverlap = kFALSE; |
108615fc | 64 | } |
65 | ||
1f79afc0 | 66 | AliL3Evaluate::AliL3Evaluate(Int_t *slice) |
67 | { | |
78127e35 | 68 | //ctor to use if you do not need any rootfile. |
69 | ||
1f79afc0 | 70 | |
71 | fMinSlice = slice[0]; | |
72 | fMaxSlice = slice[1]; | |
1f79afc0 | 73 | |
74 | } | |
75 | ||
108615fc | 76 | AliL3Evaluate::~AliL3Evaluate() |
77 | { | |
c0217812 | 78 | if(fDigitsTree) fDigitsTree->Delete(); |
e1a32fa5 | 79 | if(fDigitsFile) { |
80 | fDigitsFile->Close(); | |
81 | delete fDigitsFile; | |
c0217812 | 82 | } |
c0217812 | 83 | if(fTracks) delete fTracks; |
84 | if(fPtRes) delete fPtRes; | |
85 | if(fNGoodTracksPt) delete fNGoodTracksPt; | |
86 | if(fNFoundTracksPt) delete fNFoundTracksPt; | |
87 | if(fNFakeTracksPt) delete fNFakeTracksPt; | |
88 | if(fTrackEffPt) delete fTrackEffPt; | |
89 | if(fFakeTrackEffPt) delete fFakeTrackEffPt; | |
90 | if(fNGoodTracksEta) delete fNGoodTracksEta; | |
91 | if(fNFoundTracksEta) delete fNFoundTracksEta; | |
92 | if(fNFakeTracksEta) delete fNFakeTracksEta; | |
93 | if(fTrackEffEta) delete fTrackEffEta; | |
94 | if(fFakeTrackEffEta) delete fFakeTrackEffEta; | |
95 | if(fMcIndex) delete [] fMcIndex; | |
96 | if(fMcId) delete [] fMcId; | |
1f79afc0 | 97 | if(fNtuppel) delete fNtuppel; |
108615fc | 98 | } |
99 | ||
bc2f4f0e | 100 | void AliL3Evaluate::Setup(Char_t *trackfile,Char_t *path) |
108615fc | 101 | { |
108615fc | 102 | //Read in the hit and track information from produced files. |
103 | ||
23908b9b | 104 | fGoodFound = 0; |
105 | fGoodGen = 0; | |
108615fc | 106 | Char_t fname[256]; |
40896357 | 107 | AliL3FileHandler *clusterfile[36][6]; |
108615fc | 108 | for(Int_t s=fMinSlice; s<=fMaxSlice; s++) |
109 | { | |
40896357 | 110 | for(Int_t p=0; p<6; p++) |
108615fc | 111 | { |
1f79afc0 | 112 | fClusters[s][p] = 0; |
108615fc | 113 | clusterfile[s][p] = new AliL3FileHandler(); |
bc2f4f0e | 114 | sprintf(fname,"%s/points_%d_%d.raw",path,s,p); |
108615fc | 115 | if(!clusterfile[s][p]->SetBinaryInput(fname)) |
116 | { | |
117 | LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open") | |
118 | <<"Inputfile "<<fname<<" does not exist"<<ENDLOG; | |
bc2f4f0e | 119 | delete clusterfile[s][p]; |
120 | clusterfile[s][p] = 0; | |
121 | continue; | |
108615fc | 122 | } |
123 | fClusters[s][p] = (AliL3SpacePointData*)clusterfile[s][p]->Allocate(); | |
124 | clusterfile[s][p]->Binary2Memory(fNcl[s][p],fClusters[s][p]); | |
125 | clusterfile[s][p]->CloseBinaryInput(); | |
126 | } | |
127 | } | |
78127e35 | 128 | |
23908b9b | 129 | |
108615fc | 130 | AliL3FileHandler *tfile = new AliL3FileHandler(); |
bc2f4f0e | 131 | if(!tfile->SetBinaryInput(trackfile)){ |
132 | LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open") | |
133 | <<"Inputfile "<<trackfile<<" does not exist"<<ENDLOG; | |
134 | return; | |
135 | } | |
108615fc | 136 | fTracks = new AliL3TrackArray(); |
137 | tfile->Binary2TrackArray(fTracks); | |
138 | tfile->CloseBinaryInput(); | |
139 | delete tfile; | |
bc2f4f0e | 140 | } |
141 | ||
40896357 | 142 | void AliL3Evaluate::SetupSlow(Char_t *trackfile,Char_t *digitsfile,Char_t *path) |
bc2f4f0e | 143 | { |
144 | //Setup for using the slow simulator. | |
145 | ||
40896357 | 146 | fDigitsFile = new TFile(digitsfile,"READ"); |
147 | if(!fDigitsFile->IsOpen()) | |
148 | { | |
149 | LOG(AliL3Log::kError,"AliL3Evaluate::SetupSlow","File Open") | |
150 | <<"Inputfile "<<digitsfile<<" does not exist"<<ENDLOG; | |
151 | return; | |
152 | } | |
bc2f4f0e | 153 | fIsSlow = true; |
154 | Setup(trackfile,path); | |
108615fc | 155 | |
e1a32fa5 | 156 | if(!InitMC()) |
bc2f4f0e | 157 | LOG(AliL3Log::kError,"AliL3Evaluation::SetupSlow","Digits Tree") |
158 | <<"Error setting up digits tree"<<ENDLOG; | |
e1a32fa5 | 159 | |
108615fc | 160 | } |
161 | ||
bc2f4f0e | 162 | void AliL3Evaluate::SetupFast(Char_t *trackfile,Char_t *mcClusterfile,Char_t *path) |
108615fc | 163 | { |
164 | //Setup for using the fast simulator. | |
165 | ||
166 | fIsSlow = false; | |
c0217812 | 167 | GetFastClusterIDs(path); |
108615fc | 168 | |
169 | fMCclusterfile = new TFile(mcClusterfile); | |
170 | if(!fMCclusterfile->IsOpen()) | |
171 | LOG(AliL3Log::kError,"AliL3Evaluation::SetupFast","File Open") | |
172 | <<"Inputfile "<<mcClusterfile<<" does not exist"<<ENDLOG; | |
173 | ||
bc2f4f0e | 174 | Setup(trackfile,path); |
e1a32fa5 | 175 | InitMC(); |
108615fc | 176 | } |
177 | ||
e1a32fa5 | 178 | Bool_t AliL3Evaluate::InitMC() |
108615fc | 179 | { |
e1a32fa5 | 180 | if(fIsSlow) |
108615fc | 181 | { |
e1a32fa5 | 182 | fDigitsFile->cd(); |
40896357 | 183 | fDigitsTree = (TTree*)fDigitsFile->Get("TreeD_75x40_100x60_0"); |
184 | if(!fDigitsTree) | |
185 | { | |
186 | LOG(AliL3Log::kError,"AliL3Evaluate::InitMC","Digits Tree") | |
187 | <<AliL3Log::kHex<<"Error getting digitstree "<<(Int_t)fDigitsTree<<ENDLOG; | |
188 | return false; | |
189 | } | |
e1a32fa5 | 190 | fDigitsTree->GetBranch("Segment")->SetAddress(&fDigits); |
191 | for(Int_t i=0; i<fDigitsTree->GetEntries(); i++) | |
192 | { | |
193 | if(!fDigitsTree->GetEvent(i)) continue; | |
194 | Int_t se,ro,slice,slicerow; | |
195 | fParam->AdjustSectorRow(fDigits->GetID(),se,ro); | |
494fad94 | 196 | AliL3Transform::Sector2Slice(slice,slicerow,se,ro); |
e1a32fa5 | 197 | fRowid[slice][slicerow] = i; |
198 | } | |
108615fc | 199 | } |
200 | ||
e1a32fa5 | 201 | fEventFile->cd(); |
40896357 | 202 | gAlice = (AliRun*)fEventFile->Get("gAlice"); |
108615fc | 203 | if (!gAlice) |
204 | { | |
40896357 | 205 | LOG(AliL3Log::kError,"AliL3Evaluate::InitMC","gAlice") |
108615fc | 206 | <<"AliRun object non existing on file"<<ENDLOG; |
207 | return false; | |
208 | } | |
209 | ||
210 | gAlice->GetEvent(0); | |
e1a32fa5 | 211 | |
108615fc | 212 | return true; |
e1a32fa5 | 213 | |
108615fc | 214 | } |
215 | ||
216 | ||
e1a32fa5 | 217 | |
40896357 | 218 | void AliL3Evaluate::DefineGoodTracks(Int_t *slice,Int_t *padrow,Int_t good_number,Char_t *fname) |
108615fc | 219 | { |
220 | //Loop over MC particles, and mark the good ones | |
221 | //(which the tracker should find...) | |
e1a32fa5 | 222 | |
40896357 | 223 | |
224 | ifstream in(fname); | |
225 | if(in) | |
226 | { | |
227 | LOG(AliL3Log::kInformational,"AliL3Evaluate::DefineGoodTracks","infile") | |
228 | <<"Reading good particles from file "<<fname<<ENDLOG; | |
229 | while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>> | |
230 | fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>> | |
231 | fGoodTracks[fGoodGen].pt) | |
232 | { | |
233 | fGoodGen++; | |
234 | if (fGoodGen==15000) | |
235 | { | |
236 | LOG(AliL3Log::kWarning,"AliL3Evaluate::DefineGoodTracks","fGoodGen") | |
237 | <<"Too many good tracks"<<ENDLOG; | |
238 | break; | |
239 | } | |
240 | } | |
241 | if (!in.eof()) cerr<<"Read error (good_tracks_tpc) !\n"; | |
242 | return; | |
243 | } | |
244 | ||
245 | ||
108615fc | 246 | AliTPC *TPC = (AliTPC*)gAlice->GetDetector("TPC"); |
e1a32fa5 | 247 | |
108615fc | 248 | TPC->SetParam(fParam); |
e1a32fa5 | 249 | |
108615fc | 250 | Int_t ver = TPC->IsVersion(); |
251 | LOG(AliL3Log::kInformational,"AliL3Evaluate::DefineGoodTracks","TPC version") | |
252 | <<"TPC version "<<ver<<" found on file"<<ENDLOG; | |
253 | ||
108615fc | 254 | Int_t zero=TPC->GetParam()->GetZeroSup(); |
255 | ||
e1a32fa5 | 256 | Int_t np = gAlice->GetNtrack(); |
108615fc | 257 | Int_t *good = new Int_t[np]; |
258 | for(Int_t ii=0; ii<np; ii++) | |
259 | good[ii] = 0; | |
260 | ||
261 | if(ver==1) | |
262 | { | |
263 | if(fIsSlow) | |
264 | LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","TPC version") | |
265 | <<"TPC version "<<ver<<" does not match."<<ENDLOG; | |
266 | fMCclusterfile->cd(); | |
267 | AliTPCClustersArray carray; | |
268 | carray.Setup(fParam); | |
269 | carray.SetClusterType("AliTPCcluster"); | |
40896357 | 270 | Char_t cname[100]; |
271 | Int_t eventn = 0; | |
272 | sprintf(cname,"TreeC_TPC_%d",eventn); | |
273 | Bool_t clusterok = carray.ConnectTree(cname); | |
108615fc | 274 | if(!clusterok) |
275 | LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","Cluster Array") | |
276 | <<"Error loading clusters from rootfile"<<ENDLOG; | |
78127e35 | 277 | |
bc2f4f0e | 278 | for(Int_t i=0; i<carray.GetTree()->GetEntries(); i++) |
108615fc | 279 | { |
280 | Int_t sec,row,sl,lr; | |
281 | AliSegmentID *s = carray.LoadEntry(i); | |
282 | fParam->AdjustSectorRow(s->GetID(),sec,row); | |
494fad94 | 283 | AliL3Transform::Sector2Slice(sl,lr,sec,row); |
108615fc | 284 | |
e1a32fa5 | 285 | if(sl != slice[0]) {carray.ClearRow(sec,row); continue;} |
bc2f4f0e | 286 | if(lr < padrow[0]) {carray.ClearRow(sec,row); continue;} |
287 | if(lr > padrow[1]) {carray.ClearRow(sec,row); continue;} | |
108615fc | 288 | AliTPCClustersRow *cRow = carray.GetRow(sec,row); |
289 | for(Int_t j=0; j<cRow->GetArray()->GetEntriesFast(); j++) | |
290 | { | |
291 | AliTPCcluster *cluster=(AliTPCcluster*)(*cRow)[j]; | |
292 | Int_t lab=cluster->GetLabel(0); | |
293 | if(lab<0) continue; | |
294 | lab=TMath::Abs(lab); | |
295 | good[lab]++; | |
296 | } | |
297 | if(carray.GetRow(sec,row)) | |
298 | carray.ClearRow(sec,row); | |
299 | } | |
300 | } | |
301 | else if(ver==2) | |
302 | { | |
303 | if(!fIsSlow) | |
304 | LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","TPC version") | |
305 | <<"TPC version "<<ver<<" does not match."<<ENDLOG; | |
306 | Int_t *count = new Int_t[np]; //np number of particles. | |
307 | Int_t i; | |
355debe1 | 308 | Float_t xyz[3]; |
108615fc | 309 | for (i=0; i<np; i++) count[i]=0; |
e1a32fa5 | 310 | for(Int_t sl=slice[0]; sl<=slice[1]; sl++) |
311 | { | |
312 | for (i=padrow[0]; i<=padrow[1]; i++) { | |
313 | Int_t index = fRowid[sl][i]; | |
314 | if (!fDigitsTree->GetEvent(index)) continue; | |
315 | Int_t sec,row; | |
316 | fParam->AdjustSectorRow(fDigits->GetID(),sec,row); | |
317 | fDigits->First(); | |
318 | do { | |
319 | Int_t it=fDigits->CurrentRow(), ip=fDigits->CurrentColumn(); | |
320 | Short_t dig = fDigits->GetDigit(it,ip); | |
355debe1 | 321 | |
322 | if(dig<=fParam->GetZeroSup()) continue; | |
323 | if(it < fParam->GetMaxTBin()-1 && it > 0) | |
324 | if(fDigits->GetDigit(it+1,ip) <= fParam->GetZeroSup() | |
325 | && fDigits->GetDigit(it-1,ip) <= fParam->GetZeroSup()) | |
326 | continue; | |
327 | ||
494fad94 | 328 | AliL3Transform::Raw2Local(xyz,sec,row,ip,it); |
355debe1 | 329 | if(fParam->GetPadRowRadii(sec,row)<230./250.*fabs(xyz[2])) |
330 | continue; | |
331 | ||
e1a32fa5 | 332 | Int_t idx0=fDigits->GetTrackID(it,ip,0); |
333 | Int_t idx1=fDigits->GetTrackID(it,ip,1); | |
334 | Int_t idx2=fDigits->GetTrackID(it,ip,2); | |
e170e0f9 | 335 | |
e1a32fa5 | 336 | if (idx0>=0 && dig>=zero) count[idx0]+=1; |
337 | if (idx1>=0 && dig>=zero) count[idx1]+=1; | |
338 | if (idx2>=0 && dig>=zero) count[idx2]+=1; | |
339 | } while (fDigits->Next()); | |
340 | ||
341 | for (Int_t j=0; j<np; j++) | |
342 | { | |
343 | if (count[j]>1) //at least two digits at this padrow | |
344 | good[j]++; | |
345 | ||
346 | count[j]=0; | |
347 | } | |
108615fc | 348 | } |
108615fc | 349 | } |
108615fc | 350 | delete[] count; |
351 | } | |
e1a32fa5 | 352 | |
108615fc | 353 | else |
354 | { | |
355 | LOG(AliL3Log::kError,"AliL3Evaluation::FillEffHistos","TPC version") | |
356 | <<"No valid TPC version found"<<ENDLOG; | |
e1a32fa5 | 357 | return; |
108615fc | 358 | } |
359 | ||
e1a32fa5 | 360 | |
361 | //The following code has been taken from offlinemacro->AliTPCComparison.C | |
362 | ||
363 | TTree *TH=gAlice->TreeH(); | |
364 | Int_t npart=(Int_t)TH->GetEntries(); | |
e1a32fa5 | 365 | Int_t max = 15000; |
366 | while (npart--) { | |
367 | AliTPChit *hit0=0; | |
368 | ||
40896357 | 369 | TPC->ResetHits(); |
370 | TH->GetEvent(npart); | |
371 | AliTPChit * hit = (AliTPChit*) TPC->FirstHit(-1); | |
372 | while (hit){ | |
373 | if (hit->fQ==0.) break; | |
374 | hit = (AliTPChit*) TPC->NextHit(); | |
375 | } | |
376 | if (hit) { | |
377 | hit0 = new AliTPChit(*hit); //Make copy of hit | |
378 | hit = hit0; | |
379 | } | |
380 | else continue; | |
381 | AliTPChit *hit1=(AliTPChit*)TPC->NextHit(); | |
382 | if (hit1==0) continue; | |
383 | if (hit1->fQ != 0.) continue; | |
384 | Int_t i=hit->Track(); | |
385 | TParticle *p = (TParticle*)gAlice->Particle(i); | |
386 | ||
387 | //printf("Checking particle %d with code %d\n",i,p->GetPdgCode()); | |
388 | if (p->GetFirstMother()>=0) continue; //secondary particle | |
355debe1 | 389 | if (good[i] < good_number || good[i]>176) continue; |
40896357 | 390 | if (p->Pt()<fMinGoodPt) continue; |
391 | if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue; | |
392 | printf("Found good particle %d, nHits %d\n",i,good[i]); | |
393 | ||
394 | fGoodTracks[fGoodGen].label=i; | |
395 | fGoodTracks[fGoodGen].code=p->GetPdgCode(); | |
396 | //**** px py pz - in global coordinate system, x y z - in local ! | |
397 | fGoodTracks[fGoodGen].px=hit->X(); fGoodTracks[fGoodGen].py=hit->Y(); fGoodTracks[fGoodGen].pz=hit->Z(); | |
398 | fGoodTracks[fGoodGen].pt = p->Pt(); | |
399 | fGoodGen++; | |
400 | //nt++; | |
401 | if (hit0) delete hit0; | |
402 | if (fGoodGen==max) {cerr<<"Too many good tracks !n"; break;} | |
403 | } | |
e170e0f9 | 404 | |
108615fc | 405 | delete [] good; |
40896357 | 406 | |
407 | LOG(AliL3Log::kInformational,"AliL3Evaluate::DefineGoodParticles","Eval") | |
408 | <<AliL3Log::kDec<<"Found "<<fGoodGen<<"good tracks"<<ENDLOG; | |
409 | ||
410 | ofstream out(fname); | |
411 | if(out) | |
412 | { | |
413 | for (Int_t ngd=0; ngd<fGoodGen; ngd++) | |
414 | out<<fGoodTracks[ngd].label<<' '<<fGoodTracks[ngd].code<<' '<< | |
415 | fGoodTracks[ngd].px<<' '<<fGoodTracks[ngd].py<<' '<<fGoodTracks[ngd].pz<<' '<< | |
416 | fGoodTracks[ngd].pt<<endl; | |
417 | ||
418 | } | |
419 | else | |
420 | LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodParticles","file") | |
421 | <<"Can not open file containing good tracks"<<ENDLOG; | |
422 | ||
423 | out.close(); | |
424 | ||
108615fc | 425 | } |
426 | ||
427 | void AliL3Evaluate::EvaluatePatch(Int_t slice,Int_t patch,Int_t min_points,Int_t good_number) | |
428 | { | |
429 | //Make efficiency plots for tracking on patch level (before any merging). | |
430 | ||
40896357 | 431 | Int_t row[6][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,175}}; |
e1a32fa5 | 432 | Int_t sl[2] ={slice,slice}; |
433 | DefineGoodTracks(sl,row[patch],good_number); | |
108615fc | 434 | SetMinPoints(min_points); |
435 | AssignIDs(); | |
bc2f4f0e | 436 | CreateHistos(); |
e1a32fa5 | 437 | FillEffHistos(); |
bc2f4f0e | 438 | CalcEffHistos(); |
108615fc | 439 | } |
440 | ||
441 | void AliL3Evaluate::EvaluateSlice(Int_t slice,Int_t min_points,Int_t good_number) | |
442 | { | |
443 | //Make efficiency plots for tracking on a slice (after merging). | |
444 | //min_points = minimum points on track to be considered for evaluation | |
445 | //good_number = minimum hits (padrows) produced by simulated track for consideration. | |
446 | ||
ae97a0b9 | 447 | Int_t row[2] = {0,175}; |
e1a32fa5 | 448 | Int_t sl[2] ={slice,slice}; |
449 | DefineGoodTracks(sl,row,good_number); | |
23908b9b | 450 | |
108615fc | 451 | SetMinPoints(min_points); |
23908b9b | 452 | |
108615fc | 453 | AssignIDs(); |
bc2f4f0e | 454 | CreateHistos(); |
e1a32fa5 | 455 | FillEffHistos(); |
bc2f4f0e | 456 | CalcEffHistos(); |
108615fc | 457 | } |
458 | ||
40896357 | 459 | void AliL3Evaluate::EvaluateGlobal(Int_t min_points,Int_t good_number,Char_t *fname) |
108615fc | 460 | { |
461 | //Make efficiency plots for tracking on several slices. | |
bc2f4f0e | 462 | |
ae97a0b9 | 463 | Int_t row[2] = {0,175}; |
bc2f4f0e | 464 | SetMinPoints(min_points); |
e1a32fa5 | 465 | Int_t slice[2] = {fMinSlice,fMaxSlice}; |
40896357 | 466 | DefineGoodTracks(slice,row,good_number,fname); |
bc2f4f0e | 467 | AssignIDs(); |
355debe1 | 468 | CreateHistos(5,0,2); |
40896357 | 469 | printf("filling histos\n"); |
e1a32fa5 | 470 | FillEffHistos(); |
40896357 | 471 | printf("done fillling\n"); |
bc2f4f0e | 472 | CalcEffHistos(); |
108615fc | 473 | } |
474 | ||
108615fc | 475 | void AliL3Evaluate::AssignIDs() |
476 | { | |
477 | //Assign MC id to the tracks. | |
478 | ||
108615fc | 479 | fTracks->QSort(); |
480 | LOG(AliL3Log::kDebug,"AliL3Evaluate::AssignIDs","Track Loop") | |
481 | <<"Assigning MC id to the found tracks...."<<ENDLOG; | |
482 | for(Int_t i=0; i<fTracks->GetNTracks(); i++) | |
483 | { | |
484 | AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i); | |
c0217812 | 485 | if(!track) continue; |
108615fc | 486 | if(track->GetNumberOfPoints() < fMinPointsOnTrack) break; |
23908b9b | 487 | |
23908b9b | 488 | fGoodFound++; |
c0217812 | 489 | Int_t tID = GetMCTrackLabel(track); |
108615fc | 490 | track->SetMCid(tID); |
e170e0f9 | 491 | printf("track %i id %d nHits %d\n",i,tID,track->GetNumberOfPoints()); |
108615fc | 492 | } |
108615fc | 493 | } |
494 | ||
495 | ||
496 | struct S {Int_t lab; Int_t max;}; | |
c0217812 | 497 | Int_t AliL3Evaluate::GetMCTrackLabel(AliL3Track *track){ |
108615fc | 498 | //Returns the MCtrackID of the belonging clusters. |
499 | //If MCLabel < 0, means that track is fake. | |
500 | //Fake track means that more than 10 percent of clusters are assigned incorrectly. | |
494fad94 | 501 | |
502 | #ifdef do_mc | |
108615fc | 503 | Int_t num_of_clusters = track->GetNumberOfPoints(); |
504 | S *s=new S[num_of_clusters]; | |
505 | Int_t i; | |
506 | for (i=0; i<num_of_clusters; i++) s[i].lab=s[i].max=0; | |
507 | UInt_t *hitnum = track->GetHitNumbers(); | |
508 | UInt_t id; | |
108615fc | 509 | |
355debe1 | 510 | //Int_t **trackID = GetClusterIDs(track); |
108615fc | 511 | |
512 | Int_t lab=123456789; | |
513 | for (i=0; i<num_of_clusters; i++) | |
514 | { | |
515 | //Tricks to get the clusters belonging to this track: | |
516 | id = hitnum[i]; | |
517 | Int_t slice = (id>>25) & 0x7f; | |
518 | Int_t patch = (id>>22) & 0x7; | |
519 | UInt_t pos = id&0x3fffff; | |
520 | ||
521 | AliL3SpacePointData *points = fClusters[slice][patch]; | |
c0217812 | 522 | if(!points) continue; |
108615fc | 523 | if(pos>=fNcl[slice][patch]) |
524 | { | |
525 | LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray") | |
526 | <<AliL3Log::kDec<<"ERROR"<<ENDLOG; | |
527 | continue; | |
528 | } | |
529 | ||
530 | //Get the label of the cluster: | |
355debe1 | 531 | |
532 | lab=abs(points[pos].fTrackID[0]); | |
108615fc | 533 | Int_t j; |
534 | for (j=0; j<num_of_clusters; j++) | |
535 | if (s[j].lab==lab || s[j].max==0) break; | |
536 | s[j].lab=lab; | |
537 | s[j].max++; | |
538 | } | |
539 | ||
540 | Int_t max=0; | |
541 | for (i=0; i<num_of_clusters; i++) | |
542 | if (s[i].max>max) {max=s[i].max; lab=s[i].lab;} | |
543 | ||
544 | delete[] s; | |
545 | ||
546 | for (i=0; i<num_of_clusters; i++) | |
547 | { | |
548 | id = hitnum[i]; | |
549 | Int_t slice = (id>>25) & 0x7f; | |
550 | Int_t patch = (id>>22) & 0x7; | |
551 | UInt_t pos = id&0x3fffff; | |
552 | ||
553 | AliL3SpacePointData *points = fClusters[slice][patch]; | |
c0217812 | 554 | if(!points) continue; |
108615fc | 555 | if(pos>=fNcl[slice][patch]) |
556 | { | |
557 | LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray") | |
558 | <<AliL3Log::kDec<<"ERROR"<<ENDLOG; | |
559 | continue; | |
560 | } | |
561 | ||
355debe1 | 562 | if (abs(points[pos].fTrackID[1]) == lab || |
563 | abs(points[pos].fTrackID[2]) == lab ) max++; | |
108615fc | 564 | } |
565 | ||
566 | //check if more than 10% of the clusters are incorrectly assigned (fake track): | |
567 | if (1.-Float_t(max)/num_of_clusters > 0.10) | |
568 | { | |
569 | return -lab; | |
570 | } | |
571 | ||
108615fc | 572 | return lab; |
494fad94 | 573 | #else //If we are running with mc_ids or not |
574 | return 0; | |
575 | #endif | |
576 | ||
108615fc | 577 | } |
578 | ||
579 | ||
c0217812 | 580 | Int_t **AliL3Evaluate::GetClusterIDs(AliL3Track *track) |
108615fc | 581 | { |
582 | //Return the MC information of all clusters belonging to track. | |
c0217812 | 583 | |
108615fc | 584 | Int_t num_of_clusters = track->GetNumberOfPoints(); |
585 | Int_t **trackID = new Int_t*[num_of_clusters]; | |
586 | ||
587 | UInt_t *hitnum = track->GetHitNumbers(); | |
588 | UInt_t id; | |
108615fc | 589 | |
590 | Float_t xyz[3]; | |
c0217812 | 591 | Int_t padrow; |
108615fc | 592 | for(Int_t i=0; i<num_of_clusters; i++) |
593 | { | |
594 | id = hitnum[i]; | |
595 | Int_t slice = (id>>25) & 0x7f; | |
596 | Int_t patch = (id>>22) & 0x7; | |
597 | UInt_t pos = id&0x3fffff; | |
598 | ||
599 | AliL3SpacePointData *points = fClusters[slice][patch]; | |
600 | ||
601 | if(!points) | |
602 | { | |
603 | LOG(AliL3Log::kError,"AliL3Evaluate::GetClusterIDs","Clusterarray") | |
604 | <<"No points at slice "<<slice<<" patch "<<patch<<" pos "<<pos<<ENDLOG; | |
605 | continue; | |
606 | } | |
607 | if(pos>=fNcl[slice][patch]) | |
608 | { | |
609 | LOG(AliL3Log::kError,"AliL3Evaluate::GetClusterIDs","Clusterarray") | |
610 | <<AliL3Log::kDec<<"ERROR"<<ENDLOG; | |
611 | continue; | |
612 | } | |
613 | ||
614 | xyz[0] = points[pos].fX; | |
615 | xyz[1] = points[pos].fY; | |
616 | xyz[2] = points[pos].fZ; | |
617 | //sector = points[pos].fSector; | |
618 | padrow = points[pos].fPadRow; | |
619 | Int_t se,ro; | |
494fad94 | 620 | AliL3Transform::Slice2Sector(slice,padrow,se,ro); |
621 | AliL3Transform::Global2Raw(xyz,se,ro); | |
108615fc | 622 | |
623 | if(fIsSlow) | |
624 | { | |
625 | Int_t p = fRowid[slice][padrow]; | |
108615fc | 626 | if(!fDigitsTree->GetEvent(p)) |
627 | LOG(AliL3Log::kError,"AliL3Evaluate::GetClusterIDs","Digits Tree") | |
628 | <<"Error reading digits tree"<<ENDLOG; | |
629 | ||
630 | trackID[i] = new Int_t[3]; | |
e170e0f9 | 631 | trackID[i][0] = fDigits->GetTrackID((Int_t)rint(xyz[2]),(Int_t)rint(xyz[1]),0); |
632 | trackID[i][1] = fDigits->GetTrackID((Int_t)rint(xyz[2]),(Int_t)rint(xyz[1]),1); | |
633 | trackID[i][2] = fDigits->GetTrackID((Int_t)rint(xyz[2]),(Int_t)rint(xyz[1]),2); | |
634 | //if(trackID[i][0]==6 || trackID[i][0]==32) | |
635 | //printf("trackID %d, padrow %d pad %d time %d\n",fDigits->GetTrackID((Int_t)rint(xyz[2]),(Int_t)rint(xyz[1]),0),padrow,(int)rint(xyz[1]),(int)rint(xyz[2])); | |
636 | /*if(trackID[i][0]<0) | |
637 | { | |
638 | ||
639 | printf("trackID %d, padrow %d pad %d time %d\n",trackID[i][0],padrow,(int)rint(xyz[1]),(int)rint(xyz[2])); | |
640 | printf("on the side %d %d %d %d\n",fDigits->GetTrackID(((int)rint(xyz[2])-1),((int)rint(xyz[1])),0),fDigits->GetTrackID(((int)rint(xyz[2])+1),((int)rint(xyz[1])),0),fDigits->GetTrackID(((int)rint(xyz[2])),((int)rint(xyz[1])-1),0),fDigits->GetTrackID(((int)rint(xyz[2])),((int)rint(xyz[1])+1),0)); | |
641 | }*/ | |
108615fc | 642 | } |
643 | else | |
644 | { | |
645 | Int_t tmp_pid=0; | |
c0217812 | 646 | for(Int_t ii=0; ii<fNFastPoints; ii++) |
108615fc | 647 | { |
c0217812 | 648 | tmp_pid = fMcId[ii]; |
649 | if(fMcIndex[ii] == id) break; | |
108615fc | 650 | } |
651 | trackID[i] = new Int_t[3]; | |
652 | trackID[i][0] = tmp_pid; | |
653 | trackID[i][1] = -1; | |
654 | trackID[i][2] = -1; | |
655 | } | |
656 | } | |
e170e0f9 | 657 | |
108615fc | 658 | return trackID; |
659 | } | |
660 | ||
c0217812 | 661 | void AliL3Evaluate::GetFastClusterIDs(Char_t *path) |
108615fc | 662 | { |
663 | //Get the MC id of space points in case of using the fast simulator. | |
c0217812 | 664 | char fname[256]; |
665 | sprintf(fname,"%s/point_mc.dat",path); | |
666 | FILE *infile = fopen(fname,"r"); | |
667 | if(!infile) return; | |
108615fc | 668 | Int_t hitid,hitmc,i; |
669 | ||
670 | for(i=0; ; i++) | |
671 | if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break; | |
108615fc | 672 | rewind(infile); |
c0217812 | 673 | fNFastPoints = i; |
674 | fMcId = new Int_t[fNFastPoints]; | |
675 | fMcIndex = new UInt_t[fNFastPoints]; | |
108615fc | 676 | |
c0217812 | 677 | for(i=0; i<fNFastPoints; i++) |
108615fc | 678 | { |
679 | if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break; | |
c0217812 | 680 | fMcId[i] = hitmc; |
681 | fMcIndex[i] = hitid; | |
108615fc | 682 | } |
683 | fclose(infile); | |
c0217812 | 684 | } |
685 | ||
686 | void AliL3Evaluate::CreateHistos(Int_t nbin,Int_t xlow,Int_t xup) | |
687 | { | |
688 | //Create the histograms | |
689 | ||
40896357 | 690 | LOG(AliL3Log::kInformational,"AliL3Evaluate::CreateHistos","Allocating") |
691 | <<"Creating histograms..."<<ENDLOG; | |
692 | ||
1f79afc0 | 693 | fNtuppel = new TNtuple("fNtuppel","Pt resolution","pt_gen:pt_found:nHits"); |
694 | ||
c0217812 | 695 | fPtRes = new TH1F("fPtRes","Relative Pt resolution",30,-10.,10.); |
696 | fNGoodTracksPt = new TH1F("fNGoodTracksPt","Good tracks vs pt",nbin,xlow,xup); | |
697 | fNFoundTracksPt = new TH1F("fNFoundTracksPt","Found tracks vs pt",nbin,xlow,xup); | |
698 | fNFakeTracksPt = new TH1F("fNFakeTracksPt","Fake tracks vs pt",nbin,xlow,xup); | |
699 | fTrackEffPt = new TH1F("fTrackEffPt","Tracking efficiency vs pt",nbin,xlow,xup); | |
700 | fFakeTrackEffPt = new TH1F("fFakeTrackEffPt","Efficiency for fake tracks vs pt",nbin,xlow,xup); | |
701 | ||
702 | fNGoodTracksEta = new TH1F("fNGoodTracksEta","Good tracks vs eta",20,-50,50); | |
703 | fNFoundTracksEta = new TH1F("fNFoundTracksEta","Found tracks vs eta",20,-50,50); | |
704 | fNFakeTracksEta = new TH1F("fNFakeTracksEta","Fake tracks vs eta",20,-50,50); | |
705 | fTrackEffEta = new TH1F("fTrackEffEta","Tracking efficienct vs eta",20,-50,50); | |
706 | fFakeTrackEffEta = new TH1F("fFakeTrackEffEta","Efficiency for fake tracks vs eta",20,-50,50); | |
40896357 | 707 | printf("finished creating histos\n"); |
c0217812 | 708 | } |
355debe1 | 709 | /* |
e1a32fa5 | 710 | void AliL3Evaluate::FillEffHistos() |
c0217812 | 711 | { |
712 | //Fill the efficiency histograms. | |
355debe1 | 713 | |
e1a32fa5 | 714 | for(Int_t i=0; i<fGoodGen; i++) |
c0217812 | 715 | { |
e1a32fa5 | 716 | Double_t ptg=fGoodTracks[i].pt,pzg=fGoodTracks[i].pz; |
c0217812 | 717 | Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi(); |
355debe1 | 718 | //printf("filling particle with pt %f and dipangle %f\n",ptg,dipangle); |
c0217812 | 719 | fNGoodTracksPt->Fill(ptg); |
720 | fNGoodTracksEta->Fill(dipangle); | |
721 | Int_t found = 0; | |
722 | for(Int_t k=0; k<fTracks->GetNTracks(); k++) | |
723 | { | |
724 | AliL3Track *track = fTracks->GetCheckedTrack(k); | |
725 | if(!track) continue; | |
726 | Int_t nHits = track->GetNumberOfPoints(); | |
727 | if(nHits < fMinPointsOnTrack) break; | |
728 | ||
729 | Int_t tracklabel; | |
730 | tracklabel = track->GetMCid(); | |
355debe1 | 731 | //printf("evaluating track id %d\n",tracklabel); |
e1a32fa5 | 732 | if(TMath::Abs(tracklabel) != fGoodTracks[i].label) continue; |
c0217812 | 733 | found=1; |
e1a32fa5 | 734 | if(tracklabel == fGoodTracks[i].label) {fNFoundTracksPt->Fill(ptg); fNFoundTracksEta->Fill(dipangle);} |
c0217812 | 735 | else {fNFakeTracksPt->Fill(ptg); fNFakeTracksEta->Fill(dipangle);} |
736 | Float_t pt=track->GetPt(); | |
737 | fPtRes->Fill((pt-ptg)/ptg*100.); | |
1f79afc0 | 738 | fNtuppel->Fill(ptg,pt,nHits); |
c0217812 | 739 | break; |
740 | ||
741 | } | |
742 | } | |
743 | } | |
355debe1 | 744 | */ |
745 | void AliL3Evaluate::FillEffHistos() | |
746 | { | |
747 | //Fill the efficiency histograms. | |
748 | ||
749 | cout<<endl<<"Note: Doing NAIVE evaluation "<<endl; | |
750 | for(Int_t i=0; i<fGoodGen; i++) | |
751 | { | |
752 | Double_t ptg=fGoodTracks[i].pt,pzg=fGoodTracks[i].pz; | |
753 | Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi(); | |
754 | //printf("filling particle with pt %f and dipangle %f\n",ptg,dipangle); | |
755 | fNGoodTracksPt->Fill(ptg); | |
756 | fNGoodTracksEta->Fill(dipangle); | |
757 | ||
758 | } | |
759 | ||
760 | for(Int_t k=0; k<fTracks->GetNTracks(); k++) | |
761 | { | |
762 | AliL3Track *track = fTracks->GetCheckedTrack(k); | |
763 | if(!track) continue; | |
764 | Int_t nHits = track->GetNumberOfPoints(); | |
765 | if(nHits < fMinPointsOnTrack) break; | |
766 | if(track->GetPt()<fMinGoodPt) continue; | |
767 | if(fabs(track->GetPseudoRapidity())>0.9) continue; | |
768 | ||
769 | fNFoundTracksPt->Fill(track->GetPt()); fNFoundTracksEta->Fill(track->GetPseudoRapidity()); | |
770 | //Float_t pt=track->GetPt(); | |
771 | //fPtRes->Fill((pt-ptg)/ptg*100.); | |
772 | //fNtuppel->Fill(ptg,pt,nHits); | |
773 | ||
774 | } | |
775 | } | |
c0217812 | 776 | |
777 | void AliL3Evaluate::CalcEffHistos(){ | |
778 | ||
779 | Stat_t ngood=fNGoodTracksPt->GetEntries(); | |
780 | Stat_t nfound=fNFoundTracksPt->GetEntries(); | |
781 | Stat_t nfake=fNFakeTracksPt->GetEntries(); | |
782 | ||
783 | LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency") | |
784 | <<AliL3Log::kDec<<"There was "<<ngood<<" generated good tracks"<<ENDLOG; | |
785 | LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency") | |
786 | <<AliL3Log::kDec<<"Found "<<nfound<<" tracks"<<ENDLOG; | |
787 | LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency") | |
788 | <<AliL3Log::kDec<<"Integral efficiency is about "<<nfound/ngood*100<<ENDLOG; | |
789 | LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency") | |
790 | <<AliL3Log::kDec<<"Fake tracks relative is about "<<nfake/ngood*100<<ENDLOG; | |
791 | ||
355debe1 | 792 | LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency") |
793 | <<AliL3Log::kDec<<"Naive efficiency "<<(Double_t)fGoodFound/(Double_t)fGoodGen<<ENDLOG; | |
794 | ||
795 | ||
c0217812 | 796 | fNFoundTracksPt->Sumw2(); fNGoodTracksPt->Sumw2(); |
797 | fTrackEffPt->Divide(fNFoundTracksPt,fNGoodTracksPt,1,1.,"b"); | |
798 | fFakeTrackEffPt->Divide(fNFakeTracksPt,fNGoodTracksPt,1,1.,"b"); | |
799 | fTrackEffPt->SetMaximum(1.4); | |
800 | fTrackEffPt->SetXTitle("P_{T} [GeV]"); | |
801 | fTrackEffPt->SetLineWidth(2); | |
802 | fFakeTrackEffPt->SetFillStyle(3013); | |
803 | fTrackEffPt->SetLineColor(4); | |
804 | fFakeTrackEffPt->SetFillColor(2); | |
805 | ||
806 | fNFoundTracksEta->Sumw2(); fNGoodTracksEta->Sumw2(); | |
807 | fTrackEffEta->Divide(fNFoundTracksEta,fNGoodTracksEta,1,1.,"b"); | |
808 | fFakeTrackEffEta->Divide(fNFakeTracksEta,fNGoodTracksEta,1,1.,"b"); | |
809 | fTrackEffEta->SetMaximum(1.4); | |
810 | fTrackEffEta->SetXTitle("#lambda [degrees]"); | |
811 | fTrackEffEta->SetLineWidth(2); | |
812 | fFakeTrackEffEta->SetFillStyle(3013); | |
813 | fTrackEffEta->SetLineColor(4); | |
814 | fFakeTrackEffEta->SetFillColor(2); | |
815 | ||
108615fc | 816 | } |
817 | ||
818 | void AliL3Evaluate::Write2File(Char_t *outputfile) | |
819 | { | |
820 | //Write histograms to file: | |
821 | ||
822 | TFile *of = new TFile(outputfile,"RECREATE"); | |
823 | if(!of->IsOpen()) | |
824 | { | |
825 | LOG(AliL3Log::kError,"AliL3Evaluate::Write2File","File Open") | |
826 | <<"Problems opening rootfile"<<ENDLOG; | |
827 | return; | |
828 | } | |
829 | ||
830 | of->cd(); | |
1f79afc0 | 831 | fNtuppel->Write(); |
108615fc | 832 | fPtRes->Write(); |
833 | fNGoodTracksPt->Write(); | |
834 | fNFoundTracksPt->Write(); | |
835 | fNFakeTracksPt->Write(); | |
836 | fTrackEffPt->Write(); | |
837 | fFakeTrackEffPt->Write(); | |
838 | fNGoodTracksEta->Write(); | |
839 | fNFoundTracksEta->Write(); | |
840 | fNFakeTracksEta->Write(); | |
841 | fTrackEffEta->Write(); | |
842 | fFakeTrackEffEta->Write(); | |
843 | ||
844 | of->Close(); | |
845 | delete of; | |
846 | ||
847 | } | |
848 | ||
1f79afc0 | 849 | TNtuple *AliL3Evaluate::CalculateResiduals() |
850 | { | |
851 | ||
852 | TNtuple *ntuppel=new TNtuple("ntuppel","Residuals","residual_trans:residual_long:zHit:pt:dipangle:beta:padrow:nHits"); | |
78127e35 | 853 | |
1f79afc0 | 854 | for(int f=fMinSlice; f<=fMaxSlice; f++) |
855 | { | |
856 | AliL3FileHandler *tfile = new AliL3FileHandler(); | |
857 | char fname[256]; | |
858 | sprintf(fname,"tracks_tr_%d_0.raw",f); | |
859 | if(!tfile->SetBinaryInput(fname)){ | |
860 | LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open") | |
861 | <<"Inputfile "<<fname<<" does not exist"<<ENDLOG; | |
862 | return 0; | |
863 | } | |
864 | fTracks = new AliL3TrackArray(); | |
865 | tfile->Binary2TrackArray(fTracks); | |
866 | tfile->CloseBinaryInput(); | |
867 | delete tfile; | |
78127e35 | 868 | printf("Looking in slice %d\n",f); |
1f79afc0 | 869 | for(Int_t i=0; i<fTracks->GetNTracks(); i++) |
870 | { | |
871 | ||
872 | AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i); | |
873 | if(!track) continue; | |
874 | ||
875 | track->CalculateHelix(); | |
876 | UInt_t *hitnum = track->GetHitNumbers(); | |
877 | UInt_t id; | |
878 | ||
879 | Float_t xyz[3]; | |
880 | Int_t padrow; | |
78127e35 | 881 | for(Int_t j=0; j<track->GetNumberOfPoints()-1; j++) |
1f79afc0 | 882 | { |
883 | id = hitnum[j]; | |
884 | Int_t slice = (id>>25) & 0x7f; | |
885 | Int_t patch = (id>>22) & 0x7; | |
886 | UInt_t pos = id&0x3fffff; | |
887 | ||
78127e35 | 888 | //if(slice!=1) continue; |
889 | ||
1f79afc0 | 890 | AliL3SpacePointData *points = fClusters[slice][patch]; |
891 | ||
892 | if(!points) | |
893 | { | |
894 | LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray") | |
895 | <<"No points at slice "<<slice<<" patch "<<patch<<" pos "<<pos<<ENDLOG; | |
896 | continue; | |
897 | } | |
898 | if(pos>=fNcl[slice][patch]) | |
899 | { | |
900 | LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray") | |
901 | <<AliL3Log::kDec<<"ERROR"<<ENDLOG; | |
902 | continue; | |
903 | } | |
904 | ||
905 | xyz[0] = points[pos].fX; | |
906 | xyz[1] = points[pos].fY; | |
907 | xyz[2] = points[pos].fZ; | |
908 | padrow = points[pos].fPadRow; | |
494fad94 | 909 | //AliL3Transform::Global2Local(xyz,slice); |
1f79afc0 | 910 | |
911 | Float_t xyz_cross[3]; | |
912 | track->GetCrossingPoint(padrow,xyz_cross); | |
e1a32fa5 | 913 | Double_t beta = track->GetCrossingAngle(padrow); |
1f79afc0 | 914 | |
915 | Double_t yres = xyz_cross[1] - xyz[1]; | |
916 | Double_t zres = xyz_cross[2] - xyz[2]; | |
917 | ||
918 | Double_t dipangle = atan(track->GetTgl()); | |
919 | ntuppel->Fill(yres,zres,xyz_cross[2],track->GetPt(),dipangle,beta,padrow,track->GetNumberOfPoints()); | |
920 | ||
921 | } | |
922 | } | |
923 | if(fTracks) | |
924 | delete fTracks; | |
925 | } | |
926 | return ntuppel; | |
927 | } | |
928 | ||
e1a32fa5 | 929 | TNtuple *AliL3Evaluate::EvaluatePoints(Char_t *rootfile) |
78127e35 | 930 | { |
e1a32fa5 | 931 | //Compare points to the exact crossing points of track and padrows. |
1b76fd18 | 932 | //The input file to this function, contains the exact clusters calculated |
933 | //in AliTPC::Hits2ExactClusters. | |
e1a32fa5 | 934 | |
355debe1 | 935 | cout<<"Evaluating points"<<endl; |
e1a32fa5 | 936 | TNtuple *ntuppel = new TNtuple("ntuppel","residuals","slice:padrow:resy:resz:zHit:pt"); |
78127e35 | 937 | |
355debe1 | 938 | TFile *exfile = TFile::Open(rootfile); |
939 | if(!exfile) | |
940 | { | |
941 | cerr<<"Error opening rootfile "<<rootfile<<endl; | |
942 | return 0; | |
943 | } | |
944 | ||
e1a32fa5 | 945 | AliTPCParam *param = (AliTPCParam*)exfile->Get("75x40_100x60"); |
946 | ||
947 | //Get the exact clusters from file: | |
948 | AliTPCClustersArray *arr = new AliTPCClustersArray; | |
949 | arr->Setup(param); | |
950 | arr->SetClusterType("AliComplexCluster"); | |
951 | char treeName[500]; | |
952 | sprintf(treeName,"TreeCExact_%s",param->GetTitle()); | |
294ee877 | 953 | Bool_t clusterok = arr->ConnectTree(treeName);//Segment Tree (for offline clusters) |
e1a32fa5 | 954 | if(!clusterok) {printf("AliL3Evaluate::EvaluatePoints : Error in clusterloading\n"); return 0;} |
955 | ||
355debe1 | 956 | cout<<"Entering loop with "<<(Int_t)arr->GetTree()->GetEntries()<<endl; |
e1a32fa5 | 957 | for(Int_t i=0; i<arr->GetTree()->GetEntries(); i++) |
78127e35 | 958 | { |
e1a32fa5 | 959 | //Get the exact clusters for this row: |
960 | Int_t cursec,currow; | |
961 | AliSegmentID *s = arr->LoadEntry(i); | |
962 | param->AdjustSectorRow(s->GetID(),cursec,currow); | |
78127e35 | 963 | |
e1a32fa5 | 964 | AliTPCClustersRow *ro = (AliTPCClustersRow *)arr->GetRow(cursec,currow); |
965 | TClonesArray *clusters = ro->GetArray(); | |
966 | int num_of_offline=clusters->GetEntriesFast(); | |
78127e35 | 967 | |
e1a32fa5 | 968 | //Get the found clusters: |
969 | Int_t slice,padrow; | |
494fad94 | 970 | AliL3Transform::Sector2Slice(slice,padrow,cursec,currow); |
e1a32fa5 | 971 | if(slice<fMinSlice || slice>fMaxSlice) continue; |
972 | AliL3SpacePointData *points = fClusters[slice][0]; | |
1b76fd18 | 973 | |
974 | Int_t index = fRowid[slice][padrow]; | |
975 | if(!fDigitsTree->GetEvent(index)) | |
976 | printf("AliL3Evaluate::EvaluatePoints : ERROR IN DIGITSTREE\n"); | |
294ee877 | 977 | printf("Checking slice %d padrow %d with %d clusters\n",slice,padrow,num_of_offline); |
1b76fd18 | 978 | |
e1a32fa5 | 979 | for(UInt_t c=0; c<fNcl[slice][0]; c++) |
78127e35 | 980 | { |
e1a32fa5 | 981 | if(points[c].fPadRow!=padrow) continue; |
e1a32fa5 | 982 | for(Int_t m=0; m<num_of_offline; m++) |
78127e35 | 983 | { |
e1a32fa5 | 984 | AliComplexCluster *cluster = (AliComplexCluster *)clusters->UncheckedAt(m); |
985 | Int_t mcId = cluster->fTracks[0]; | |
294ee877 | 986 | //Int_t mcId = cluster->GetLabel(0); |
987 | if(mcId <0) continue; | |
e1a32fa5 | 988 | TParticle *part = gAlice->Particle(mcId); |
294ee877 | 989 | |
e1a32fa5 | 990 | Float_t xyz_cl[3] = {points[c].fX,points[c].fY,points[c].fZ}; |
294ee877 | 991 | Float_t xyz_ex[3]; |
494fad94 | 992 | AliL3Transform::Global2Raw(xyz_cl,cursec,currow); |
e170e0f9 | 993 | if(fDigits->GetTrackID((Int_t)rint(xyz_cl[2]),(Int_t)rint(xyz_cl[1]),0)!=mcId && |
994 | fDigits->GetTrackID((Int_t)rint(xyz_cl[2]),(Int_t)rint(xyz_cl[1]),1)!=mcId && | |
995 | fDigits->GetTrackID((Int_t)rint(xyz_cl[2]),(Int_t)rint(xyz_cl[1]),2)!=mcId) | |
e1a32fa5 | 996 | continue; |
494fad94 | 997 | AliL3Transform::Raw2Local(xyz_ex,cursec,currow,cluster->fY,cluster->fX); |
998 | AliL3Transform::Raw2Local(xyz_cl,cursec,currow,xyz_cl[1],xyz_cl[2]); | |
294ee877 | 999 | Float_t resy = xyz_cl[1] - xyz_ex[1];//cluster->GetY() |
1000 | Float_t resz = xyz_cl[2] - xyz_ex[2];//cluster->GetZ() | |
e1a32fa5 | 1001 | |
294ee877 | 1002 | ntuppel->Fill(slice,padrow,resy,resz,xyz_ex[2],part->Pt()); |
78127e35 | 1003 | } |
e1a32fa5 | 1004 | } |
294ee877 | 1005 | arr->ClearRow(cursec,currow); |
78127e35 | 1006 | } |
e1a32fa5 | 1007 | |
78127e35 | 1008 | return ntuppel; |
1009 | } | |
355debe1 | 1010 | |
1011 | void AliL3Evaluate::GetCFeff(Char_t *outfile) | |
78127e35 | 1012 | { |
78127e35 | 1013 | |
355debe1 | 1014 | TNtuple *ntuppel = new TNtuple("ntuppel","Cluster finder efficiency","row:ncrossings:nclusters"); |
1015 | ||
1016 | AliTPC *TPC = (AliTPC*)gAlice->GetDetector("TPC"); | |
1017 | ||
1018 | TPC->SetParam(fParam); | |
78127e35 | 1019 | |
355debe1 | 1020 | Int_t ver = TPC->IsVersion(); |
494fad94 | 1021 | LOG(AliL3Log::kInformational,"AliL3Evaluate::GetCFeff","TPC version") |
355debe1 | 1022 | <<"TPC version "<<ver<<" found on file"<<ENDLOG; |
1023 | ||
1024 | Int_t zero=TPC->GetParam()->GetZeroSup(); | |
1025 | ||
1026 | Int_t np = gAlice->GetNtrack(); | |
1027 | ||
1028 | Int_t crossed,recs; | |
1029 | Int_t *count = new Int_t[np]; //np number of particles. | |
1030 | Int_t i; | |
1031 | Float_t xyz[3]; | |
1032 | for (i=0; i<np; i++) count[i]=0; | |
1033 | for(Int_t sl=fMinSlice; sl<=fMaxSlice; sl++) | |
78127e35 | 1034 | { |
355debe1 | 1035 | for (i=0; i<=175; i++) |
78127e35 | 1036 | { |
355debe1 | 1037 | crossed=0; |
1038 | recs=0; | |
1039 | Int_t index = fRowid[sl][i]; | |
1040 | if (!fDigitsTree->GetEvent(index)) continue; | |
1041 | Int_t sec,row; | |
1042 | fParam->AdjustSectorRow(fDigits->GetID(),sec,row); | |
1043 | fDigits->First(); | |
1044 | do { | |
1045 | Int_t it=fDigits->CurrentRow(), ip=fDigits->CurrentColumn(); | |
1046 | Short_t dig = fDigits->GetDigit(it,ip); | |
1047 | ||
1048 | if(dig<=fParam->GetZeroSup()) continue; | |
1049 | if(it < fParam->GetMaxTBin()-1 && it > 0) | |
1050 | if(fDigits->GetDigit(it+1,ip) <= fParam->GetZeroSup() | |
1051 | && fDigits->GetDigit(it-1,ip) <= fParam->GetZeroSup()) | |
1052 | continue; | |
1053 | ||
494fad94 | 1054 | AliL3Transform::Raw2Local(xyz,sec,row,ip,it); |
355debe1 | 1055 | if(fParam->GetPadRowRadii(sec,row)<230./250.*fabs(xyz[2])) |
1056 | continue; | |
1057 | ||
1058 | ||
1059 | Int_t idx0=fDigits->GetTrackID(it,ip,0); | |
1060 | Int_t idx1=fDigits->GetTrackID(it,ip,1); | |
1061 | Int_t idx2=fDigits->GetTrackID(it,ip,2); | |
1062 | ||
1063 | if (idx0>=0 && dig>=zero) count[idx0]+=1; | |
1064 | if (idx1>=0 && dig>=zero) count[idx1]+=1; | |
1065 | if (idx2>=0 && dig>=zero) count[idx2]+=1; | |
1066 | } while (fDigits->Next()); | |
1067 | ||
1068 | for (Int_t j=0; j<np; j++) | |
78127e35 | 1069 | { |
355debe1 | 1070 | if (count[j]>1) //at least two digits at this padrow |
494fad94 | 1071 | { |
1072 | crossed++; | |
1073 | count[j]=0; | |
1074 | } | |
355debe1 | 1075 | } |
1076 | AliL3SpacePointData *points = fClusters[sl][0]; | |
494fad94 | 1077 | for(UInt_t k=0; k<fNcl[sl][0]; k++) |
355debe1 | 1078 | { |
494fad94 | 1079 | if(points[k].fPadRow!=i) continue; |
355debe1 | 1080 | recs++; |
78127e35 | 1081 | } |
355debe1 | 1082 | ntuppel->Fill(i,crossed,recs); |
78127e35 | 1083 | } |
494fad94 | 1084 | |
78127e35 | 1085 | } |
355debe1 | 1086 | delete[] count; |
1087 | ||
1088 | TFile *file = TFile::Open(outfile,"RECREATE"); | |
1089 | ntuppel->Write(); | |
1090 | file->Close(); | |
1091 | ||
78127e35 | 1092 | } |
355debe1 | 1093 | |
78127e35 | 1094 | Bool_t AliL3Evaluate::GetParticleCrossingPoint(TParticle *part,Int_t slice,Int_t padrow,Float_t *xyz) |
1095 | { | |
1096 | //Calcluate the crossing point between a generated particle and given padrow. | |
1097 | ||
1098 | ||
1099 | Double_t kappa = 0.2*0.0029980/part->Pt(); | |
1100 | ||
1101 | Double_t radius = 1/fabs(kappa); | |
1102 | if(part->GetPdgCode() > 0) kappa = -kappa; | |
1103 | ||
23908b9b | 1104 | Int_t charg; |
1105 | if(kappa>0) | |
1106 | charg=-1; | |
1107 | else | |
1108 | charg=1; | |
1109 | ||
78127e35 | 1110 | Float_t angl[1] = {part->Phi()}; |
1111 | ||
494fad94 | 1112 | AliL3Transform::Global2LocalAngle(angl,slice); |
78127e35 | 1113 | |
1114 | Double_t charge = -1.*kappa; | |
84122e15 | 1115 | Double_t trackPhi0 = angl[0] + charge*0.5*AliL3Transform::Pi()/fabs(charge); |
78127e35 | 1116 | |
1117 | Double_t x0=0; | |
1118 | Double_t y0=0; | |
1119 | Double_t xc = x0 - radius * cos(trackPhi0); | |
1120 | Double_t yc = y0 - radius * sin(trackPhi0); | |
1121 | ||
1122 | //printf("radius %f xc %f yc %f\n",radius,xc,yc); | |
1123 | ||
494fad94 | 1124 | Double_t xHit = AliL3Transform::Row2X(padrow); |
78127e35 | 1125 | xyz[0] = xHit; |
1126 | Double_t aa = (xHit - xc)*(xHit - xc); | |
1127 | Double_t r2 = radius*radius; | |
1128 | if(aa > r2) | |
1129 | return false; | |
1130 | ||
1131 | Double_t aa2 = sqrt(r2 - aa); | |
1132 | Double_t y1 = yc + aa2; | |
1133 | Double_t y2 = yc - aa2; | |
1134 | xyz[1] = y1; | |
1135 | if(fabs(y2) < fabs(y1)) xyz[1] = y2; | |
1136 | ||
23908b9b | 1137 | Double_t tgl = part->Pz()/part->Pt(); |
1138 | ||
1139 | Double_t yHit = xyz[1]; | |
1140 | Double_t angle1 = atan2((yHit - yc),(xHit - xc)); | |
84122e15 | 1141 | if(angle1 < 0) angle1 += 2.*AliL3Transform::Pi(); |
23908b9b | 1142 | Double_t angle2 = atan2((0 - yc),(0 - xc)); |
84122e15 | 1143 | if(angle2 < 0) angle2 += 2.*AliL3Transform::Pi(); |
23908b9b | 1144 | Double_t diff_angle = angle1 - angle2; |
84122e15 | 1145 | diff_angle = fmod(diff_angle,2*AliL3Transform::Pi()); |
1146 | if((charg*diff_angle) > 0) diff_angle = diff_angle - charg*2.*AliL3Transform::Pi(); | |
23908b9b | 1147 | Double_t s_tot = fabs(diff_angle)*radius; |
1148 | Double_t zHit = 0 + s_tot*tgl; | |
1149 | xyz[2] = zHit; | |
1150 | ||
78127e35 | 1151 | return true; |
1152 | } | |
1153 |