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