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