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