03a14309b94e0ebe4305eac9bb79c3a0f205bda6
[u/mrichter/AliRoot.git] / HLT / src / AliL3Evaluate.cxx
1 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
2 //*-- Copyright &copy ASV 
3
4 #include <stdio.h>
5 #include <math.h>
6 #include <fstream.h>
7 #include <TFile.h>
8 #include <TH1.h>
9 #include <TParticle.h>
10 #include <TTree.h>
11 #include <TClonesArray.h>
12
13 #include "AliRun.h"
14 #include "AliSimDigits.h"
15 #include "AliTPC.h"
16 #include "AliTPCClustersArray.h"
17 #include "AliTPCClustersRow.h"
18 #include "AliTPCcluster.h"
19 #include "AliTPCParam.h"
20 #include "AliComplexCluster.h"
21
22 #include "AliL3Defs.h"
23 #include "AliL3Transform.h"
24 #include "AliL3SpacePointData.h"
25 #include "AliL3Track.h"
26 #include "AliL3FileHandler.h"
27 #include "AliL3TrackArray.h"
28 #include "AliL3Evaluate.h"
29 #include "AliL3Logging.h"
30
31 //_____________________________________________________________
32 // AliL3Evaluate
33 //
34 // Evaluation class for tracking. Plots efficiencies etc..
35
36 ClassImp(AliL3Evaluate)
37
38 AliL3Evaluate::AliL3Evaluate()
39 {
40   fDigitsFile = NULL;
41   fTracks = NULL;
42   fMCclusterfile = NULL;
43   fNFastPoints = 0;
44   fMcIndex = 0;
45   fMcId = 0;
46   fMinSlice=0;
47   fMaxSlice=0;
48   fTransform = new AliL3Transform();
49 }
50
51 AliL3Evaluate::AliL3Evaluate(Char_t *mcfile,Int_t *slice)
52 {
53   //Normal constructor. Input are the rootfile containing the 
54   //original MC information of the simulated event. 
55
56   fEventFile = new TFile(mcfile,"READ");
57     
58   fParam = (AliTPCParam*)fEventFile->Get("75x40_100x60");
59   fTransform = new AliL3Transform();
60   
61   fMinSlice = slice[0];
62   fMaxSlice = slice[1];
63   fMinGoodPt = 0.1;
64   fNoOverlap = kFALSE;
65 }
66
67 AliL3Evaluate::AliL3Evaluate(Int_t *slice)
68 {
69   //ctor to use if you do not need any rootfile.
70
71
72   fMinSlice = slice[0];
73   fMaxSlice = slice[1];
74   fTransform = new AliL3Transform();
75
76 }
77
78 AliL3Evaluate::~AliL3Evaluate()
79 {
80   if(fDigitsTree) fDigitsTree->Delete();
81   if(fDigitsFile) {
82     fDigitsFile->Close();
83     delete fDigitsFile;
84   }
85   if(fTransform) delete fTransform;
86   if(fTracks) delete fTracks;
87   if(fPtRes) delete fPtRes;
88   if(fNGoodTracksPt) delete fNGoodTracksPt;
89   if(fNFoundTracksPt) delete fNFoundTracksPt;
90   if(fNFakeTracksPt) delete fNFakeTracksPt;
91   if(fTrackEffPt) delete fTrackEffPt;
92   if(fFakeTrackEffPt) delete fFakeTrackEffPt;
93   if(fNGoodTracksEta) delete fNGoodTracksEta;
94   if(fNFoundTracksEta) delete fNFoundTracksEta;
95   if(fNFakeTracksEta) delete fNFakeTracksEta;
96   if(fTrackEffEta) delete fTrackEffEta;
97   if(fFakeTrackEffEta) delete fFakeTrackEffEta;
98   if(fMcIndex) delete [] fMcIndex;
99   if(fMcId)    delete [] fMcId;
100   if(fNtuppel) delete fNtuppel;
101 }
102
103 void AliL3Evaluate::Setup(Char_t *trackfile,Char_t *path)
104 {
105   //Read in the hit and track information from produced files.
106   
107   fGoodFound = 0;
108   fGoodGen = 0;
109   Char_t fname[256];
110   AliL3FileHandler *clusterfile[36][6];
111   for(Int_t s=fMinSlice; s<=fMaxSlice; s++)
112     {
113       for(Int_t p=0; p<6; p++)
114         {
115           fClusters[s][p] = 0;
116           clusterfile[s][p] = new AliL3FileHandler();
117           sprintf(fname,"%s/points_%d_%d.raw",path,s,p);
118           if(!clusterfile[s][p]->SetBinaryInput(fname))
119             {
120               LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
121                 <<"Inputfile "<<fname<<" does not exist"<<ENDLOG; 
122               delete clusterfile[s][p];
123               clusterfile[s][p] = 0; 
124               continue;
125             }
126           fClusters[s][p] = (AliL3SpacePointData*)clusterfile[s][p]->Allocate();
127           clusterfile[s][p]->Binary2Memory(fNcl[s][p],fClusters[s][p]);
128           clusterfile[s][p]->CloseBinaryInput();
129         }
130     }
131
132   
133   AliL3FileHandler *tfile = new AliL3FileHandler();
134   if(!tfile->SetBinaryInput(trackfile)){
135     LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
136     <<"Inputfile "<<trackfile<<" does not exist"<<ENDLOG; 
137     return;
138   }
139   fTracks = new AliL3TrackArray();
140   tfile->Binary2TrackArray(fTracks);
141   tfile->CloseBinaryInput();
142   delete tfile;
143 }
144
145 void AliL3Evaluate::SetupSlow(Char_t *trackfile,Char_t *digitsfile,Char_t *path)
146 {
147   //Setup for using the slow simulator.
148   
149   fDigitsFile = new TFile(digitsfile,"READ");
150   if(!fDigitsFile->IsOpen())
151     {
152       LOG(AliL3Log::kError,"AliL3Evaluate::SetupSlow","File Open")
153         <<"Inputfile "<<digitsfile<<" does not exist"<<ENDLOG;
154       return;
155     }
156   fIsSlow = true;
157   Setup(trackfile,path);
158   
159   if(!InitMC())
160     LOG(AliL3Log::kError,"AliL3Evaluation::SetupSlow","Digits Tree")
161     <<"Error setting up digits tree"<<ENDLOG;
162   
163 }
164
165 void AliL3Evaluate::SetupFast(Char_t *trackfile,Char_t *mcClusterfile,Char_t *path)
166 {
167   //Setup for using the fast simulator.
168
169   fIsSlow = false;
170   GetFastClusterIDs(path);
171   
172   fMCclusterfile = new TFile(mcClusterfile);
173   if(!fMCclusterfile->IsOpen())
174     LOG(AliL3Log::kError,"AliL3Evaluation::SetupFast","File Open")
175       <<"Inputfile "<<mcClusterfile<<" does not exist"<<ENDLOG; 
176
177   Setup(trackfile,path);
178   InitMC();
179 }
180
181 Bool_t AliL3Evaluate::InitMC()
182 {
183   if(fIsSlow)
184     {
185       fDigitsFile->cd();
186       fDigitsTree = (TTree*)fDigitsFile->Get("TreeD_75x40_100x60_0");
187       if(!fDigitsTree) 
188         {
189           LOG(AliL3Log::kError,"AliL3Evaluate::InitMC","Digits Tree")
190             <<AliL3Log::kHex<<"Error getting digitstree "<<(Int_t)fDigitsTree<<ENDLOG;
191           return false;
192         }
193       fDigitsTree->GetBranch("Segment")->SetAddress(&fDigits);
194       for(Int_t i=0; i<fDigitsTree->GetEntries(); i++)
195         {
196           if(!fDigitsTree->GetEvent(i)) continue;
197           Int_t se,ro,slice,slicerow;
198           fParam->AdjustSectorRow(fDigits->GetID(),se,ro);
199           fTransform->Sector2Slice(slice,slicerow,se,ro);
200           fRowid[slice][slicerow] = i;
201         }
202     }
203   
204   fEventFile->cd();
205   gAlice = (AliRun*)fEventFile->Get("gAlice");
206   if (!gAlice) 
207     {
208       LOG(AliL3Log::kError,"AliL3Evaluate::InitMC","gAlice")
209         <<"AliRun object non existing on file"<<ENDLOG;
210       return false;
211     }
212   
213   gAlice->GetEvent(0);
214   
215   return true;
216   
217 }
218
219
220
221 void AliL3Evaluate::DefineGoodTracks(Int_t *slice,Int_t *padrow,Int_t good_number,Char_t *fname)
222 {
223   //Loop over MC particles, and mark the good ones
224   //(which the tracker should find...)
225   
226   
227   ifstream in(fname);
228   if(in)
229     {
230       LOG(AliL3Log::kInformational,"AliL3Evaluate::DefineGoodTracks","infile")
231         <<"Reading good particles from file "<<fname<<ENDLOG;
232       while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>>
233              fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>>
234              fGoodTracks[fGoodGen].pt) 
235         {
236           fGoodGen++;
237           if (fGoodGen==15000) 
238             {
239               LOG(AliL3Log::kWarning,"AliL3Evaluate::DefineGoodTracks","fGoodGen")
240                 <<"Too many good tracks"<<ENDLOG;
241               break;
242             }
243         }
244       if (!in.eof()) cerr<<"Read error (good_tracks_tpc) !\n";
245       return;
246     }
247   
248   
249   AliTPC *TPC = (AliTPC*)gAlice->GetDetector("TPC");
250   
251   TPC->SetParam(fParam);
252   
253   Int_t ver = TPC->IsVersion();
254   LOG(AliL3Log::kInformational,"AliL3Evaluate::DefineGoodTracks","TPC version")
255     <<"TPC version "<<ver<<" found on file"<<ENDLOG;
256   
257   Int_t zero=TPC->GetParam()->GetZeroSup();
258   
259   Int_t np = gAlice->GetNtrack();
260   Int_t *good = new Int_t[np];
261   for(Int_t ii=0; ii<np; ii++)
262     good[ii] = 0;
263
264   if(ver==1)
265     {
266       if(fIsSlow)
267         LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","TPC version")
268           <<"TPC version "<<ver<<" does not match."<<ENDLOG;
269       fMCclusterfile->cd();
270       AliTPCClustersArray carray;
271       carray.Setup(fParam);
272       carray.SetClusterType("AliTPCcluster");
273       Char_t cname[100];
274       Int_t eventn = 0;
275       sprintf(cname,"TreeC_TPC_%d",eventn);
276       Bool_t clusterok = carray.ConnectTree(cname);
277       if(!clusterok) 
278         LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","Cluster Array")
279           <<"Error loading clusters from rootfile"<<ENDLOG;
280          
281       for(Int_t i=0; i<carray.GetTree()->GetEntries(); i++)  
282         {
283           Int_t sec,row,sl,lr;
284           AliSegmentID *s = carray.LoadEntry(i);
285           fParam->AdjustSectorRow(s->GetID(),sec,row);
286           fTransform->Sector2Slice(sl,lr,sec,row);
287           
288           if(sl != slice[0]) {carray.ClearRow(sec,row); continue;}
289           if(lr < padrow[0]) {carray.ClearRow(sec,row); continue;}
290           if(lr > padrow[1]) {carray.ClearRow(sec,row); continue;}
291           AliTPCClustersRow *cRow = carray.GetRow(sec,row);
292           for(Int_t j=0; j<cRow->GetArray()->GetEntriesFast(); j++)
293             {
294               AliTPCcluster *cluster=(AliTPCcluster*)(*cRow)[j];
295               Int_t lab=cluster->GetLabel(0);
296               if(lab<0) continue;
297               lab=TMath::Abs(lab);
298               good[lab]++;
299             }
300           if(carray.GetRow(sec,row)) 
301             carray.ClearRow(sec,row);
302         }
303     }
304   else if(ver==2)
305     {
306       if(!fIsSlow)
307         LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","TPC version")
308           <<"TPC version "<<ver<<" does not match."<<ENDLOG;
309       Int_t *count = new Int_t[np]; //np number of particles.
310       Int_t i;
311       Float_t xyz[3];
312       for (i=0; i<np; i++) count[i]=0;
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);
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               
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);
338
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               }
351           }
352         }
353       delete[] count;
354     }
355   
356   else 
357     {
358       LOG(AliL3Log::kError,"AliL3Evaluation::FillEffHistos","TPC version")
359         <<"No valid TPC version found"<<ENDLOG;
360       return;
361     }
362   
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();
368   Int_t max = 15000;
369   while (npart--) {
370     AliTPChit *hit0=0;
371     
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
392     if (good[i] < good_number || good[i]>176) continue;
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   }
407   
408   delete [] good;
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   
428 }
429
430 void 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   
434   Int_t row[6][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,175}};
435   Int_t sl[2] ={slice,slice};
436   DefineGoodTracks(sl,row[patch],good_number);
437   SetMinPoints(min_points);
438   AssignIDs();
439   CreateHistos();
440   FillEffHistos();
441   CalcEffHistos();
442 }
443
444 void 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
450   Int_t row[2] = {0,175};
451   Int_t sl[2] ={slice,slice};
452   DefineGoodTracks(sl,row,good_number);
453
454   SetMinPoints(min_points);
455   
456   AssignIDs();
457   CreateHistos();
458   FillEffHistos();
459   CalcEffHistos();
460 }
461
462 void AliL3Evaluate::EvaluateGlobal(Int_t min_points,Int_t good_number,Char_t *fname)
463 {
464   //Make efficiency plots for tracking on several slices.
465   
466   Int_t row[2] = {0,175};
467   SetMinPoints(min_points);
468   Int_t slice[2] = {fMinSlice,fMaxSlice};
469   DefineGoodTracks(slice,row,good_number,fname);
470   AssignIDs();
471   CreateHistos(5,0,2);
472   printf("filling histos\n");
473   FillEffHistos();
474   printf("done fillling\n");
475   CalcEffHistos();
476 }
477
478 void AliL3Evaluate::AssignIDs()
479 {
480   //Assign MC id to the tracks.
481   
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);
488       if(!track) continue; 
489       if(track->GetNumberOfPoints() < fMinPointsOnTrack) break;
490       
491       fGoodFound++;
492       Int_t tID = GetMCTrackLabel(track);
493       track->SetMCid(tID);
494       printf("track %i id %d nHits %d\n",i,tID,track->GetNumberOfPoints());
495     }
496 }
497
498
499 struct S {Int_t lab; Int_t max;};
500 Int_t AliL3Evaluate::GetMCTrackLabel(AliL3Track *track){ 
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.
504
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;
511     
512   //Int_t **trackID = GetClusterIDs(track);
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];
524       if(!points) continue; 
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:
533       
534       lab=abs(points[pos].fTrackID[0]);
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];
556       if(!points) continue; 
557       if(pos>=fNcl[slice][patch]) 
558         {
559           LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
560             <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
561           continue;
562         }
563       
564       if (abs(points[pos].fTrackID[1]) == lab || 
565           abs(points[pos].fTrackID[2]) == lab ) max++;
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   
574   return lab;
575 }
576
577
578 Int_t **AliL3Evaluate::GetClusterIDs(AliL3Track *track)
579 {
580   //Return the MC information of all clusters belonging to track.
581
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;
587   
588   Float_t xyz[3];
589   Int_t padrow;
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];
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];
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             }*/
640         }
641       else
642         {
643           Int_t tmp_pid=0;
644           for(Int_t ii=0; ii<fNFastPoints; ii++)
645             {
646               tmp_pid = fMcId[ii];
647               if(fMcIndex[ii] == id) break;
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     }
655
656   return trackID;
657 }
658
659 void AliL3Evaluate::GetFastClusterIDs(Char_t *path)
660 {
661   //Get the MC id of space points in case of using the fast simulator. 
662   char fname[256];
663   sprintf(fname,"%s/point_mc.dat",path);
664   FILE *infile = fopen(fname,"r");
665   if(!infile) return;
666   Int_t hitid,hitmc,i;
667   
668   for(i=0; ; i++)
669     if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
670   rewind(infile);
671   fNFastPoints = i;
672   fMcId = new Int_t[fNFastPoints];
673   fMcIndex = new UInt_t[fNFastPoints];
674   
675   for(i=0; i<fNFastPoints; i++)
676     {
677       if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
678       fMcId[i] = hitmc;
679       fMcIndex[i] = hitid;
680     }
681   fclose(infile);
682 }
683
684 void AliL3Evaluate::CreateHistos(Int_t nbin,Int_t xlow,Int_t xup)
685 {
686   //Create the histograms 
687   
688   LOG(AliL3Log::kInformational,"AliL3Evaluate::CreateHistos","Allocating")
689     <<"Creating histograms..."<<ENDLOG;
690   
691   fNtuppel = new TNtuple("fNtuppel","Pt resolution","pt_gen:pt_found:nHits");
692
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);
705   printf("finished creating histos\n");
706 }
707 /*
708 void AliL3Evaluate::FillEffHistos()
709 {  
710   //Fill the efficiency histograms.
711   
712   for(Int_t i=0; i<fGoodGen; i++)
713     {
714       Double_t ptg=fGoodTracks[i].pt,pzg=fGoodTracks[i].pz;
715       Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
716       //printf("filling particle with pt %f and dipangle %f\n",ptg,dipangle);
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();
729           //printf("evaluating track id %d\n",tracklabel);
730           if(TMath::Abs(tracklabel) != fGoodTracks[i].label) continue;
731           found=1;
732           if(tracklabel == fGoodTracks[i].label) {fNFoundTracksPt->Fill(ptg); fNFoundTracksEta->Fill(dipangle);}
733           else {fNFakeTracksPt->Fill(ptg); fNFakeTracksEta->Fill(dipangle);}
734           Float_t pt=track->GetPt();
735           fPtRes->Fill((pt-ptg)/ptg*100.);
736           fNtuppel->Fill(ptg,pt,nHits);
737           break;
738           
739         }
740     }
741 }
742 */
743 void 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 }
774
775 void 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   
790   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
791     <<AliL3Log::kDec<<"Naive efficiency "<<(Double_t)fGoodFound/(Double_t)fGoodGen<<ENDLOG;
792
793   
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        
814 }
815
816 void 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();
829   fNtuppel->Write();
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
847 TNtuple *AliL3Evaluate::CalculateResiduals()
848 {
849
850   TNtuple *ntuppel=new TNtuple("ntuppel","Residuals","residual_trans:residual_long:zHit:pt:dipangle:beta:padrow:nHits");
851
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;
866       printf("Looking in slice %d\n",f);
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;
879           for(Int_t j=0; j<track->GetNumberOfPoints()-1; j++)
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               
886               //if(slice!=1) continue;
887               
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;
907               //fTransform->Global2Local(xyz,slice);
908               
909               Float_t xyz_cross[3];
910               track->GetCrossingPoint(padrow,xyz_cross);
911               Double_t beta = track->GetCrossingAngle(padrow);
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
927 TNtuple *AliL3Evaluate::EvaluatePoints(Char_t *rootfile)
928 {
929   //Compare points to the exact crossing points of track and padrows.
930   //The input file to this function, contains the exact clusters calculated
931   //in AliTPC::Hits2ExactClusters.
932     
933   cout<<"Evaluating points"<<endl;
934   TNtuple *ntuppel = new TNtuple("ntuppel","residuals","slice:padrow:resy:resz:zHit:pt");
935   
936   TFile *exfile = TFile::Open(rootfile);
937   if(!exfile)
938     {
939       cerr<<"Error opening rootfile "<<rootfile<<endl;
940       return 0;
941     }
942
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());
951   Bool_t clusterok = arr->ConnectTree(treeName);//Segment Tree (for offline clusters)
952   if(!clusterok) {printf("AliL3Evaluate::EvaluatePoints : Error in clusterloading\n"); return 0;}
953   
954   cout<<"Entering loop with "<<(Int_t)arr->GetTree()->GetEntries()<<endl;
955   for(Int_t i=0; i<arr->GetTree()->GetEntries(); i++)
956     {
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);
961       
962       AliTPCClustersRow *ro = (AliTPCClustersRow *)arr->GetRow(cursec,currow);
963       TClonesArray *clusters = ro->GetArray();
964       int num_of_offline=clusters->GetEntriesFast();
965       
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];
971       
972       Int_t index = fRowid[slice][padrow];
973       if(!fDigitsTree->GetEvent(index))
974         printf("AliL3Evaluate::EvaluatePoints : ERROR IN DIGITSTREE\n");
975       printf("Checking slice %d padrow %d with %d clusters\n",slice,padrow,num_of_offline);
976       
977       for(UInt_t c=0; c<fNcl[slice][0]; c++)
978         {
979           if(points[c].fPadRow!=padrow) continue;
980           for(Int_t m=0; m<num_of_offline; m++)
981             {
982               AliComplexCluster *cluster = (AliComplexCluster *)clusters->UncheckedAt(m);
983               Int_t mcId = cluster->fTracks[0];
984               //Int_t mcId = cluster->GetLabel(0);
985               if(mcId <0) continue;
986               TParticle *part = gAlice->Particle(mcId);
987               
988               Float_t xyz_cl[3] = {points[c].fX,points[c].fY,points[c].fZ};
989               Float_t xyz_ex[3];
990               fTransform->Global2Raw(xyz_cl,cursec,currow);
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)
994                 continue;
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()
999               
1000               ntuppel->Fill(slice,padrow,resy,resz,xyz_ex[2],part->Pt());
1001             }
1002         }      
1003       arr->ClearRow(cursec,currow);
1004     }
1005
1006   return ntuppel;
1007 }
1008
1009 void AliL3Evaluate::GetCFeff(Char_t *outfile)
1010 {
1011   
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);
1017   
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++)
1032     {
1033       for (i=0; i<=175; i++) 
1034         {
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++) 
1067             {
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++;
1077             }
1078           ntuppel->Fill(i,crossed,recs);
1079         }
1080     }
1081   delete[] count;
1082   
1083   TFile *file = TFile::Open(outfile,"RECREATE");
1084   ntuppel->Write();
1085   file->Close();
1086   
1087 }
1088
1089 Bool_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   
1099   Int_t charg;
1100   if(kappa>0)
1101     charg=-1;
1102   else
1103     charg=1;
1104   
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
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
1146   return true;
1147 }
1148