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