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