Use new to create arrays with variable size (Sun)
[u/mrichter/AliRoot.git] / TOF / AliTOFReconstructioner.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //_________________________________________________________________________
17 // Manager class for TOF reconstruction.
18 // 
19 //
20 //-- Authors: Bologna-ITEP-Salerno Group
21 //
22 // Description: Manager class for TOF reconstruction (derived from TTask)
23 // Summary of the main methods:
24 // - extraction of the TPC (assumed to be) reconstructed tracks 
25 //   comment: it has to me moved as soon as possible into a separate
26 //   class AliTOFTrackReader (K. Safarik suggestion)
27 // - geometrical propagation of the above tracks till TOF detector
28 // - matching of the tracks with the TOF signals
29 // 
30 // Remark: the GEANT3.21 geometry is used during the geometrical propagation
31 // of the tracks in order to know the current volume reached by the track.
32 //
33 //////////////////////////////////////////////////////////////////////////////
34
35
36 #include "TTask.h"
37 #include "TBenchmark.h"
38 #include "TTree.h"
39 #include "TSystem.h"
40 #include "TFile.h"
41 #include "TParticle.h"
42
43 #include "AliConst.h"
44 #include "AliRun.h"
45 #include "AliTOFConstants.h"
46 #include "AliTOFHitMap.h"
47 #include "AliTOFSDigit.h"
48 #include "AliTOFhit.h"
49 #include "AliTOFRecHit.h"
50 #include "AliTOFPad.h"
51 #include "AliTOFTrack.h"
52 #include "AliTOF.h"
53 #include "AliTOFv1.h"
54 #include "AliTOFv2.h"
55 #include "AliTOFv2FHoles.h"
56 #include "AliTOFv3.h"
57 #include "AliTOFv4.h"
58 #include "AliTOFv4T0.h"
59 #include "AliTOFReconstructioner.h"
60 // this line has to be commented till TPC will provide fPx fPy fPz and fL in
61 // AliTPChit class or somewhere
62 // #include "../TPC/AliTPC.h"
63 #include "AliRun.h"
64 #include "AliDetector.h"
65 #include "AliMC.h"
66
67 #include <TClonesArray.h>
68 #include "../TGeant3/TGeant3.h"
69 #include "TFile.h"
70 #include <TF1.h>
71 #include <TF2.h>
72 #include "TTask.h"
73 #include "TTree.h"
74 #include "TSystem.h"
75 #include "TROOT.h"
76 #include "TFolder.h"
77 #include "TNtuple.h"
78 #include <stdlib.h>
79 #include <iostream.h>
80 #include <fstream.h>
81
82 ClassImp(AliTOFReconstructioner)
83
84 //____________________________________________________________________________ 
85   AliTOFReconstructioner::AliTOFReconstructioner():TTask("AliTOFReconstructioner","") 
86 {
87   // default ctor
88   fNevents = 0 ; 
89   fg3      = 0;
90   foutputfile  = 0; 
91   foutputntuple= 0;
92   fZnoise  = 0;
93   ftail    = 0;
94 }
95            
96 //____________________________________________________________________________ 
97   AliTOFReconstructioner::AliTOFReconstructioner(char* headerFile, Option_t* opt, char *RecFile ):TTask("AliTOFReconstructioner","") 
98 {
99   //
100   // ctor
101   //
102   fNevents = 0 ;     // Number of events to reconstruct, 0 means all evens in current file
103   fg3      = 0;
104   foutputfile  = 0; 
105   foutputntuple= 0;
106   fZnoise  = 0;
107   ftail    = 0;
108
109   Init(opt);
110
111   // create output file
112   if (RecFile){
113     foutputfile= new TFile(RecFile,"RECREATE","root file for matching");
114   } else {
115     char outFileName[100];
116     strcpy(outFileName,"match");
117     strcat(outFileName,headerFile);
118     foutputfile= new TFile(outFileName,"RECREATE","root file for matching");
119   }
120   
121   // initialize the ALIROOT geometry 
122   gAlice->Init();
123   gAlice->Print(); 
124
125   // set the fg3 pointer to geometry used by IsInsideThePad method
126   fg3 = (TGeant3*) gMC;  
127
128   CreateNTuple();  
129
130   // add Task to //root/Tasks folder
131   TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
132   roottasks->Add(this) ; 
133 }
134 //____________________________________________________________________________ 
135 void AliTOFReconstructioner::Init(Option_t* opt)
136 {
137   // Initialize the AliTOFReconstructioner setting parameters for
138   // reconstruction.
139   // Option values: Pb-Pb for Pb-Pb events
140   //                pp    for pp    events
141
142   // set common parameters
143   fdbg=1;
144   fNevents    = 1;
145   fFirstEvent = 1;
146   fLastEvent  = 1;
147   fTimeResolution =0.120;
148   fpadefficiency  =0.99 ;
149   fEdgeEffect     = 2   ;
150   fEdgeTails      = 0   ;
151   fHparameter     = 0.4 ;
152   fH2parameter    = 0.15;
153   fKparameter     = 0.5 ;
154   fK2parameter    = 0.35;
155   fEffCenter      = fpadefficiency;
156   fEffBoundary    = 0.65;
157   fEff2Boundary   = 0.90;
158   fEff3Boundary   = 0.08;
159   fResCenter      = 50. ;
160   fResBoundary    = 70. ;
161   fResSlope       = 40. ;
162   fTimeWalkCenter = 0.  ;
163   fTimeWalkBoundary=0.  ;
164   fTimeWalkSlope  = 0.  ;
165   fTimeDelayFlag  = 1   ;
166   fPulseHeightSlope=2.0 ;
167   fTimeDelaySlope =0.060;
168   // was fMinimumCharge = TMath::Exp(fPulseHeightSlope*fKparameter/2.);
169   fMinimumCharge = TMath::Exp(-fPulseHeightSlope*fHparameter);
170   fChargeSmearing=0.0   ;
171   fLogChargeSmearing=0.13;
172   fTimeSmearing   =0.022;
173   fAverageTimeFlag=0    ;
174   fChargeFactorForMatching=1;
175   fTrackingEfficiency=1.0; // 100% TPC tracking efficiency assumed
176   fSigmavsp = 1.        ;
177   fSigmaZ   = 0.        ;
178   fSigmarphi= 0.        ;
179   fSigmap   = 0.        ;
180   fSigmaPhi = 0.        ;
181   fSigmaTheta=0.        ;
182   fField    = 0.2       ;
183   // fRadLenTPC : 0.2 includes TRD / 0.03 TPC only
184   fRadLenTPC=0.06        ; // last value
185   fCorrectionTRD=0.     ;
186   fLastTPCRow=111       ;
187   fRadiusvtxBound=50.   ; // expressed in [cm]
188   fStep     = 0.1       ; // expressed in [cm] step during propagation of the
189                           // track inside TOF volumes 
190   fMatchingStyle=2      ;
191   /* previous values default
192   fMaxPixels=70000      ;
193   fMaxAllTracks=70000   ;
194   fMaxTracks=15000      ;
195   */
196   fMaxPixels=165000      ;
197   fMaxAllTracks=500000   ;
198   fMaxTracks=15000      ;
199
200   fMaxTOFHits=35000     ;
201   fPBound      =0.0     ; // bending effect: P_t=0.3*z*B*R , z particle charge
202   fNoiseSlope=20.       ;
203   // set parameters as specified in opt
204   //pp case
205   if(strstr(opt,"pp")){
206   fMaxTestTracks=500    ; 
207   fNoise    = 26.       ;
208   fNoiseMeanTof= 26.4       ; // to check
209   }
210   //Pb-Pb case
211   if(strstr(opt,"Pb-Pb")){
212   fMaxTestTracks=20     ;
213   fNoise    = 9400.     ;
214   fNoiseMeanTof= 26.4       ;
215   }
216 }
217
218 //____________________________________________________________________________ 
219   AliTOFReconstructioner::~AliTOFReconstructioner()
220 {
221   //
222   // dtor
223   //
224   if (fg3)
225     {
226       delete fg3;
227       fg3 = 0;
228     }
229   if (foutputfile)
230     {
231       delete foutputfile;
232       foutputfile = 0;
233     }
234   if (foutputntuple)
235     {
236       delete foutputntuple;
237       foutputntuple = 0;
238     }
239
240   if (fZnoise)
241     {
242       delete fZnoise;
243       fZnoise = 0;
244     }
245
246   if (ftail)
247     {
248       delete ftail;
249       ftail = 0;
250     }
251 }
252
253 //____________________________________________________________________________
254 void AliTOFReconstructioner::CreateNTuple()
255 {
256   //
257   // Create a Ntuple where information about reconstructed charged particles 
258   // (both primaries and secondaries) are stored 
259   // Variables: event ipart imam xvtx yvtx zvtx pxvtx pyvtx pzvtx time leng matc text mext
260   // Meaning:
261   // event - event number (0, 1, ...)
262   // ipart - PDG code of particles 
263   // imam  - PDG code for the parent
264   // =0 for primary particle
265   // xvtx  - x-coordinate of the vertex (cm)
266   // yvtx  - y-coordinate of the vertex (cm)
267   // zvtx  - z-coordinate of the vertex (cm)
268   // pxvtx - x-coordinate of the momentum in the vertex (GeV)
269   // pyvtx - y-coordinate of the momentum in the vertex (GeV)
270   // pzvtx - z-coordinate of the momentum in the vertex (GeV)
271   // time  - time of flight from TOF for given track (ps) - TOF time for the
272   //         first TOF hit of the track
273   // leng  - track length to the TOF pixel (cm), evaluate as a sum of the
274   // track length from the track vertex to TPC and the average
275   // length of the extrapolated track from TPC to TOF.
276   // for the track without TOF hits leng=-abs(leng)
277   // matc  - index of the (TPC track) - (TOF pixel) matching
278   // =0 for tracks which are not tracks for matching, i.e. 
279   // there is not hit on the TPC or Rvxt>200 cm
280   // >0 for tracks with positive matching procedure:
281   //   =1 or 2 for non-identified tracks:
282   //     =1, if the corresponding pixel is not fired,
283   //     =2, if the corresponding pixel is also matched to the 
284   //         other track,
285   //   =3 or 4 for identified tracks:
286   //     =3, if identified with true time,
287   //     =4, if identified with wrong time.
288   // <0 for tracks with negative mathing procedure:
289   //   =-1, if track do not reach the pixel plate (curved in the 
290   //        magnetic field),
291   //   =-2, if track is out of z-size of the TOF,
292   //   =-3, if track is or into the RICH hole, or into the PHOS hole, or in the space between the plates,
293   //   =-4, if track is into the dead space of the TOF.
294   // text  - time of fligth from the matching procedure = time of the 
295   //         pixel corresponding to the track (ps)
296   //         =0 for the tracks with matc<=1
297   // mext  - mass of the track from the matching procedure
298   //           =p*sqrt(900*(text/leng)**2-1), if 900*(text/leng)**2-1>=0
299   //           =-p*sqrt(abs(900*(text/leng)**2-1)), if 900*(text/leng)**2-1<0
300
301   foutputntuple= new TNtuple("Ntuple","matching","event:ipart:imam:xvtx:yvtx:zvtx:pxvtx:pyvtx:pzvtx:time:leng:matc:text:mext",2000000); // buffersize set for 25 Pb-Pb events
302 }
303
304 //__________________________________________________________________
305 Double_t TimeWithTailR(Double_t* x, Double_t* par)
306 {
307   // sigma - par[0], alpha - par[1], part - par[2]
308   //  at x<part*sigma - gauss
309   //  at x>part*sigma - TMath::Exp(-x/alpha)
310   Float_t xx =x[0];
311   Double_t f;
312   if(xx<par[0]*par[2]) {
313     f = TMath::Exp(-xx*xx/(2*par[0]*par[0]));
314   } else {
315     f = TMath::Exp(-(xx-par[0]*par[2])/par[1]-0.5*par[2]*par[2]);
316   }
317   return f;
318 }
319
320 //____________________________________________________________________________
321 void AliTOFReconstructioner::Exec(const char* datafile, Option_t *option) 
322
323   //
324   // Performs reconstruction for TOF detector
325   // 
326   gBenchmark->Start("TOFReconstruction");
327
328   TFile *file = TFile::Open(datafile);
329
330   // Get AliRun object from file or create it if not on file
331   gAlice = (AliRun*)file->Get("gAlice");
332
333   AliTOF* TOF = (AliTOF *) gAlice->GetDetector ("TOF");
334   AliDetector* TPC = gAlice->GetDetector("TPC");
335
336   if (!TOF) {
337     Error("AliTOFReconstructioner","TOF not found");
338     return;
339   }
340   if (!TPC) {
341     Error("AliTOFReconstructioner","TPC Detector not found");
342     return;
343   }
344
345   if (fEdgeTails) ftail = new TF1("tail",TimeWithTailR,-2,2,3);
346
347   if (fNevents == 0) fNevents = (Int_t) gAlice->TreeE()->GetEntries();
348   // You have to set the number of event with the ad hoc setter
349   // see testrecon.C
350
351   for (Int_t ievent = 0; ievent < fNevents; ievent++) { // start loop on events
352
353     Int_t nparticles=gAlice->GetEvent(ievent);
354     if (nparticles <= 0) return;
355
356     TClonesArray* tofhits=0;
357     TClonesArray* tpchits=0;
358
359     if (TOF) tofhits = TOF->Hits();
360     if (TPC) tpchits = TPC->Hits();
361
362     TTree *TH = gAlice->TreeH();
363     if (!TH) return;
364     Int_t ntracks    = (Int_t) (TH->GetEntries()); // primary tracks
365     cout << "number of primary tracked tracks in current event " << ntracks << endl; // number of primary tracked tracks
366     // array declaration and initialization
367     // TOF arrays
368     //    Int_t mapPixels[AliTOFConstants::fgkNSectors*AliTOFConstants::fgkNPlates][AliTOFConstants::fgkNStripC][AliTOFConstants::fgkNpadZ*AliTOFConstants::fgkNpadX];
369
370     Int_t *** mapPixels = new Int_t**[AliTOFConstants::fgkNSectors*AliTOFConstants::fgkNPlates];
371     for (Int_t i=0; i<AliTOFConstants::fgkNSectors*AliTOFConstants::fgkNPlates; i++) mapPixels[i] = new Int_t*[AliTOFConstants::fgkNStripC];
372     for (Int_t i=0; i<AliTOFConstants::fgkNSectors*AliTOFConstants::fgkNPlates; i++) {
373       for (Int_t j=0; j<AliTOFConstants::fgkNStripC; j++) {
374         mapPixels[i][j]= new Int_t[AliTOFConstants::fgkNpadZ*AliTOFConstants::fgkNpadX];
375       }
376     }
377
378
379     // initializing the previous array
380     for (Int_t i=0;i<AliTOFConstants::fgkNSectors*AliTOFConstants::fgkNPlates;i++) {
381       for (Int_t j=0;j<AliTOFConstants::fgkNStripC;j++) {
382         for (Int_t l=0;l<AliTOFConstants::fgkNpadZ*AliTOFConstants::fgkNpadX;l++) {
383           mapPixels[i][j][l]=0;
384         }
385       } 
386     }
387
388     Float_t * toftime = new Float_t toftime[fMaxAllTracks]; 
389     InitArray(toftime, fMaxAllTracks);
390     AliTOFPad* pixelArray = new AliTOFPad[fMaxPixels];
391     Int_t* iTOFpixel        = new Int_t[fMaxAllTracks];
392     InitArray(iTOFpixel   , fMaxAllTracks);
393     Int_t* kTOFhitFirst     = new Int_t[fMaxAllTracks];
394     InitArray(kTOFhitFirst, fMaxAllTracks);
395     AliTOFRecHit* hitArray  = new AliTOFRecHit[fMaxTOFHits];
396     Int_t isHitOnFiredPad=0; // index used to fill hitArray (array used to store informations
397                              // about pads that contains an hit)
398     Int_t ntotFiredPads=0;   // index used to fill array -> total number of fired pads (at least one time)
399
400     // TPC arrays
401     AliTOFTrack* trackArray = new AliTOFTrack[fMaxTracks];
402     Int_t * iparticle = new Int_t[fMaxAllTracks];
403     InitArray(iparticle,fMaxAllTracks); 
404     Int_t * iTrackPt  = new Int_t[fMaxTracks];
405     InitArray(iTrackPt, fMaxTracks);  // array 
406     Float_t * ptTrack = new Float_t[fMaxTracks];
407     InitArray( ptTrack, fMaxTracks);  // array for selected track pt  
408     Int_t   ntotTPCtracks=0; // total number of selected TPC tracks
409
410     
411     // reading TOF hits
412     if(TOF) ReadTOFHits(ntracks, TH, tofhits, mapPixels, kTOFhitFirst, pixelArray, iTOFpixel, toftime, hitArray,isHitOnFiredPad,ntotFiredPads);
413     cout << "isHitOnFiredPad " << isHitOnFiredPad << " for event " << ievent << endl;
414
415     // start debug for adding noise
416     // adding noise
417     Int_t nHitsNoNoise=isHitOnFiredPad;
418
419     
420     if(fNoise) AddNoiseFromOuter(option,mapPixels,pixelArray,hitArray,isHitOnFiredPad,ntotFiredPads);
421     cout << "ntotFiredPads after adding noise  " << ntotFiredPads   << " for event " << ievent << endl;
422     // set the hitArray distance to nearest hit
423     SetMinDistance(hitArray,nHitsNoNoise);
424
425     // these lines has to be commented till TPC will provide fPx fPy fPz 
426     // and fL in AliTPChit class
427     // reading TPC hits
428     /*
429     if(TPC) ReadTPCHits(ntracks, TH, tpchits, iTrackPt, iparticle, ptTrack, trackArray,ntotTPCtracks);
430     */
431     
432     // geometrical matching
433     if(TOF && TPC) Matching(trackArray,hitArray,mapPixels,pixelArray,kTOFhitFirst,ntotFiredPads,iTrackPt,iTOFpixel,ntotTPCtracks);
434     
435     // fill ntuple with reconstructed particles from current event
436     FillNtuple(ntracks,trackArray,hitArray,pixelArray,iTOFpixel,iparticle,toftime,ntotFiredPads,ntotTPCtracks);
437     
438
439     // free used memory
440     delete [] toftime;
441     delete [] pixelArray;
442     delete [] iTOFpixel;
443     delete [] kTOFhitFirst;
444     delete [] hitArray;
445     delete [] trackArray;
446     delete [] iparticle;
447     delete [] iTrackPt;
448     delete [] ptTrack;
449
450    for (Int_t i=0; i<AliTOFConstants::fgkNSectors*AliTOFConstants::fgkNPlates; i++) {
451       for (Int_t j=0; j<AliTOFConstants::fgkNStripC; j++) {
452         delete [] mapPixels[i][j];
453       }
454     }
455     for (Int_t i=0; i<AliTOFConstants::fgkNSectors*AliTOFConstants::fgkNPlates; i++) delete [] mapPixels[i];
456
457     delete [] mapPixels;
458
459   }//event loop
460
461   // free used memory for ftail
462   if (ftail)
463     {
464       delete ftail;
465       ftail = 0;
466     }
467
468   // writing ntuple on output file
469   foutputfile->cd();
470   //foutputntuple->Write(0,TObject::kOverwrite);
471   foutputntuple->Write();
472   foutputfile->Write();
473   foutputfile->Close();
474
475   gBenchmark->Stop("TOFReconstruction");
476   cout << "AliTOFReconstructioner:" << endl ;
477   cout << "   took " << gBenchmark->GetCpuTime("TOFReconstruction") << " seconds in order to make the reconstruction for " <<  fNevents << " events " << endl;
478   cout <<  gBenchmark->GetCpuTime("TOFReconstruction")/fNevents << " seconds per event " << endl ;
479   cout << endl ;
480   
481 }
482  
483 //__________________________________________________________________
484 void AliTOFReconstructioner::SetRecFile(char * file )
485 {
486   //
487   // Set the file name for reconstruction output 
488   //
489   if(!fRecFile.IsNull())
490     cout << "Changing destination file for TOF reconstruction from " <<(char *)fRecFile.Data() << " to " << file << endl ;
491   fRecFile=file ;
492 }
493 //__________________________________________________________________
494 void AliTOFReconstructioner::Print(Option_t* option)const
495 {
496   //
497   // Print reconstruction output file name
498   //
499   cout << "------------------- "<< GetName() << " -------------" << endl ;
500   if(fRecFile.IsNull())
501     cout << " Writing reconstructed particles to file galice.root "<< endl ;
502   else
503     cout << "    Writing reconstructed particle to file  " << (char*) fRecFile.Data() << endl ;
504
505 }
506
507 //__________________________________________________________________
508 void AliTOFReconstructioner::PrintParameters()const
509 {
510   //
511   // Print parameters used for reconstruction
512   //
513   cout << " ------------------- "<< GetName() << " -------------" << endl ;
514   cout << " Parameters used for TOF reconstruction " << endl ;
515   //  Printing the parameters
516   
517   cout << " Number of events:                        " << fNevents << endl; 
518   cout << " Recostruction from event                 "<< fFirstEvent << "  to event "<< fLastEvent << endl;
519   cout << " TOF geometry parameters                  " << endl;
520   cout << " Min. radius of the TOF (cm)              "<< AliTOFConstants::fgkrmin << endl;
521   cout << " Max. radius of the TOF (cm)              "<< AliTOFConstants::fgkrmax << endl;
522   cout << " Number of TOF geom. levels               "<< AliTOFConstants::fgkmaxtoftree<< endl;
523   cout << " Number of TOF sectors                    "<< AliTOFConstants::fgkNSectors << endl;
524   cout << " Number of TOF modules                    "<< AliTOFConstants::fgkNPlates << endl;
525   cout << " Max. Number of strips in a module        "<< AliTOFConstants::fgkNStripC << endl;
526   cout << " Number of pads per strip                 "<< AliTOFConstants::fgkNpadX*AliTOFConstants::fgkNpadZ << endl;
527   cout << " Number of strips in central module       "<< AliTOFConstants::fgkNStripA << endl;
528   cout << " Number of strips in intermediate modules "<< AliTOFConstants::fgkNStripB << endl;
529   cout << " Number of strips in outer modules        "<< AliTOFConstants::fgkNStripC << endl;
530   cout << " Number of MRPC in x strip direction      "<< AliTOFConstants::fgkNpadX<< endl;
531   cout << " Size of MRPC (cm) along X                "<< AliTOFConstants::fgkXPad<< endl;
532   cout << " Number of MRPC in z strip direction      "<< AliTOFConstants::fgkNpadZ<<endl;
533   cout << " Size of MRPC (cm) along Z                "<< AliTOFConstants::fgkZPad<<endl;
534   cout << " Module Lengths (cm)" << endl;
535   cout << " A Module: "<< AliTOFConstants::fgkzlenA<< "  B Modules: "<< AliTOFConstants::fgkzlenB<< "  C Modules: "<< AliTOFConstants::fgkzlenC<< endl;
536   cout << " Inner radius of the TOF detector (cm): "<<AliTOFConstants::fgkrmin << endl;
537   cout << " Outer radius of the TOF detector (cm): "<<AliTOFConstants::fgkrmax << endl;
538   cout << " Max. half z-size of TOF (cm)         : "<<AliTOFConstants::fgkMaxhZtof << endl;
539   cout << " TOF Pad parameters   " << endl;
540   cout << " Time Resolution (ns) "<< fTimeResolution <<" Pad Efficiency: "<< fpadefficiency << endl;
541   cout << " Edge Effect option:  "<<  fEdgeEffect<< endl;
542
543   cout << " Boundary Effect Simulation Parameters " << endl;
544   cout << " Hparameter: "<< fHparameter<<"  H2parameter:"<< fH2parameter <<"  Kparameter:"<< fKparameter<<"  K2parameter: "<< fK2parameter << endl;
545   cout << " Efficiency in the central region of the pad: "<< fEffCenter << endl;
546   cout << " Efficiency at the boundary region of the pad: "<< fEffBoundary << endl;
547   cout << " Efficiency value at H2parameter "<< fEff2Boundary << endl;
548   cout << " Efficiency value at K2parameter "<< fEff3Boundary << endl;
549   cout << " Resolution (ps) in the central region of the pad: "<< fResCenter << endl;
550   cout << " Resolution (ps) at the boundary of the pad      : "<< fResBoundary << endl;
551   cout << " Slope (ps/K) for neighbouring pad               : "<< fResSlope <<endl;
552   cout << " Time walk (ps) in the central region of the pad : "<< fTimeWalkCenter << endl;
553   cout << " Time walk (ps) at the boundary of the pad       : "<< fTimeWalkBoundary<< endl;
554   cout << " Slope (ps/K) for neighbouring pad               : "<< fTimeWalkSlope<<endl;
555   cout << " Pulse Heigth Simulation Parameters " << endl;
556   cout << " Flag for delay due to the PulseHeightEffect: "<< fTimeDelayFlag <<endl;
557   cout << " Pulse Height Slope                           : "<< fPulseHeightSlope<<endl;
558   cout << " Time Delay Slope                             : "<< fTimeDelaySlope<<endl;
559   cout << " Minimum charge amount which could be induced : "<< fMinimumCharge<<endl;
560   cout << " Smearing in charge in (q1/q2) vs x plot      : "<< fChargeSmearing<<endl;
561   cout << " Smearing in log of charge ratio              : "<< fLogChargeSmearing<<endl;
562   cout << " Smearing in time in time vs log(q1/q2) plot  : "<< fTimeSmearing<<endl;
563   cout << " Flag for average time                        : "<< fAverageTimeFlag<<endl;
564   cout << " Charge factor flag for matching              : "<< fChargeFactorForMatching<<endl;
565   cout << " Edge tails option                            : "<< fEdgeTails << endl;
566   cout << " TPC tracking  parameters " << endl;
567   cout << " TPC tracking efficiency                      : "<< fTrackingEfficiency<< endl;
568   cout << " Sigma vs momentum dependency flag            : "<< fSigmavsp << endl;
569   cout << " Space uncertainties (cm). sigma(z) (cm): "<< fSigmaZ << " sigma(R(phi)) (cm): "<< fSigmarphi << endl;
570   cout << " Momentum uncertainties.   sigma(delta(P)/P): "<< fSigmap <<" sigma(phi) (rad): "<< fSigmaPhi <<" sigma(theta) (rad): "<< fSigmaTheta << endl;   
571   cout << " Parameters for additional noise hits " << endl;
572   cout << " Number of noise hits : " << fNoise <<" Slope parameter (ns) in the time distribution: " << fNoiseSlope << endl;
573   cout << " Mean TOF for noise from outer regions (ns)" <<  fNoiseMeanTof << endl;
574   cout << " Physical parameters " << endl;
575   cout << " Magnetic Field (tesla)                   : "<< fField <<endl;
576   cout << " Radiation length of the outer wall of TPC: "<< fRadLenTPC << endl;
577   cout << " (TPC tracks)-(TOF pads) matching parameters " << endl;
578   cout << " TRD Correction flag       : "<< fCorrectionTRD <<endl;
579   cout << " Number of the last TPC row: "<< fLastTPCRow <<" Vertex radius (cm) for selected tracks: "<<fRadiusvtxBound<<endl;
580   cout << " Max. number of test tracks: "<<fMaxTestTracks << endl;
581   cout << " Space step (cm)           : "<< fStep <<endl;
582   cout << " Matching style option     : "<< fMatchingStyle <<endl;
583   cout << " Array parameters " << endl;
584   cout << " Max.number of pads involved in the matching procedure: "<< fMaxPixels << endl;
585   cout << " Max.number of TOF hits per event                     : "<< fMaxTOFHits<< endl;
586   cout << " Max.number of tracks selected for matching           : "<< fMaxTracks << endl;
587   cout << " Max.number of all tracks including the neutral ones  : "<< fMaxAllTracks<< endl;
588   cout << " Debug Flag                                           : "<< fdbg << endl;
589   cout << " Cut on momentum for selecting tracks                 : "<< fPBound << endl;
590   
591 }
592
593 //__________________________________________________________________
594 void AliTOFReconstructioner::IsInsideThePad(TGeant3 *g3, Float_t x, Float_t y, Float_t z, Int_t *nGeom, Float_t& zPad, Float_t& xPad) 
595 {
596   //   input: x,y,z - coordinates of a hit
597   //   output: array  nGeom[]
598   //          nGeom[0] - the TOF sector number, 1,2,...,18 along azimuthal direction starting from -90 deg.!!!
599   //          nGeom[1] - the TOF module number, 1,2,3,4,5=C,B,A,B,C along z-direction
600   //          nGeom[2] - the TOF strip  number, 1,2,... along z-direction
601   //          nGeom[3] - the TOF padz  number,  1,2=NPZ across a strip
602   //          nGeom[4] - the TOF padx  number,  1,2,...,48=NPX along a strip
603   //          zPad, xPad - coordinates of the hit in the pad frame
604   //  numbering is adopted for the version 3.05 of AliRoot
605   //  example:
606   //   from Hits: sec,pla,str,padz,padx=4,2,14,2,35
607   //  Vol. n.0: ALIC, copy number 1
608   //  Vol. n.1: B077, copy number 1
609   //  Vol. n.2: B074, copy number 5
610   //  Vol. n.3: BTO2, copy number 1
611   //  Vol. n.4: FTOB, copy number 2
612   //  Vol. n.5: FLTB, copy number 0
613   //  Vol. n.6: FSTR, copy number 14
614   //  Vol. n.7: FSEN, copy number 0
615   //  Vol. n.8: FSEZ, copy number 2
616   //  Vol. n.9: FSEX, copy number 35
617   //  Vol. n.10: FPAD, copy number 0
618
619
620   Float_t xTOF[3];
621   Int_t sector=0,module=0,strip=0,padz=0,padx=0;
622   Int_t i,numed,nLevel,copyNumber;
623   Gcvolu_t* gcvolu;
624   char name[5];
625   name[4]=0;
626   
627   for (i=0; i<AliTOFConstants::fgkmaxtoftree; i++) nGeom[i]=0;
628   zPad=100.;
629   xPad=100.;
630   
631   xTOF[0]=x;
632   xTOF[1]=y;
633   xTOF[2]=z;
634   
635   g3->Gmedia(xTOF, numed);
636   gcvolu=g3->Gcvolu();
637   nLevel=gcvolu->nlevel;
638   if(fdbg) {
639     for (Int_t i=0; i<nLevel; i++) {
640       strncpy(name,(char*) (&gcvolu->names[i]),4);
641       cout<<"Vol. n."<<i<<": "<<name<<", copy number "<<gcvolu->number[i]<<endl;
642     }
643   }
644   if(nLevel>=2) {
645     // sector type name: B071(1,2,...,10),B074(1,2,3,4,5-PHOS),B075(1,2,3-RICH)
646     strncpy(name,(char*) (&gcvolu->names[2]),4);
647     // volume copy: 1,2,...,10 for B071, 1,2,3,4,5 for B074, 1,2,3 for B075
648     copyNumber=gcvolu->number[2];
649    if(!strcmp(name,"B071")) {
650      if (copyNumber>=6 && copyNumber<=8) {
651        sector=copyNumber+10;
652      } else if (copyNumber>=1 && copyNumber<=5){
653        sector=copyNumber+7;
654      } else {
655        sector=copyNumber-8;
656      }
657    } else if(!strcmp(name,"B075")) {
658      sector=copyNumber+12;
659    } else if(!strcmp(name,"B074")) {
660      if (copyNumber>=1 && copyNumber<=3){
661        sector=copyNumber+4;
662      } else {
663        sector=copyNumber-1;
664      }
665    }
666   }
667   if(sector) {
668     nGeom[0]=sector;
669     if(nLevel>=4) {
670       // we'll use the module value in z-direction:
671       //                                    1    2    3    4    5
672       // the module order in z-direction: FTOC,FTOB,FTOA,FTOB,FTOC
673       // the module copy:                   2    2    0    1    1
674       // module type name: FTOA, FTOB, FTOC
675       strncpy(name,(char*) (&gcvolu->names[4]),4);
676       // module copy:  
677       copyNumber=gcvolu->number[4];
678       if(!strcmp(name,"FTOC")) {
679         if (copyNumber==2) {
680           module=1;
681         } else {
682           module=5;
683         }
684       } else if(!strcmp(name,"FTOB")) {
685         if (copyNumber==2) {
686           module=2;
687         } else {
688           module=4;
689         }
690       } else if(!strcmp(name,"FTOA")) {
691         module=3;
692       }
693     }
694   }
695   
696   if(module) {
697     nGeom[1]=module;
698     if(nLevel>=6) {
699       // strip type name: FSTR
700       strncpy(name,(char*) (&gcvolu->names[6]),4);
701       // strip copy:  
702       copyNumber=gcvolu->number[6];
703       if(!strcmp(name,"FSTR")) strip=copyNumber; 
704     }
705   }
706   
707   if(strip) {
708     nGeom[2]=strip;
709     if(nLevel>=8) {
710       // padz type name: FSEZ
711       strncpy(name,(char*) (&gcvolu->names[8]),4);
712       // padz copy:  
713       copyNumber=gcvolu->number[8];
714       if(!strcmp(name,"FSEZ")) padz=copyNumber; 
715     }
716   }
717   if(padz) {
718     nGeom[3]=padz;
719     if(nLevel>=9) {
720       // padx type name: FSEX
721       strncpy(name,(char*) (&gcvolu->names[9]),4);
722       // padx copy:  
723       copyNumber=gcvolu->number[9];
724       if(!strcmp(name,"FSEX")) padx=copyNumber; 
725     }
726   }
727   
728   if(padx) {
729     nGeom[4]=padx;
730     zPad=gcvolu->glx[2];  // check here
731     xPad=gcvolu->glx[0];  // check here
732   }
733   
734   //   printf(" nGeom[0,1,2,3,4]=%i,%i,%i,%i,%i\n",nGeom[0],nGeom[1],nGeom[2],nGeom[3],nGeom[4]); 
735 }
736
737 //__________________________________________________________________
738 void AliTOFReconstructioner::EpMulScatt(Float_t& px, Float_t& py, Float_t& pz, Float_t& p, Float_t& theta)
739 {
740   //   Momentum p  - before mult.scat.
741   //   Momentum p2 - after mult.scat.
742   //   THE0 - r.m.s. of deviation angle in plane
743   //           (see RPP'96: Phys.Rev.D54 (1996) 134)
744   
745   Float_t pt,thex,they,tantx,tanty,p2px,p2py,p2pz,costhe,sinthe,cospsi,sinpsi,p2x,p2y,p2z,p2,g;
746   
747   pt=TMath::Sqrt(px*px+py*py);
748   //   angles for p in the ' frame with Z'along p
749   if(fMatchingStyle==1) {
750     thex=theta*gRandom->Gaus();
751     they=theta*gRandom->Gaus();
752   } else {
753     thex=3*(-theta+2*theta*gRandom->Rndm());
754     they=3*(-theta+2*theta*gRandom->Rndm());
755   }
756   tantx=TMath::Tan(thex);
757   tanty=TMath::Tan(they);
758   
759   //   p2p - p2 in the ' frame
760   p2pz=p/TMath::Sqrt(1.+tantx*tantx+tanty*tanty);
761   p2py=p2pz*tanty;
762   p2px=p2pz*tantx;
763   //   choose X'so that PHI=0 (see Il'in, Pozdnyak Analiticheskaya geometriya, 1968, c.88
764   //   for Euler angles PSI, THETA (PHI=0)
765   costhe=pz/p;
766   sinthe=pt/p;
767   cospsi=-py/pt;
768   sinpsi=px/pt;
769   //
770   g=p2py*costhe-p2pz*sinthe;
771   p2x=p2px*cospsi-g*sinpsi;
772   p2y=p2px*sinpsi+g*cospsi;
773   p2z=p2py*sinthe+p2pz*costhe;
774   p2=TMath::Sqrt(p2x*p2x+p2y*p2y+p2z*p2z);
775   
776   //   Test angle
777   g=(px*p2x+py*p2y+pz*p2z)/(p*p2);
778   if(g>1) g=1;
779   theta=TMath::ACos(g);
780   px=p2x;
781   py=p2y;
782   pz=p2z;
783   p=p2;
784   
785 }
786
787 // std border effect algorithm
788 //__________________________________________________________________
789 void AliTOFReconstructioner::BorderEffect(Float_t z0, Float_t x0, Float_t geantTime, Int_t& nActivatedPads, Int_t& nFiredPads, Bool_t* isFired, Int_t* nPlace, Float_t* qInduced, Float_t* tofTime, Float_t& averageTime)
790 {
791   // Input:  z0, x0 - hit position in the strip system (0,0 - center of the strip), cm
792   //         geantTime - time generated by Geant, ns
793   // Output: nActivatedPads - the number of pads activated by the hit (1 || 2 || 4)
794   //         nFiredPads - the number of pads fired (really activated) by the hit (nFiredPads <= nActivatedPads)
795   //         qInduced[iPad]- charge induced on pad, arb. units
796   //                         this array is initialized at zero by the caller
797   //         tofAfterSimul[iPad] - time calculated with edge effect algorithm, ns
798   //                                   this array is initialized at zero by the caller
799   //         averageTime - time given by pad hited by the Geant track taking into account the times (weighted) given by the pads fired for edge effect also.
800   //                       The weight is given by the qInduced[iPad]/qCenterPad
801   //                                   this variable is initialized at zero by the caller
802   //         nPlace[iPad] - the number of the pad place, iPad = 0, 1, 2, 3
803   //                                   this variable is initialized at zero by the caller
804   //
805   // Description of used variables:
806   //         eff[iPad] - efficiency of the pad
807   //         res[iPad] - resolution of the pad, ns
808   //         timeWalk[iPad] - time walk of the pad, ns
809   //         timeDelay[iPad] - time delay for neighbouring pad to hited pad, ns
810   //         PadId[iPad] - Pad Identifier
811   //                    E | F    -->   PadId[iPad] = 5 | 6
812   //                    A | B    -->   PadId[iPad] = 1 | 2
813   //                    C | D    -->   PadId[iPad] = 3 | 4
814   //         nTail[iPad] - the tail number, = 1 for tailA, = 2 for tailB
815   //         qCenterPad - charge extimated for each pad, arb. units
816   //         weightsSum - sum of weights extimated for each pad fired, arb. units
817   
818   const Float_t kSigmaForTail[2] = {AliTOFConstants::fgkSigmaForTail1,AliTOFConstants::fgkSigmaForTail2}; //for tail                                                   
819   Int_t iz = 0, ix = 0;
820   Float_t dX = 0., dZ = 0., x = 0., z = 0.;
821   Float_t h = fHparameter, h2 = fH2parameter, k = fKparameter, k2 = fK2parameter;
822   Float_t effX = 0., effZ = 0., resX = 0., resZ = 0., timeWalkX = 0., timeWalkZ = 0.;
823   Float_t logOfqInd = 0.;
824   Float_t weightsSum = 0.;
825   Int_t nTail[4]  = {0,0,0,0};
826   Int_t padId[4]  = {0,0,0,0};
827   Float_t eff[4]  = {0.,0.,0.,0.};
828   Float_t res[4]  = {0.,0.,0.,0.};
829   //  Float_t qCenterPad = fMinimumCharge * fMinimumCharge;
830   Float_t qCenterPad = 1.;
831   Float_t timeWalk[4]  = {0.,0.,0.,0.};
832   Float_t timeDelay[4] = {0.,0.,0.,0.};
833   
834   nActivatedPads = 0;
835   nFiredPads = 0;
836   
837   (z0 <= 0) ? iz = 0 : iz = 1;
838   dZ = z0 + (0.5 * AliTOFConstants::fgkNpadZ - iz - 0.5) * AliTOFConstants::fgkZPad; // hit position in the pad frame, (0,0) - center of the pad
839   z = 0.5 * AliTOFConstants::fgkZPad - TMath::Abs(dZ);                               // variable for eff., res. and timeWalk. functions
840   iz++;                                                                              // z row: 1, ..., AliTOFConstants::fgkNpadZ = 2
841   ix = (Int_t)((x0 + 0.5 * AliTOFConstants::fgkNpadX * AliTOFConstants::fgkXPad) / AliTOFConstants::fgkXPad);
842   dX = x0 + (0.5 * AliTOFConstants::fgkNpadX - ix - 0.5) * AliTOFConstants::fgkXPad; // hit position in the pad frame, (0,0) - center of the pad
843   x = 0.5 * AliTOFConstants::fgkXPad - TMath::Abs(dX);                               // variable for eff., res. and timeWalk. functions;
844   ix++;                                                                              // x row: 1, ..., AliTOFConstants::fgkNpadX = 48
845   
846   ////// Pad A:
847   nActivatedPads++;
848   nPlace[nActivatedPads-1] = (iz - 1) * AliTOFConstants::fgkNpadX + ix;
849   qInduced[nActivatedPads-1] = qCenterPad;
850   padId[nActivatedPads-1] = 1;
851   
852   if (fEdgeEffect == 0) {
853     eff[nActivatedPads-1] = fEffCenter;
854     if (gRandom->Rndm() < eff[nActivatedPads-1]) {
855       nFiredPads = 1;
856       res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + fResCenter * fResCenter); // 10400=30^2+20^2+40^2+50^2+50^2+50^2  ns;
857       isFired[nActivatedPads-1] = kTRUE;
858       tofTime[nActivatedPads-1] = gRandom->Gaus(geantTime + fTimeWalkCenter, res[0]);
859       averageTime = tofTime[nActivatedPads-1];
860     }
861   } else {
862      
863     if(z < h) {
864       if(z < h2) {
865         effZ = fEffBoundary + (fEff2Boundary - fEffBoundary) * z / h2;
866       } else {
867         effZ = fEff2Boundary + (fEffCenter - fEff2Boundary) * (z - h2) / (h - h2);
868       }
869       resZ = fResBoundary + (fResCenter - fResBoundary) * z / h;
870       timeWalkZ = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * z / h;
871       nTail[nActivatedPads-1] = 1;
872     } else {
873       effZ = fEffCenter;
874       resZ = fResCenter;
875       timeWalkZ = fTimeWalkCenter;
876     }
877     
878     if(x < h) {
879       if(x < h2) {
880         effX = fEffBoundary + (fEff2Boundary - fEffBoundary) * x / h2;
881       } else {
882         effX = fEff2Boundary + (fEffCenter - fEff2Boundary) * (x - h2) / (h - h2);
883       }
884       resX = fResBoundary + (fResCenter - fResBoundary) * x / h;
885       timeWalkX = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * x / h;
886       nTail[nActivatedPads-1] = 1;
887     } else {
888       effX = fEffCenter;
889       resX = fResCenter;
890       timeWalkX = fTimeWalkCenter;
891     }
892     
893     (effZ<effX) ? eff[nActivatedPads-1] = effZ : eff[nActivatedPads-1] = effX;
894     (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2  ns
895     (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 *  timeWalkZ : timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
896
897
898     ////// Pad B:
899     if(z < k2) {
900       effZ = fEffBoundary - (fEffBoundary - fEff3Boundary) * (z / k2);
901     } else {
902       effZ = fEff3Boundary * (k - z) / (k - k2);
903     }
904     resZ = fResBoundary + fResSlope * z / k;
905     timeWalkZ = fTimeWalkBoundary + fTimeWalkSlope * z / k;
906     
907     if(z < k && z > 0) {
908       if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
909         nActivatedPads++;
910         nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX;
911         eff[nActivatedPads-1] = effZ;
912         res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns 
913         timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ; // ns
914         nTail[nActivatedPads-1] = 2;
915         if (fTimeDelayFlag) {
916           //      qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
917           //      qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
918           qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * z);
919           logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
920           timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
921         } else {
922           timeDelay[nActivatedPads-1] = 0.;
923         }
924         padId[nActivatedPads-1] = 2;
925       }
926     }
927
928     
929     ////// Pad C, D, E, F:
930     if(x < k2) {
931       effX = fEffBoundary - (fEffBoundary - fEff3Boundary) * (x / k2);
932     } else {
933       effX = fEff3Boundary * (k - x) / (k - k2);
934     }
935     resX = fResBoundary + fResSlope*x/k;
936     timeWalkX = fTimeWalkBoundary + fTimeWalkSlope*x/k;
937     
938     if(x < k && x > 0) {
939       //   C:
940       if(ix > 1 && dX < 0) {
941         nActivatedPads++;
942         nPlace[nActivatedPads-1] = nPlace[0] - 1;
943         eff[nActivatedPads-1] = effX;
944         res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns 
945         timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
946         nTail[nActivatedPads-1] = 2;
947         if (fTimeDelayFlag) {
948           //      qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
949           //      qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
950           qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
951           logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
952           timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
953         } else {
954           timeDelay[nActivatedPads-1] = 0.;
955         }
956         padId[nActivatedPads-1] = 3;
957
958         //     D:
959         if(z < k && z > 0) {
960           if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
961             nActivatedPads++;
962             nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX - 1;
963             eff[nActivatedPads-1] = effX * effZ;
964             (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
965             (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
966             
967             nTail[nActivatedPads-1] = 2;
968             if (fTimeDelayFlag) {
969               if (TMath::Abs(x) < TMath::Abs(z)) {
970                 //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
971                 //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
972                 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * z);
973                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
974               } else {
975                 //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
976                 //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
977                 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
978                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
979               }
980               timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
981             } else {
982               timeDelay[nActivatedPads-1] = 0.;
983             }
984             padId[nActivatedPads-1] = 4;
985           }
986         }  // end D
987       }  // end C
988       
989       //   E:
990       if(ix < AliTOFConstants::fgkNpadX && dX > 0) {
991         nActivatedPads++;
992         nPlace[nActivatedPads-1] = nPlace[0] + 1;
993         eff[nActivatedPads-1] = effX;
994         res[nActivatedPads-1] = 0.001 * (TMath::Sqrt(10400 + resX * resX)); // ns
995         timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
996         nTail[nActivatedPads-1] = 2;
997         if (fTimeDelayFlag) {
998           //      qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
999           //      qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
1000           qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
1001           logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
1002           timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1003         } else {
1004           timeDelay[nActivatedPads-1] = 0.;
1005         }
1006         padId[nActivatedPads-1] = 5;
1007
1008
1009         //     F:
1010         if(z < k && z > 0) {
1011           if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
1012             nActivatedPads++;
1013             nPlace[nActivatedPads - 1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX + 1;
1014             eff[nActivatedPads - 1] = effX * effZ;
1015             (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
1016             (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001*timeWalkX; // ns
1017             nTail[nActivatedPads-1] = 2;
1018             if (fTimeDelayFlag) {
1019               if (TMath::Abs(x) < TMath::Abs(z)) {
1020                 //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
1021                 //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
1022                 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * z);
1023                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
1024               } else {
1025                 //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
1026                 //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
1027                 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
1028                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
1029               }
1030               timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1031             } else {
1032               timeDelay[nActivatedPads-1] = 0.;
1033             }
1034             padId[nActivatedPads-1] = 6;
1035           }
1036         }  // end F
1037       }  // end E
1038     } // end if(x < k)
1039
1040
1041     for (Int_t iPad = 0; iPad < nActivatedPads; iPad++) {
1042       if (res[iPad] < fTimeResolution) res[iPad] = fTimeResolution;
1043       if(gRandom->Rndm() < eff[iPad]) {
1044         isFired[iPad] = kTRUE;
1045         nFiredPads++;
1046         if(fEdgeTails) {
1047           if(nTail[iPad] == 0) {
1048             tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]);
1049           } else {
1050             ftail->SetParameters(res[iPad], 2. * res[iPad], kSigmaForTail[nTail[iPad]-1]);
1051             Double_t timeAB = ftail->GetRandom();
1052             tofTime[iPad] = geantTime + timeWalk[iPad] + timeDelay[iPad] + timeAB;
1053           }
1054         } else {
1055           tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]);
1056         }
1057         if (fAverageTimeFlag) {
1058           averageTime += tofTime[iPad] * qInduced[iPad];
1059           weightsSum += qInduced[iPad];
1060         } else {
1061           averageTime += tofTime[iPad];
1062           weightsSum += 1.;
1063         }
1064       }
1065     }
1066     if (weightsSum!=0) averageTime /= weightsSum;
1067   } // end else (fEdgeEffect != 0)
1068 }
1069
1070
1071 /* new algorithm (to be checked)
1072 //__________________________________________________________________
1073 void AliTOFReconstructioner::BorderEffect(Float_t z0, Float_t x0, Float_t geantTime, Int_t& nActivatedPads, Int_t& nFiredPads, Bool_t* isFired, Int_t* nPlace, Float_t* qInduced, Float_t* tofTime, Float_t& averageTime)
1074 {
1075   // Input:  z0, x0 - hit position in the strip system (0,0 - center of the strip), cm
1076   //         geantTime - time generated by Geant, ns
1077   // Output: nActivatedPads - the number of pads activated by the hit (1 || 2 || 4)
1078   //         nFiredPads - the number of pads fired (really activated) by the hit (nFiredPads <= nActivatedPads)
1079   //         qInduced[iPad]- charge induced on pad, arb. units
1080   //                         this array is initialized at zero by the caller
1081   //         tofAfterSimul[iPad] - time calculated with edge effect algorithm, ns
1082   //                                   this array is initialized at zero by the caller
1083   //         averageTime - time given by pad hited by the Geant track taking into account the times (weighted) given by the pads fired for edge effect also.
1084   //                       The weight is given by the qInduced[iPad]/qCenterPad
1085   //                                   this variable is initialized at zero by the caller
1086   //         nPlace[iPad] - the number of the pad place, iPad = 0, 1, 2, 3
1087   //                                   this variable is initialized at zero by the caller
1088   //
1089   // Description of used variables:
1090   //         eff[iPad] - efficiency of the pad
1091   //         res[iPad] - resolution of the pad, ns
1092   //         timeWalk[iPad] - time walk of the pad, ns
1093   //         timeDelay[iPad] - time delay for neighbouring pad to hited pad, ns
1094   //         PadId[iPad] - Pad Identifier
1095   //                    E | F    -->   PadId[iPad] = 5 | 6
1096   //                    A | B    -->   PadId[iPad] = 1 | 2
1097   //                    C | D    -->   PadId[iPad] = 3 | 4
1098   //         nTail[iPad] - the tail number, = 1 for tailA, = 2 for tailB
1099   //         qCenterPad - charge extimated for each pad, arb. units
1100   //         weightsSum - sum of weights extimated for each pad fired, arb. units
1101   
1102   const Float_t kSigmaForTail[2] = {AliTOFConstants::fgkSigmaForTail1,AliTOFConstants::fgkSigmaForTail2}; //for tail                                                   
1103   Int_t iz = 0, ix = 0;
1104   Float_t dX = 0., dZ = 0., x = 0., z = 0.;
1105   Float_t h = fHparameter, h2 = fH2parameter, k = fKparameter, k2 = fK2parameter;
1106   Float_t effX = 0., effZ = 0., resX = 0., resZ = 0., timeWalkX = 0., timeWalkZ = 0.;
1107   Float_t logOfqInd = 0.;
1108   Float_t weightsSum = 0.;
1109   Int_t nTail[4]  = {0,0,0,0};
1110   Int_t padId[4]  = {0,0,0,0};
1111   Float_t eff[4]  = {0.,0.,0.,0.};
1112   Float_t res[4]  = {0.,0.,0.,0.};
1113   Float_t qCenterPad = fMinimumCharge * fMinimumCharge;
1114   Float_t timeWalk[4]  = {0.,0.,0.,0.};
1115   Float_t timeDelay[4] = {0.,0.,0.,0.};
1116   
1117   nActivatedPads = 0;
1118   nFiredPads = 0;
1119   
1120   (z0 <= 0) ? iz = 0 : iz = 1;
1121   dZ = z0 + (0.5 * AliTOFConstants::fgkNpadZ - iz - 0.5) * AliTOFConstants::fgkZPad; // hit position in the pad frame, (0,0) - center of the pad
1122   z = 0.5 * AliTOFConstants::fgkZPad - TMath::Abs(dZ);                               // variable for eff., res. and timeWalk. functions
1123   iz++;                                                                              // z row: 1, ..., AliTOFConstants::fgkNpadZ = 2
1124   ix = (Int_t)((x0 + 0.5 * AliTOFConstants::fgkNpadX * AliTOFConstants::fgkXPad) / AliTOFConstants::fgkXPad);
1125   dX = x0 + (0.5 * AliTOFConstants::fgkNpadX - ix - 0.5) * AliTOFConstants::fgkXPad; // hit position in the pad frame, (0,0) - center of the pad
1126   x = 0.5 * AliTOFConstants::fgkXPad - TMath::Abs(dX);                               // variable for eff., res. and timeWalk. functions;
1127   ix++;                                                                              // x row: 1, ..., AliTOFConstants::fgkNpadX = 48
1128   
1129   ////// Pad A:
1130   nActivatedPads++;
1131   nPlace[nActivatedPads-1] = (iz - 1) * AliTOFConstants::fgkNpadX + ix;
1132   qInduced[nActivatedPads-1] = qCenterPad;
1133   padId[nActivatedPads-1] = 1;
1134   
1135   if (fEdgeEffect == 0) {
1136     eff[nActivatedPads-1] = fEffCenter;
1137     if (gRandom->Rndm() < eff[nActivatedPads-1]) {
1138       nFiredPads = 1;
1139       res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + fResCenter * fResCenter); // 10400=30^2+20^2+40^2+50^2+50^2+50^2  ns;
1140       isFired[nActivatedPads-1] = kTRUE;
1141       tofTime[nActivatedPads-1] = gRandom->Gaus(geantTime + fTimeWalkCenter, res[0]);
1142       averageTime = tofTime[nActivatedPads-1];
1143     }
1144   } else {
1145      
1146     if(z < h) {
1147       if(z < h2) {
1148         effZ = fEffBoundary + (fEff2Boundary - fEffBoundary) * z / h2;
1149       } else {
1150         effZ = fEff2Boundary + (fEffCenter - fEff2Boundary) * (z - h2) / (h - h2);
1151       }
1152       resZ = fResBoundary + (fResCenter - fResBoundary) * z / h;
1153       timeWalkZ = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * z / h;
1154       nTail[nActivatedPads-1] = 1;
1155     } else {
1156       effZ = fEffCenter;
1157       resZ = fResCenter;
1158       timeWalkZ = fTimeWalkCenter;
1159     }
1160     
1161     if(x < h) {
1162       if(x < h2) {
1163         effX = fEffBoundary + (fEff2Boundary - fEffBoundary) * x / h2;
1164       } else {
1165         effX = fEff2Boundary + (fEffCenter - fEff2Boundary) * (x - h2) / (h - h2);
1166       }
1167       resX = fResBoundary + (fResCenter - fResBoundary) * x / h;
1168       timeWalkX = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * x / h;
1169       nTail[nActivatedPads-1] = 1;
1170     } else {
1171       effX = fEffCenter;
1172       resX = fResCenter;
1173       timeWalkX = fTimeWalkCenter;
1174     }
1175     
1176     (effZ<effX) ? eff[nActivatedPads-1] = effZ : eff[nActivatedPads-1] = effX;
1177     (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2  ns
1178     (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 *  timeWalkZ : timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
1179
1180
1181     ////// Pad B:
1182     if(z < k2) {
1183       effZ = fEffBoundary - (fEffBoundary - fEff3Boundary) * (z / k2);
1184     } else {
1185       effZ = fEff3Boundary * (k - z) / (k - k2);
1186     }
1187     resZ = fResBoundary + fResSlope * z / k;
1188     timeWalkZ = fTimeWalkBoundary + fTimeWalkSlope * z / k;
1189     
1190     if(z < k && z > 0) {
1191       if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
1192         nActivatedPads++;
1193         nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX;
1194         eff[nActivatedPads-1] = effZ;
1195         res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns 
1196         timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ; // ns
1197         nTail[nActivatedPads-1] = 2;
1198         if (fTimeDelayFlag) {
1199           qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
1200           qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
1201           logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
1202           timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1203         } else {
1204           timeDelay[nActivatedPads-1] = 0.;
1205         }
1206         padId[nActivatedPads-1] = 2;
1207       }
1208     }
1209
1210     
1211     ////// Pad C, D, E, F:
1212     if(x < k2) {
1213       effX = fEffBoundary - (fEffBoundary - fEff3Boundary) * (x / k2);
1214     } else {
1215       effX = fEff3Boundary * (k - x) / (k - k2);
1216     }
1217     resX = fResBoundary + fResSlope*x/k;
1218     timeWalkX = fTimeWalkBoundary + fTimeWalkSlope*x/k;
1219     
1220     if(x < k && x > 0) {
1221       //   C:
1222       if(ix > 1 && dX < 0) {
1223         nActivatedPads++;
1224         nPlace[nActivatedPads-1] = nPlace[0] - 1;
1225         eff[nActivatedPads-1] = effX;
1226         res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns 
1227         timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
1228         nTail[nActivatedPads-1] = 2;
1229         if (fTimeDelayFlag) {
1230           qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
1231           qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
1232           logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
1233           timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1234         } else {
1235           timeDelay[nActivatedPads-1] = 0.;
1236         }
1237         padId[nActivatedPads-1] = 3;
1238
1239         //     D:
1240         if(z < k && z > 0) {
1241           if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
1242             nActivatedPads++;
1243             nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX - 1;
1244             eff[nActivatedPads-1] = effX * effZ;
1245             (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
1246             (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
1247             
1248             nTail[nActivatedPads-1] = 2;
1249             if (fTimeDelayFlag) {
1250               if (TMath::Abs(x) < TMath::Abs(z)) {
1251                 qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
1252                 qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
1253                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
1254               } else {
1255                 qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
1256                 qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
1257                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
1258               }
1259               timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1260             } else {
1261               timeDelay[nActivatedPads-1] = 0.;
1262             }
1263             padId[nActivatedPads-1] = 4;
1264           }
1265         }  // end D
1266       }  // end C
1267       
1268       //   E:
1269       if(ix < AliTOFConstants::fgkNpadX && dX > 0) {
1270         nActivatedPads++;
1271         nPlace[nActivatedPads-1] = nPlace[0] + 1;
1272         eff[nActivatedPads-1] = effX;
1273         res[nActivatedPads-1] = 0.001 * (TMath::Sqrt(10400 + resX * resX)); // ns
1274         timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
1275         nTail[nActivatedPads-1] = 2;
1276         if (fTimeDelayFlag) {
1277           qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
1278           qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
1279           logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
1280           timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1281         } else {
1282           timeDelay[nActivatedPads-1] = 0.;
1283         }
1284         padId[nActivatedPads-1] = 5;
1285
1286
1287         //     F:
1288         if(z < k && z > 0) {
1289           if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
1290             nActivatedPads++;
1291             nPlace[nActivatedPads - 1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX + 1;
1292             eff[nActivatedPads - 1] = effX * effZ;
1293             (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
1294             (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001*timeWalkX; // ns
1295             nTail[nActivatedPads-1] = 2;
1296             if (fTimeDelayFlag) {
1297               if (TMath::Abs(x) < TMath::Abs(z)) {
1298                 qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
1299                 qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
1300                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
1301               } else {
1302                 qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
1303                 qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
1304                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
1305               }
1306               timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1307             } else {
1308               timeDelay[nActivatedPads-1] = 0.;
1309             }
1310             padId[nActivatedPads-1] = 6;
1311           }
1312         }  // end F
1313       }  // end E
1314     } // end if(x < k)
1315
1316
1317     for (Int_t iPad = 0; iPad < nActivatedPads; iPad++) {
1318       if (res[iPad] < fTimeResolution) res[iPad] = fTimeResolution;
1319       if(gRandom->Rndm() < eff[iPad]) {
1320         isFired[iPad] = kTRUE;
1321         nFiredPads++;
1322         if(fEdgeTails) {
1323           if(nTail[iPad] == 0) {
1324             tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]);
1325           } else {
1326             ftail->SetParameters(res[iPad], 2. * res[iPad], kSigmaForTail[nTail[iPad]-1]);
1327             Double_t timeAB = ftail->GetRandom();
1328             tofTime[iPad] = geantTime + timeWalk[iPad] + timeDelay[iPad] + timeAB;
1329           }
1330         } else {
1331           tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]);
1332         }
1333         if (fAverageTimeFlag) {
1334           averageTime += tofTime[iPad] * qInduced[iPad];
1335           weightsSum += qInduced[iPad];
1336         } else {
1337           averageTime += tofTime[iPad];
1338           weightsSum += 1.;
1339         }
1340       }
1341     }
1342     if (weightsSum!=0) averageTime /= weightsSum;
1343
1344   } // end else (fEdgeEffect != 0)
1345   
1346   //cout << "timedelay " << timeDelay[0] << endl;
1347   //cout << "timedelay " << timeDelay[1] << endl;
1348   //cout << "timedelay " << timeDelay[2] << endl;
1349   //cout << "timedelay " << timeDelay[3] << endl;
1350   
1351 }
1352 */
1353
1354
1355 //__________________________________________________________________
1356 Int_t AliTOFReconstructioner::PDGtoGeantCode(Int_t pdgcode) 
1357 {
1358   //
1359   // Gives the GEANT code from KF code of LUND JETSET
1360   //
1361   Int_t geantCode=0; // default value
1362   switch (pdgcode) {
1363   case 22:
1364     geantCode=1;            // GAMMA
1365     break ;
1366   case -11:
1367     geantCode=2;            // E+
1368     break ;
1369   case 11:
1370     geantCode=3;            // E-
1371     break ;
1372   case 12:
1373     geantCode=4;            // NUE
1374     break ;
1375   case 14:
1376     geantCode=4;            // NUMU
1377     break ;
1378   case -13:
1379     geantCode=5;            // MU+
1380     break ;
1381   case 13:
1382     geantCode=6;            // MU-
1383     break ;
1384   case 111:
1385     geantCode=7;            // PI0
1386     break ;
1387   case 211:
1388     geantCode=8;            // PI+
1389     break ;
1390   case -211:
1391     geantCode=9;            // PI-
1392     break ;
1393   case 130:
1394     geantCode=10;           // K_L0
1395     break ;
1396   case 321:
1397     geantCode=11;           // K+
1398     break ;
1399   case -321:
1400     geantCode=12;           // K-
1401     break ;
1402   case 2112:
1403     geantCode=13;           // N0
1404     break ;
1405   case 2212:
1406     geantCode=14;           // P+
1407     break ;
1408   case -2212:
1409     geantCode=15;           // P~-
1410     break ;
1411   case 310:
1412     geantCode=16;           // K_S0
1413     break ;
1414   case 221:
1415     geantCode=17;           // ETA
1416     break ;
1417   case 3122:
1418     geantCode=18;           // LAMBDA0
1419     break ;
1420   case 3222:
1421     geantCode=19;           // SIGMA+
1422     break ;
1423   case 3212:
1424     geantCode=20;           // SIGMA0
1425     break ;
1426   case 3112:
1427     geantCode=21;           // SIGMA-
1428     break ;
1429   case 3322:
1430     geantCode=22;           // XI0
1431     break ;
1432   case 3312:
1433     geantCode=23;           // XI-
1434     break ;
1435   case 3334:
1436     geantCode=24;           // OMEGA-
1437     break ;
1438   case -2112:
1439     geantCode=25;           // N~0
1440     break ;
1441   case -3122:
1442     geantCode=26;           // LAMBDA~0
1443     break ;
1444   case -3112:
1445     geantCode=27;           // SIGMA~+
1446     break ;
1447   case -3212:
1448     geantCode=28;           // SIGMA~0
1449     break ;
1450   case -3222:
1451     geantCode=29;           // SIGMA~-
1452     break ;
1453   case -3322:
1454     geantCode=30;           // XI~0
1455     break ;
1456   case -3312:
1457     geantCode=31;           // XI~+
1458     break ;
1459   case -3334:
1460     geantCode=32;           // OMEGA~+
1461     break ;
1462   case 223:
1463     geantCode=33;           // OMEGA(782)
1464     break ;
1465   case 333:
1466     geantCode=34;           // PHI(1020)
1467     break ;
1468   case 411:
1469     geantCode=35;           // D+
1470     break ;
1471   case -411:
1472     geantCode=36;           // D-
1473     break ;
1474   case 421:
1475     geantCode=37;           // D0
1476     break ;
1477   case -421:
1478     geantCode=38;           // D~0
1479     break ;
1480   case 431:
1481     geantCode=39;           // D_S+
1482     break ;
1483   case -431:
1484     geantCode=40;           // D_S~-
1485     break ;
1486   case 4122:
1487     geantCode=41;           // LAMBDA_C+
1488     break ;
1489   case 213:
1490     geantCode=42;           // RHP(770)+
1491     break ;
1492   case -213:
1493     geantCode=43;           // RHO(770)-
1494     break ;
1495   case 113:
1496     geantCode=44;           // RHO(770)0
1497     break ;
1498   default:
1499     geantCode=45;
1500     break;
1501   }
1502
1503   return geantCode;
1504 }
1505
1506 //__________________________________________________________________
1507 Bool_t AliTOFReconstructioner::operator==( AliTOFReconstructioner const & tofrec)const
1508 {
1509   // Equal operator.
1510   // Reconstructioners are equal if their parameters are equal
1511
1512   // split the member variables in analogous categories
1513   
1514   // time resolution and edge effect parameters
1515   Bool_t dummy0=(fTimeResolution==tofrec.fTimeResolution)&&(fpadefficiency==tofrec.fpadefficiency)&&(fEdgeEffect==tofrec.fEdgeEffect)&&(fEdgeTails==tofrec.fEdgeTails)&&(fHparameter==tofrec.fHparameter)&&(fH2parameter==tofrec.fH2parameter)&&(fKparameter==tofrec.fKparameter)&&(fK2parameter==tofrec.fK2parameter);
1516   
1517   // pad efficiency parameters
1518   Bool_t dummy1=(fEffCenter==tofrec.fEffCenter)&&(fEffBoundary==tofrec.fEffBoundary)&&(fEff2Boundary==tofrec.fEff2Boundary)&&(fEff3Boundary==tofrec.fEff3Boundary)&&(fResCenter==tofrec.fResCenter)&&(fResBoundary==tofrec.fResBoundary)&&(fResSlope==tofrec.fResSlope);
1519
1520   // time walk parameters
1521   Bool_t dummy2=(fTimeWalkCenter==tofrec.fTimeWalkCenter)&&(fTimeWalkBoundary==tofrec.fTimeWalkBoundary)&&(fTimeWalkSlope==tofrec.fTimeWalkSlope)&&(fTimeDelayFlag==tofrec.fTimeDelayFlag)&&(fPulseHeightSlope==tofrec.fPulseHeightSlope)&&(fTimeDelaySlope==tofrec.fTimeDelaySlope);
1522
1523   // ADC-TDC correlation parameters
1524   Bool_t dummy3=(fMinimumCharge==tofrec.fMinimumCharge)&&(fChargeSmearing==tofrec.fChargeSmearing )&&(fLogChargeSmearing==tofrec.fLogChargeSmearing )&&(fTimeSmearing==tofrec.fTimeSmearing )&&(fAverageTimeFlag==tofrec.fAverageTimeFlag)&&(fChargeFactorForMatching==tofrec.fChargeFactorForMatching)&&(fMatchingStyle==tofrec.fMatchingStyle);
1525   
1526   Bool_t dummy4=(fTrackingEfficiency==tofrec.fTrackingEfficiency)&&(fSigmavsp==tofrec.fSigmavsp)&&(fSigmaZ==tofrec.fSigmaZ)&&(fSigmarphi==tofrec.fSigmarphi)&&(fSigmap==tofrec.fSigmap)&&(fSigmaPhi==tofrec.fSigmaPhi)&&(fSigmaTheta==tofrec.fSigmaTheta)&&(fNoise==tofrec.fNoise)&&(fNoiseSlope==tofrec.fNoiseSlope)&&(fField==tofrec.fField)&&(fRadLenTPC==tofrec.fRadLenTPC)&&(fCorrectionTRD==tofrec.fCorrectionTRD)&&(fLastTPCRow==tofrec.fLastTPCRow)&&(fRadiusvtxBound==tofrec.fRadiusvtxBound)&&(fMaxTestTracks==tofrec.fMaxTestTracks)&&(fStep==tofrec.fStep)&&(fMaxPixels==tofrec.fMaxPixels)&&(fMaxAllTracks==tofrec.fMaxAllTracks)&&(fMaxTracks==tofrec.fMaxTracks)&&(fMaxTOFHits==tofrec.fMaxTOFHits)&&(fPBound==tofrec.fPBound);
1527
1528   if( dummy0 && dummy1 && dummy2 && dummy3 && dummy4)
1529     return kTRUE ;
1530   else
1531     return kFALSE ;
1532
1533 }
1534 //____________________________________________________________________________ 
1535 void AliTOFReconstructioner::UseHitsFrom(const char * filename)
1536 {
1537   SetTitle(filename) ; 
1538 }
1539
1540 //____________________________________________________________________________ 
1541 void AliTOFReconstructioner::InitArray(Float_t array[], Int_t nlocations)
1542 {
1543   //
1544   // Initialize the array of Float_t
1545   // 
1546   for (Int_t i = 0; i < nlocations; i++) {
1547     array[i]=0.;
1548   }                             // end loop
1549
1550 }
1551
1552 //____________________________________________________________________________ 
1553 void AliTOFReconstructioner::InitArray(Int_t array[], Int_t nlocations)
1554 {
1555   //
1556   // Initialize the array of Int_t
1557   // 
1558   for (Int_t i = 0; i < nlocations; i++) {
1559     array[i]=0;
1560   }                             // end loop
1561
1562 }
1563
1564
1565 //____________________________________________________________________________
1566 void AliTOFReconstructioner::ReadTOFHits(Int_t ntracks, TTree* treehits, TClonesArray* tofhits, Int_t ***MapPixels, Int_t* kTOFhitFirst, AliTOFPad* pixelArray , Int_t* iTOFpixel, Float_t* toftime, AliTOFRecHit* hitArray, Int_t& isHitOnFiredPad, Int_t& ipixel)
1567 {
1568   //
1569   // Read TOF hits for the current event and fill arrays
1570   // 
1571   // Start loop on primary tracks in the hits containers
1572   //
1573   // Noise meaning in ReadTOFHits: we use the word 'noise' in the  
1574   // following cases
1575   // - signals produced by secondary particles
1576   // - signals produced by the next hits (out of the first) of a given track
1577   //   (both primary and secondary)
1578   // - signals produced by edge effect
1579
1580
1581   TParticle *particle;
1582   Int_t nHitOutofTofVolumes; // number of hits out of TOF GEANT volumes (it happens in very
1583                              // few cases)
1584   Int_t * npixel = new Int_t[AliTOFConstants::fgkmaxtoftree]; // array used by TOFRecon for check on TOF geometry
1585   Int_t npions=0;    // number of pions for the current event
1586   Int_t nkaons=0;    // number of kaons for the current event
1587   Int_t nprotons=0;  // number of protons for the current event
1588   Int_t nelectrons=0;// number of electrons for the current event
1589   Int_t nmuons=0;    // number of muons for the current event
1590   Float_t tofpos[3];     // TOF hit position and GEANT time
1591   Float_t zPad,xPad; 
1592   Int_t nbytes = 0;
1593   Int_t ipart, nhits=0, nHitsFromPrimaries=0;
1594   Int_t ntotalTOFhits=0; // total number of TOF hits for the current event
1595   Int_t ipartLast=-1;    // last track identifier
1596   Int_t iFirstHit;       // flag to check if the current hit is the first hit on TOF for the 
1597                          // current track
1598   Int_t iNoiseHit=0;     // flag used to tag noise hits (the noise meaning is reported in the
1599                          // header of the ReadTOFHits method)
1600   Int_t nhitWithoutNoise;// number of hits not due to noise
1601   Int_t inoise=0,inoise2=0;
1602   Int_t nMultipleSignOnSamePad=0; // number of cases where a pad is fired more than one time
1603   Int_t nPixEdge=0;      // additional pads fired due to edge effect in ReadTOFHits (local var)
1604   // array used for counting different types of primary particles
1605   Int_t particleTypeGEANT[50]={0,4,4,0,5,5,0,3,3,0,
1606                    2,2,0,1,1,0,0,0,0,0,
1607                    0,0,0,0,0,0,0,0,0,0,
1608                    0,0,0,0,0,0,0,0,0,0,
1609                    0,0,0,0,0,0,0,0,0,0};
1610   Int_t particleType,particleInTOFtype[6][3];
1611   for (Int_t i=0;i<6;i++) {
1612     for (Int_t j=0;j<3;j++) {
1613       particleInTOFtype[i][j]=0;
1614     }
1615   }
1616
1617   // speed-up the code
1618   treehits->SetBranchStatus("*",0); // switch off all branches
1619   treehits->SetBranchStatus("TOF*",1); // switch on only TOF
1620
1621   for (Int_t track=0; track<ntracks;track++) { // starting loop on primary tracks for the current event
1622
1623     gAlice->ResetHits();
1624     nbytes += treehits->GetEvent(track);
1625     nhits = tofhits->GetEntriesFast();
1626
1627     ntotalTOFhits+=nhits;      
1628
1629     // Start loop on hits connected to the current primary tracked particle
1630     // (including hits produced by secondary particles generaterd by the
1631     // current ptimary tracked particle)
1632     for (Int_t hit=0;hit<nhits;hit++) {
1633       AliTOFhit* tofHit   = (AliTOFhit*)tofhits->UncheckedAt(hit);
1634       ipart    = tofHit->GetTrack();
1635       if(ipart>=fMaxAllTracks) break;
1636       Float_t geantTime= tofHit->GetTof(); // it is given in [s]
1637       particle = (TParticle*)gAlice->Particle(ipart);
1638
1639       Int_t pdgCode=particle->GetPdgCode();
1640       // Only high momentum tracks (see fPBound value)
1641       // momentum components at vertex
1642       Float_t pxvtx = particle->Px();
1643       Float_t pyvtx = particle->Py();
1644       Float_t pzvtx = particle->Pz();
1645       Float_t pvtx = TMath::Sqrt(pxvtx*pxvtx+pyvtx*pyvtx+pzvtx*pzvtx);
1646       if(pvtx>fPBound) {
1647         
1648         if(particle->GetFirstMother() < 0) nHitsFromPrimaries++; // count primaries
1649
1650         // x and y coordinates of the particle production vertex
1651         Float_t vx = particle->Vx();
1652         Float_t vy = particle->Vy();
1653         Float_t vr = TMath::Sqrt(vx*vx+vy*vy); // cylindrical radius of the particle production vertex
1654         
1655         Float_t x = tofHit->X(); tofpos[0]=x;
1656         Float_t y = tofHit->Y(); tofpos[1]=y;
1657         Float_t z = tofHit->Z(); tofpos[2]=z;
1658         
1659         Float_t tofradius = TMath::Sqrt(x*x+y*y); // radius cilindrical coordinate of the TOF hit
1660         // momentum components (cosine) when striking the TOF
1661         Float_t pxtof = tofHit->GetPx();
1662         Float_t pytof = tofHit->GetPy();
1663         Float_t pztof = tofHit->GetPz();
1664         // scalar product indicating the direction of the particle when striking the TOF 
1665         // (>0 for outgoing particles)
1666         Float_t isGoingOut = (x*pxtof+y*pytof+z*pztof)/TMath::Sqrt(x*x+y*y+z*z);
1667         Float_t momtof = tofHit->GetMom();
1668         // now momentum components when striking the TOF
1669         pxtof *= momtof;
1670         pytof *= momtof;
1671         pztof *= momtof;
1672         particleType=particleTypeGEANT[PDGtoGeantCode(pdgCode)-1];
1673         if(particleType) {
1674           particleInTOFtype[5][2]++;
1675           particleInTOFtype[particleType-1][2]++;
1676         }
1677         iFirstHit=0;
1678         //  without noise hits
1679         
1680         if(ipart!=ipartLast) {
1681           iFirstHit=1;
1682           toftime[ipart]=geantTime;    //time [s]
1683           // tofMom[ipart]=momtof;
1684           ipartLast=ipart;
1685           if(particle->GetFirstMother() < 0) {
1686             Int_t abspdgCode=TMath::Abs(pdgCode);
1687             switch (abspdgCode) {
1688             case 211:
1689               npions++;
1690               break ;
1691             case 321:
1692               nkaons++;
1693               break ;
1694             case 2212:
1695               nprotons++;
1696               break ;
1697             case 11:
1698               nelectrons++;
1699               break ;
1700             case 13:
1701               nmuons++;
1702               break ;
1703             }
1704           }
1705           if(vr>fRadiusvtxBound) {
1706             if(particleType) { 
1707               particleInTOFtype[5][1]++;
1708               particleInTOFtype[particleType-1][1]++;
1709             }
1710             inoise++;
1711             inoise2++;
1712           } else {
1713             if(particleType) {
1714               particleInTOFtype[5][0]++;
1715               particleInTOFtype[particleType-1][0]++;
1716             }
1717           }
1718         } else {
1719           inoise++;
1720           if(particleType) {
1721             particleInTOFtype[5][1]++;
1722             particleInTOFtype[particleType-1][1]++;
1723           }
1724         } //end if(ipart!=ipartLast)
1725
1726         IsInsideThePad(fg3,x,y,z,npixel,zPad,xPad);
1727
1728         Int_t sec  = tofHit->GetSector();
1729         Int_t pla  = tofHit->GetPlate();
1730         Int_t str  = tofHit->GetStrip();
1731         if(sec!=npixel[0] || pla!=npixel[1] || str!=npixel[2]){// check on volume 
1732           cout << "sector" << sec << " npixel[0] " << npixel[0] << endl;
1733           cout << "plate " << pla << " npixel[1] " << npixel[1] << endl;
1734           cout << "strip " << str << " npixel[2] " << npixel[2] << endl;
1735         } // close check on volume
1736         
1737         Int_t padz = tofHit->GetPadz();
1738         Int_t padx = tofHit->GetPadx();
1739         Float_t Zpad = tofHit->GetDz();
1740         Float_t Xpad = tofHit->GetDx();
1741         
1742         
1743         if (npixel[4]==0){
1744           IsInsideThePad(fg3,x,y,z,npixel,zPad,xPad);
1745           if (npixel[4]==0){          
1746             nHitOutofTofVolumes++;
1747           }
1748         } else {
1749           Float_t zStrip=AliTOFConstants::fgkZPad*(padz-0.5-0.5*AliTOFConstants::fgkNpadZ)+Zpad; 
1750           if(padz!=npixel[3]) printf("            : Zpad=%f, padz=%i, npixel[3]=%i, zStrip=%f\n",Zpad,padz,npixel[3],zStrip);
1751           Float_t xStrip=AliTOFConstants::fgkXPad*(padx-0.5-0.5*AliTOFConstants::fgkNpadX)+Xpad;
1752           
1753           Int_t nPlace[4]={0,0,0,0};
1754           nPlace[0]=(padz-1)*AliTOFConstants::fgkNpadX+padx;
1755           
1756           Int_t   nActivatedPads=0;
1757           Int_t   nFiredPads=0;
1758           Bool_t  isFired[4]={kFALSE,kFALSE,kFALSE,kFALSE};
1759           Float_t tofAfterSimul[4]={0.,0.,0.,0.};
1760           Float_t qInduced[4]={0.,0.,0.,0.};
1761           Float_t averageTime=0.;    
1762
1763
1764           BorderEffect(zStrip,xStrip,geantTime*1.0e+09,nActivatedPads,nFiredPads,isFired,nPlace,qInduced,tofAfterSimul,averageTime); // simulate edge effect
1765
1766
1767           if(nFiredPads) {
1768             for(Int_t indexOfPad=0; indexOfPad<nActivatedPads; indexOfPad++) {
1769               if(isFired[indexOfPad]){// the pad has fired
1770                 if(indexOfPad==0) {// the hit belongs to a fired pad
1771                   isHitOnFiredPad++;
1772                   hitArray[isHitOnFiredPad-1].SetHit(ipart,pdgCode,tofpos,momtof,vr,iFirstHit);
1773                   iNoiseHit=0;
1774
1775                   if(vr>fRadiusvtxBound || iFirstHit==0) iNoiseHit=1;
1776                   
1777                   hitArray[isHitOnFiredPad-1].SetNoise(iNoiseHit);
1778                   if(iFirstHit) kTOFhitFirst[ipart]=isHitOnFiredPad;      
1779
1780                 }// close - the hit belongs to a fired pad
1781                 
1782                 Int_t iMapFirstIndex=AliTOFConstants::fgkNSectors*(npixel[1]-1)+npixel[0]-1;
1783                 Int_t iMapValue=MapPixels[iMapFirstIndex][npixel[2]-1][nPlace[indexOfPad]-1];
1784         
1785                 if(iMapValue==0) {
1786                   ipixel++;
1787                   if(indexOfPad) {
1788                     iNoiseHit=1;
1789                     nPixEdge++;
1790                   } else { 
1791                     iTOFpixel[ipart]=ipixel;
1792                   }
1793                   
1794                   if(ipixel>fMaxPixels){ // check on the total number of activated pads
1795                     cout << "ipixel=" << ipixel << " > fMaxPixels=" << fMaxPixels << endl;
1796                     return;
1797                   } // close check on the number of activated pads
1798                   
1799                   MapPixels[iMapFirstIndex][npixel[2]-1][nPlace[indexOfPad]-1]=ipixel;
1800                   pixelArray[ipixel-1].SetGeom(npixel[0],npixel[1],npixel[2],nPlace[indexOfPad]);
1801                   pixelArray[ipixel-1].SetTrack(ipart);
1802                   if(iNoiseHit) {
1803                     pixelArray[ipixel-1].AddState(1);
1804                   } else {
1805                     if(tofAfterSimul[indexOfPad]<0) cout << "Time of Flight after detector simulation is negative" << endl;
1806                     pixelArray[ipixel-1].AddState(10);
1807                   }
1808                   
1809                   pixelArray[ipixel-1].SetTofChargeHit(tofAfterSimul[indexOfPad],qInduced[indexOfPad],geantTime*1.0e+09,isHitOnFiredPad);
1810                 } else { //else if(iMapValue==0)
1811                   if(indexOfPad==0) iTOFpixel[ipart]=iMapValue;
1812                   nMultipleSignOnSamePad++;
1813                   
1814                   if(tofAfterSimul[indexOfPad] < pixelArray[iMapValue-1].GetRealTime() ) {
1815                     pixelArray[iMapValue-1].SetTrack(ipart);
1816                     //                   if(indexOfPad==0) pixelArray[iMapValue-1].SetTrack(ipart);
1817                     if(indexOfPad) iNoiseHit=1;
1818                     if(iNoiseHit) {
1819                       pixelArray[iMapValue-1].AddState(1);
1820                     } else {
1821                       pixelArray[iMapValue-1].AddState(10);
1822                     }
1823                     pixelArray[iMapValue-1].SetRealTime(tofAfterSimul[indexOfPad]);
1824                     pixelArray[iMapValue-1].SetGeantTime(geantTime*1.0e+09);
1825                     pixelArray[iMapValue-1].SetHit(isHitOnFiredPad);
1826                   } // close if(tofAfterSimul[indexOfPad] < pixelArray[iMapValue-1].GetTime() )
1827                 }  //end of Pixel filling
1828               }  // close if(isFired[indexOfPad])       
1829             }  //end loop on activated pads indexOfPad
1830           } // close if(nFiredPads)               
1831         }  //end of hit with npixel[3]!=0
1832       }  //high momentum tracks
1833     }  //end on TOF hits
1834   }  //end on primary tracks
1835   
1836   
1837   if(fdbg) {
1838     cout << ntotalTOFhits << " - total number of TOF hits   " << nHitsFromPrimaries << " -  primary     " <<  endl; 
1839     cout << inoise << " - noise hits, " << inoise2<< " - first crossing of a track with Rvtx>" << fRadiusvtxBound << endl;
1840     //   cout << inoise << " - noise hits (" << 100.*inoise/ihit << " %), " << inoise2 
1841     //<< " - first crossing of a track with Rvtx>" << RVTXBOUND << endl;
1842     nhitWithoutNoise=isHitOnFiredPad;
1843     
1844     cout << ipixel << " fired pixels (" << nMultipleSignOnSamePad << " multiple fired pads, " << endl;
1845     //j << " fired by noise, " << j1 << " noise+track)" <<  endl;
1846     printf(" %i additional pads are fired due to edge effect\n",nPixEdge);
1847     cout << npions <<   "   primary pions     reached TOF" << endl;
1848     cout << nkaons <<   "   primary kaons     reached TOF" << endl;
1849     cout << nprotons << "   primary protons   reached TOF" << endl;
1850     cout << nelectrons<<"   primary electrons reached TOF" << endl;
1851     cout << nmuons <<   "   primary muons     reached TOF" << endl;
1852     cout << "number of TOF hits for different species: 1-p, 2-K, 3-pi, 4-e, 5-mu, 6-all" << endl;
1853     cout << "   first number - track hits, second - noise ones, third - all" << endl;
1854     for (Int_t i=0;i<6;i++) cout << i+1 << "  " << particleInTOFtype[i][0] << "  " << particleInTOFtype[i][1] << "  " << particleInTOFtype[i][2] << endl; 
1855
1856     Int_t primaryReachedTOF[6];
1857     primaryReachedTOF[0]=npions;
1858     primaryReachedTOF[1]=nkaons;
1859     primaryReachedTOF[2]=nprotons;
1860     primaryReachedTOF[3]=nelectrons;
1861     primaryReachedTOF[4]=nmuons;
1862     primaryReachedTOF[5]=npions+nkaons+nprotons+nelectrons+nmuons;
1863     
1864     cout << " Reading TOF hits done" << endl;
1865   }
1866
1867   delete [] npixel;
1868 }
1869
1870 //____________________________________________________________________________
1871 void AliTOFReconstructioner::AddNoiseFromOuter(Option_t *option, Int_t ***MapPixels, AliTOFPad* pixelArray , AliTOFRecHit* hitArray, Int_t& isHitOnFiredPad, Int_t& ipixel)
1872 {
1873   //
1874   // Add noise hits from outer regions (forward and backward) according
1875   // to parameterized fZNoise distribution (to be used with events 
1876   // generated in the barrel region)
1877
1878   Float_t * zLen = new Float_t[AliTOFConstants::fgkNPlates+1];
1879   Float_t * zStrips = new Float_t[AliTOFConstants::fgkNPlates];
1880   zStrips[0]=(Float_t) (AliTOFConstants::fgkNStripC);
1881   zStrips[1]=(Float_t) (AliTOFConstants::fgkNStripB);
1882   zStrips[2]=(Float_t) (AliTOFConstants::fgkNStripA);
1883   zStrips[3]=(Float_t) (AliTOFConstants::fgkNStripB);
1884   zStrips[4]=(Float_t) (AliTOFConstants::fgkNStripC);
1885
1886   zLen[5]=AliTOFConstants::fgkzlenA*0.5+AliTOFConstants::fgkzlenB+AliTOFConstants::fgkzlenC;
1887   zLen[4]=zLen[5]-AliTOFConstants::fgkzlenC;
1888   zLen[3]=zLen[4]-AliTOFConstants::fgkzlenB;
1889   zLen[2]=zLen[3]-AliTOFConstants::fgkzlenA;
1890   zLen[1]=zLen[2]-AliTOFConstants::fgkzlenB;
1891   zLen[0]=zLen[1]-AliTOFConstants::fgkzlenC;
1892
1893   
1894   Int_t isector;    // random sector number 
1895   Int_t iplate;     // random plate number
1896   Int_t istrip;     // random strip number in the plate
1897   Int_t ipadAlongX; // random pad number along x direction
1898   Int_t ipadAlongZ; // random pad number along z direction
1899   Int_t ipad;
1900   Int_t nPixEdge=0; // additional pads fired due to edge effect when adding noise from outer
1901                     // regions
1902
1903   // x -> time of flight given in ns
1904   TF1 *noiseTof = new TF1("noiseTof","exp(-x/20)",0,100);
1905
1906   if(strstr(option,"pp")){
1907     fZnoise = new TF1("fZnoise","257.8-0.178*x-0.000457*x*x",-AliTOFConstants::fgkMaxhZtof,AliTOFConstants::fgkMaxhZtof);
1908   }
1909   if(strstr(option,"Pb-Pb")){
1910     fZnoise = new TF1("fZnoise","182.2-0.09179*x-0.0001931*x*x",-AliTOFConstants::fgkMaxhZtof,AliTOFConstants::fgkMaxhZtof);
1911   }
1912
1913   if(fNoise) {
1914     if(fdbg) cout << " Start adding additional noise  hits from outer regions" << endl;
1915
1916     for(Int_t i=0;i<fNoise;i++) {
1917
1918       isector=(Int_t) (AliTOFConstants::fgkNSectors*gRandom->Rndm())+1; //the sector number
1919       //  non-flat z-distribution of additional hits
1920       Float_t zNoise=fZnoise->GetRandom();
1921
1922       // holes for PHOS and HMPID
1923       if(((AliTOF *) gAlice->GetDetector("TOF"))->IsVersion()==2) {
1924         // to be checked the holes case
1925         if(isector>12 && isector<16) { // sectors 13,14,15 - RICH
1926           do {
1927             iplate=(Int_t) (AliTOFConstants::fgkNPlates*gRandom->Rndm())+1;
1928           } while (iplate==2 || iplate==3 || iplate==4);
1929           //         } else if(isector>11 && isector<17) { // sectors 12,13,14,15,16 - PHOS
1930         } else if(isector>2 && isector<8) { // sectors 3,4,5,6,7 - PHOS
1931           do {
1932             iplate=(Int_t) (AliTOFConstants::fgkNPlates*gRandom->Rndm())+1;
1933           } while (iplate==3);
1934         } else {
1935           iplate=(Int_t) (AliTOFConstants::fgkNPlates*gRandom->Rndm())+1;
1936         }
1937       } else {
1938         iplate=0;
1939         do {
1940           iplate++;
1941         } while(zNoise>zLen[iplate]);
1942       }
1943       // end of holes
1944
1945       if(iplate<1 || iplate>5) {
1946         printf("  iplate<1 or iplate>5, iplate=%i\n",iplate);
1947         return; 
1948       }
1949
1950       Float_t nStripes=0;
1951       if(iplate>1) {
1952         for (Int_t i=0;i<iplate-1;i++) {
1953           nStripes += zStrips[i];
1954         }
1955       }
1956
1957       istrip=(Int_t)(zNoise-zLen[iplate-1])/((zLen[iplate]-zLen[iplate-1])/zStrips[iplate-1]); //the strip number in the plate
1958       istrip++;
1959
1960       ipadAlongX = (Int_t)(AliTOFConstants::fgkNpadX*gRandom->Rndm())+1;
1961       ipadAlongZ = (Int_t)(AliTOFConstants::fgkNpadZ*gRandom->Rndm())+1;
1962       ipad=(Int_t)(ipadAlongZ-1)*AliTOFConstants::fgkNpadX+ipadAlongX;    //the pad number
1963       
1964       Float_t xStrip=(ipadAlongX-1)*AliTOFConstants::fgkXPad+AliTOFConstants::fgkXPad*gRandom->Rndm()-0.5*AliTOFConstants::fgkNpadX*AliTOFConstants::fgkXPad;//x-coor.in the strip frame
1965       Float_t zStrip=(ipadAlongZ-1)*AliTOFConstants::fgkZPad+AliTOFConstants::fgkZPad*gRandom->Rndm()-0.5*AliTOFConstants::fgkNpadZ*AliTOFConstants::fgkZPad;//z-coor.in the strip frame 
1966
1967       Int_t nPlace[4]={0,0,0,0};
1968       nPlace[0]=ipad;
1969
1970       Int_t   nActivatedPads=0;
1971       Int_t   nFiredPads=0;
1972       Bool_t  isFired[4]={kFALSE,kFALSE,kFALSE,kFALSE};
1973       Float_t tofAfterSimul[4]={0.,0.,0.,0.};
1974       Float_t qInduced[4]={0.,0.,0.,0.};
1975       Float_t averageTime=0.;    
1976       Float_t toffornoise=10.+noiseTof->GetRandom(); // 10 ns offset + parameterization [ns]
1977
1978       BorderEffect(zStrip,xStrip,toffornoise,nActivatedPads,nFiredPads,isFired,nPlace,qInduced,tofAfterSimul,averageTime); // simulate edge effect
1979
1980       if(nFiredPads) {
1981         for(Int_t indexOfPad=0; indexOfPad<nActivatedPads; indexOfPad++) {
1982           if(isFired[indexOfPad]){// the pad has fired
1983
1984             if(indexOfPad==0) {// the hit belongs to a fired pad
1985               isHitOnFiredPad++;
1986               hitArray[isHitOnFiredPad-1].SetX(0.);
1987               hitArray[isHitOnFiredPad-1].SetY(0.);
1988               hitArray[isHitOnFiredPad-1].SetZ(zNoise);
1989               hitArray[isHitOnFiredPad-1].SetNoise(1);
1990             } // close if(indexOfPad==0)
1991
1992             ipad = nPlace[indexOfPad];
1993             
1994             Int_t iMapValue=MapPixels[AliTOFConstants::fgkNSectors*(iplate-1)+isector-1][istrip-1][ipad-1];
1995             
1996             if(iMapValue==0) {
1997               ipixel++;
1998               if(indexOfPad) nPixEdge++;
1999               MapPixels[AliTOFConstants::fgkNSectors*(iplate-1)+isector-1][istrip-1][ipad-1]=ipixel;
2000               pixelArray[ipixel-1].SetGeom(isector,iplate,istrip,ipad);
2001               pixelArray[ipixel-1].AddState(1);
2002               pixelArray[ipixel-1].SetRealTime(tofAfterSimul[indexOfPad]);
2003               pixelArray[ipixel-1].SetHit(isHitOnFiredPad);
2004             } else if( tofAfterSimul[indexOfPad] < pixelArray[iMapValue-1].GetRealTime() ) {
2005               pixelArray[iMapValue-1].SetTrack(-1);
2006               pixelArray[iMapValue-1].AddState(1);
2007               pixelArray[iMapValue-1].SetRealTime(tofAfterSimul[indexOfPad]);
2008               pixelArray[iMapValue-1].SetHit(isHitOnFiredPad);
2009             }  //end of if(iMapValue==0)
2010             
2011           }// close if(isFired[indexOfPad])
2012         }  //end loop on activated pads indexOfPad
2013       } // close if(nFiredPads)
2014     }  //end of NOISE cycle
2015   }
2016
2017   // free used memory
2018   if (fZnoise)
2019     {
2020       delete fZnoise;
2021       fZnoise = 0;
2022     }
2023
2024   if (noiseTof)
2025     {
2026       delete noiseTof;
2027        noiseTof = 0;
2028     }
2029
2030   Int_t nNoiseSignals=0;
2031   Int_t nAll=0;
2032   for(Int_t idummy=1; idummy<ipixel+1; idummy++) {
2033     if(hitArray[pixelArray[idummy-1].GetHit()-1].GetNoise()==1) {
2034       nNoiseSignals++;
2035       if(pixelArray[idummy-1].GetState()>10) nAll++;
2036     }
2037   }
2038
2039   if(fdbg) {
2040     cout << " after adding " << fNoise << " noise hits: " << ipixel << " fired pixels (" << nNoiseSignals << " fired by noise, " << nAll << " noise+track)" << endl;
2041     printf(" %i additional pixels are fired by noise due to edge effect\n",nPixEdge);
2042     cout << " End of adding additional noise hits from outer regions" << endl;
2043   }
2044
2045   Float_t occupancy;
2046   // numberOfPads for AliTOFV4 (Full coverage) 
2047   // - to be upgraded checking the used TOF version -
2048   Float_t numberOfPads=AliTOFConstants::fgkPadXSector*AliTOFConstants::fgkNSectors;
2049   occupancy=100.*ipixel/numberOfPads;   // percentage of fired pads
2050   printf(" Overall TOF occupancy (percentage of fired pads after adding noise) = %f\n",occupancy); 
2051   delete [] zLen;
2052   delete [] zStrips;
2053   
2054 }
2055
2056
2057 //____________________________________________________________________________
2058 void AliTOFReconstructioner::SetMinDistance(AliTOFRecHit* hitArray, Int_t ilastEntry)
2059 {
2060   //
2061   // Set the distance to the nearest hit for hitArray
2062   // ilastEntry is the index of the last entry of hitArray
2063
2064   // starting the setting for the distance to the nearest TOFhit (cm)
2065   for(Int_t i=0; i<ilastEntry; i++) {
2066     
2067     if(hitArray[i].GetFirst()==1 && hitArray[i].GetNoise()==0) { // select the first hit of the track 
2068       // hits are not due to noise
2069       Float_t minDistance=10000.,squareDistance; // current values of the (square) distance
2070       Int_t jAtMin=0;                            // index of the hit nearest to the i-th hit
2071       Float_t xhit=hitArray[i].X(); // x coordinate for i-th hit
2072       Float_t yhit=hitArray[i].Y(); // y coordinate for i-th hit
2073       Float_t zhit=hitArray[i].Z(); // z coordinate for i-th hit
2074       //  was    for(Int_t j=0; j<isHitOnFiredPad; j++) {
2075       for(Int_t j=0; j<ilastEntry; j++) {
2076         if(i!=j) {
2077           squareDistance=(hitArray[j].X()-xhit)*(hitArray[j].X()-xhit)+
2078             (hitArray[j].Y()-yhit)*(hitArray[j].Y()-yhit)+
2079             (hitArray[j].Z()-zhit)*(hitArray[j].Z()-zhit);
2080           if(squareDistance<minDistance) {
2081             minDistance=squareDistance;
2082             jAtMin=j;
2083           }
2084         }
2085       }
2086       minDistance=TMath::Sqrt(minDistance);
2087       hitArray[i].SetRmin(minDistance);
2088       if(minDistance==0.) printf(" Rmin=0, i=%i, j=%i, x=%f,y=%f,z=%f\n",i,jAtMin,xhit,yhit,zhit);// it cannot happen
2089     }
2090   }
2091
2092 }
2093
2094 // these lines has to be commented till TPC will provide fPx fPy fPz 
2095 // and fL in AliTPChit class
2096 //____________________________________________________________________________ 
2097 /*
2098 void AliTOFReconstructioner::ReadTPCHits(Int_t ntracks, TTree* treehits, TClonesArray* tpchits, Int_t* iTrackPt, Int_t* iparticle, Float_t* ptTrack, AliTOFTrack* trackArray, Int_t& itrack)
2099 {
2100   //
2101   // Read TPC hits for the current event
2102   // 
2103   TParticle *particle=0;
2104   Int_t npions=0;    // number of pions for the current event
2105   Int_t nkaons=0;    // number of kaons for the current event
2106   Int_t nprotons=0;  // number of protons for the current event
2107   Int_t nelectrons=0;// number of electrons for the current event
2108   Int_t nmuons=0;    // number of muons for the current event
2109   Int_t ntotalTPChits=0; // total number of TPC hits for the current event
2110   Int_t idummy=-1;       // dummy var used to count double hit TPC cases
2111   Int_t nTpcDoubleHitsLastRow=0; // number of double TPC hits in the last pad row
2112   Int_t nTpcHitsLastRow=0;       // number of TPC hits in the last pad row
2113   Float_t trdpos[2]={0.,0.};
2114   Float_t pos[3];               // TPC hit position
2115   Float_t mom[3]; // momentum components in the last TPC row
2116   Float_t pt=0., tpclen; // pt: transverse momentum in the last TPC row
2117   Int_t nbytes = 0;
2118   Int_t ipart=0, nhits=0, iprim=0;
2119
2120   itrack=0; // itrack: total number of selected TPC tracks
2121
2122   // speed-up the code
2123   treehits->SetBranchStatus("*",0); // switch off all branches
2124   treehits->SetBranchStatus("TPC*",1); // switch on only TPC
2125
2126   for (Int_t track=0; track<ntracks;track++) {
2127     gAlice->ResetHits();
2128     nbytes += treehits->GetEvent(track);
2129     
2130     
2131     nhits = tpchits->GetEntriesFast();
2132     
2133     for (Int_t hit=0;hit<nhits;hit++) {
2134       ntotalTPChits++;
2135       AliTPChit* tpcHit = (AliTPChit*)tpchits->UncheckedAt(hit);
2136       Int_t row = tpcHit->fPadRow;
2137       ipart    = tpcHit->GetTrack();
2138       if(ipart>=fMaxAllTracks) break;
2139       particle = (TParticle*)gAlice->Particle(ipart);
2140       Int_t pdgCode=particle->GetPdgCode();
2141       // only high momentum tracks
2142       // momentum components at production vertex
2143       Float_t pxvtx = particle->Px();
2144       Float_t pyvtx = particle->Py();
2145       Float_t pzvtx = particle->Pz();
2146       Float_t pvtx = TMath::Sqrt(pxvtx*pxvtx+pyvtx*pyvtx+pzvtx*pzvtx);
2147       if(pvtx>fPBound && row == fLastTPCRow) {
2148         Float_t vx = particle->Vx();
2149         Float_t vy = particle->Vy();
2150         Float_t vr = TMath::Sqrt(vx*vx+vy*vy);
2151         Float_t x = tpcHit->X();
2152         Float_t y = tpcHit->Y();
2153         Float_t z = tpcHit->Z();
2154         pos[0]=x; pos[1]=y; pos[2]=z;
2155         
2156         Float_t pxtpc = tpcHit->fPx;
2157         Float_t pytpc = tpcHit->fPy;
2158         Float_t pztpc = tpcHit->fPz;
2159         mom[0]=pxtpc; mom[1]=pytpc; mom[2]=pztpc; 
2160         Float_t momtpc = TMath::Sqrt(pxtpc*pxtpc+pytpc*pytpc+pztpc*pztpc);
2161         
2162         if(x*pxtpc+y*pytpc>0) { // only tracks going out of TPC
2163           
2164           Float_t isoutgoing = x*pxtpc+y*pytpc+z*pztpc;
2165           isoutgoing /= (momtpc*TMath::Sqrt(x*x+y*y+z*z));
2166           tpclen = tpcHit->fL;
2167           
2168           
2169           if(ipart!=idummy) {
2170             if(particle->GetFirstMother() < 0) {
2171               Int_t abspdgCode=TMath::Abs(pdgCode);
2172               switch (abspdgCode) {
2173               case 211:
2174                 npions++;
2175                 break ;
2176               case 321:
2177                 nkaons++;
2178                 break ;
2179               case 2212:
2180                 nprotons++;
2181                 break ;
2182               case 11:
2183                 nelectrons++;
2184                 break ;
2185               case 13:
2186                 nmuons++;
2187                 break ;
2188               }
2189             } // close if(particle->GetFirstMother() < 0)
2190           } // close if(ipart!=idummy)
2191           
2192           if(gRandom->Rndm()<fTrackingEfficiency && vr<fRadiusvtxBound && ipart!=idummy) {
2193             
2194             itrack++;
2195             if(particle->GetFirstMother() < 0) iprim++;
2196             
2197             if(itrack>fMaxTracks) {
2198               cout << "itrack=" << itrack << " > MAXTRACKS=" << fMaxTracks << endl;
2199               return;
2200             } // close if(itrack>fMaxTracks)
2201             
2202             
2203             iparticle[ipart]=itrack;
2204             
2205             trackArray[itrack-1].SetTrack(ipart,pvtx,pdgCode,tpclen,pos,mom,trdpos);
2206             
2207             pt=TMath::Sqrt(pxtpc*pxtpc+pytpc*pytpc); // pt: transverse momentum at TPC
2208             // Filling iTrackPt[MAXTRACKS] by itrack ordering on Pt
2209             if(itrack==1) {
2210               iTrackPt[itrack-1]=itrack;
2211               ptTrack[itrack-1]=pt;
2212             } else {
2213               for (Int_t i=0; i<itrack-1; i++) {
2214                 if(pt>ptTrack[i]) {
2215                   for(Int_t j=i; j<itrack-1; j++) {
2216                     Int_t k=itrack-1+i-j;
2217                     iTrackPt[k]= iTrackPt[k-1];
2218                     ptTrack[k] = ptTrack[k-1];
2219                   }
2220                   iTrackPt[i]=itrack;
2221                   ptTrack[i]=pt;
2222                   break;
2223                 }
2224                 if(i==itrack-2) {
2225                   iTrackPt[itrack-1]=itrack;
2226                   ptTrack[itrack-1]=pt;
2227                 }
2228               }
2229             }
2230             
2231           }  //end of itrack
2232           if(vr>fRadiusvtxBound) nTpcHitsLastRow++;
2233           if(ipart==idummy) nTpcDoubleHitsLastRow++;
2234           idummy=ipart;
2235         }  // close if(x*px+y*py>0)
2236       }  // close if(pvtx>fPBound && row == fLastTPCRow)
2237     }  //end of hits  
2238   }  // close loop on tracks   
2239   
2240   
2241   if(fdbg) {
2242     cout << ntotalTPChits << " TPC hits in the last TPC row " << fLastTPCRow << endl;
2243     cout << "   " << nTpcHitsLastRow << " - hits with Rvtx>fRadiusvtxBound=" << fRadiusvtxBound << endl;
2244     cout << "   " << nTpcDoubleHitsLastRow << " double TPC hits" << endl;
2245     cout << itrack    << " - extracted TPC tracks   "     << iprim << " - primary" << endl;
2246     cout << npions    << " primary pions reached TPC"     << endl;
2247     cout << nkaons    << " primary kaons reached TPC"     << endl;
2248     cout << nprotons  << " primary protons reached TPC"   << endl;
2249     cout << nelectrons<< " primary electrons reached TPC" << endl;
2250     cout << nmuons    << " primary muons reached TPC"     << endl;
2251   } // if(fdbg)
2252   
2253   Int_t primaryInTPC[6]={0,0,0,0,0,0};
2254   primaryInTPC[0]=npions;
2255   primaryInTPC[1]=nkaons;
2256   primaryInTPC[2]=nprotons;
2257   primaryInTPC[3]=nelectrons;
2258   primaryInTPC[4]=nmuons;
2259   primaryInTPC[5]=npions+nkaons+nprotons+nelectrons+nmuons;
2260   
2261   if(fdbg) {
2262     printf("  contents of iTrackPt[MAXTRACKS],PtTrack[MAXTRACKS]\n");
2263     for (Int_t i=0; i<itrack; i++) {
2264       printf(" %i : iTrackPt=%i, PtTrack=%f\n",i+1,iTrackPt[i],ptTrack[i]); 
2265     }
2266     printf(" Check ordered transverse momentum array\n");
2267     for (Int_t i=itrack-1; i>=0; i--) {
2268       printf(" %i : iTrackPt=%i, PtTrack=%f\n",i+1,iTrackPt[i],ptTrack[i]); 
2269     }
2270   }// if(fdbg)
2271   
2272 }
2273 */
2274 //____________________________________________________________________________
2275 void cylcor(Float_t& x, Float_t& y) {
2276   Float_t rho,phi;
2277   
2278   rho=TMath::Sqrt(x*x+y*y);
2279   phi=0.;
2280   if(TMath::Abs(x)>0. || TMath::Abs(y)>0.) phi=TMath::ATan2(y,x);
2281   if(phi<0.) phi=phi+2.*TMath::Pi();
2282   x=rho;
2283   y=phi;
2284   
2285 }
2286
2287 //____________________________________________________________________________
2288 void AliTOFReconstructioner::Matching(AliTOFTrack* trackArray, AliTOFRecHit* hitArray, Int_t ***mapPixels, AliTOFPad* pixelArray, Int_t* kTOFhitFirst, Int_t& ipixel, Int_t* iTrackPt, Int_t* iTOFpixel, Int_t ntotTpcTracks)
2289 {
2290   Int_t TestTracks,iTestTrack,itest,wPixel=0,itestc;
2291   Int_t * ntest = new Int_t[fMaxTestTracks];
2292   Int_t * testPixel = new Int_t[fMaxTestTracks];
2293   Float_t wLength=0.,wRho=0.,wZ=0.;
2294   Float_t * testLength = new Float_t[fMaxTestTracks];
2295   Float_t * testRho = new Float_t[fMaxTestTracks];
2296   Float_t * testZ = new Float_t[fMaxTestTracks];
2297   Float_t weight;
2298   Float_t * testWeight = new Float_t[fMaxTestTracks];
2299   Float_t rotationFactor,phi0,coslam,sinlam,helixRadius,xHelixCenter,yHelixCenter,zHelixCenter,helixFactor;
2300   Int_t npixel[5],iMapValue,iwork1,iwork2,iwork3,iwork4,ihit=0;
2301   Int_t charge[48]={ 0, 1,-1, 0, 1,-1, 0, 1,-1, 0,
2302                      1,-1, 0, 1,-1, 0, 0, 0, 1, 0,
2303                      -1, 0,-1,-1, 0, 0,-1, 0, 1, 0,
2304                      1, 1, 0, 0, 1,-1, 0, 0, 1,-1,
2305                      1, 1,-1, 0, 1, 1, 2, 0};
2306   Float_t theta0,gpx,gpy,gpz,gp,gpt,gtheta,gx,gy,gz,gr,gxLast,gyLast,gzLast,chargeField;
2307   Float_t sumOfTheta=0.,weightTestTracksOutTof[4];
2308   Float_t s,ds,xRespectToHelixCenter,yRespectToHelixCenter,deltaRadius,fp,xp,yp,grho;
2309   Float_t mass,energy,g;
2310   Int_t itrack=0,itr,particleCharge,istep,iplate=0,iPadAlongX=0;  
2311   Int_t itra,t34=0,t32=0,t44=0,t43=0,t42=0;
2312   Int_t wstate=0,m2state=0,wPix;
2313   Int_t idelR=0,idelR1=0,idelR2=0,iRmin=0,iRmin1=0,iRmin2=0;
2314   Float_t massArray[50] = {0.0,0.00051,0.00051,0.0,0.1057,0.1057,0.135,0.1396,0.1396,0.4977,
2315                        0.4936,0.4936,0.9396,0.9383,0.9383,0.4977,0.5488,1.1156,1.1894,1.1926,1.1926,
2316                        1.3149,1.3213,1.6724,0.9396,1.1156,1.1894,1.1926,1.1974,1.3149,
2317                        0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
2318   Float_t delR;
2319   Float_t radius,area,normR,normS,cosAngl;
2320   Int_t iPlateFirst,iTestGmax=0;
2321   Int_t fstate,iPrintM1=0,iPrintM2=0;
2322   Float_t gxExtrap=0.,gyExtrap=0.,gzExtrap=0.;
2323   Float_t avSigZ=0,avSigRPHI=0,avSigP=0,avSigPHI=0,avSigTHETA=0;
2324
2325   Float_t gxW,gyW,gzW;
2326   Float_t length0;
2327   Float_t snr=0;
2328   Int_t indexOfTestTrack;
2329   Float_t zPad,xPad;
2330   Int_t istate=0,imax=0,match,iMaxTestTracksOutTof=0,matchw;
2331   Float_t w,wmax=0.,inverseOfParticleSpeed,w2,smat[9],largestWeightTracksOutTof,sw;
2332   Float_t sumWeightTracksOutTof,sGeomWeigth;
2333   Int_t imatched;
2334   Int_t m10=0,m20=0,m22=0,m23=0;
2335   Int_t PRINT=0;
2336   TParticle *particle;
2337
2338   Float_t time=0.;
2339   itr=ntotTpcTracks;
2340   printf(" itr=%i\n",itr);
2341   for (itra=1; itra<itr+1; itra++) {
2342
2343     Int_t itrack=iTrackPt[itra-1];
2344     if(itrack==0) printf("  iTrackPt[itra-1]=0 for itra=%i\n",itra);
2345     Int_t ipart=trackArray[itrack-1].GetTrack(); 
2346     Float_t pvtx=trackArray[itrack-1].GetP();
2347     Int_t pdgCode=trackArray[itrack-1].GetPdgCode();
2348     Float_t tpclength=trackArray[itrack-1].GetlTPC();
2349     Float_t x=trackArray[itrack-1].GetRxTPC();
2350     Float_t y=trackArray[itrack-1].GetRyTPC();
2351     Float_t z=trackArray[itrack-1].GetRzTPC();
2352     Float_t RxTPC=x;
2353     Float_t RyTPC=y;
2354     Float_t RzTPC=z;
2355     Float_t Wx=x;
2356     Float_t Wy=y;
2357     Float_t Wz=z;
2358     Float_t px=trackArray[itrack-1].GetPxTPC();
2359     Float_t py=trackArray[itrack-1].GetPyTPC();
2360     Float_t pz=trackArray[itrack-1].GetPzTPC();
2361     Float_t pxTPC=px;
2362     Float_t pyTPC=py;
2363     Float_t pzTPC=pz;
2364
2365     Float_t p = TMath::Sqrt(px*px+py*py+pz*pz);
2366     Float_t pTPC=p;
2367
2368     Float_t rho = TMath::Sqrt(x*x+y*y);
2369     Float_t phi=0.;
2370     if(TMath::Abs(x)>0. || TMath::Abs(y)>0.) phi=TMath::ATan2(y,x);
2371     if(phi<0.) phi=phi+2.*TMath::Pi();
2372     Float_t phiTPC=phi*kRaddeg;
2373     if(fSigmavsp) {
2374       if(p==0) printf(" p=%f in g=0.022/p\n",p);
2375       g=0.022/p;
2376       avSigRPHI += g;      // (cm)
2377       if(rho==0) printf(" rho=%f in phi += g*gRandom->Gaus()/rho\n",rho);
2378       phi += g*gRandom->Gaus()/rho; 
2379     } else {
2380       if(rho==0) printf(" rho=%f in phi += (SIGMARPHI*gRandom->Gaus()/rho\n",rho);
2381       phi += (fSigmarphi*gRandom->Gaus()/rho);
2382     }
2383     x=rho*TMath::Cos(phi);
2384     y=rho*TMath::Sin(phi);
2385     Float_t zTPC=z;
2386     if(fSigmavsp) {
2387       if(p==0) printf(" p=%f in g=0.0275/p\n",p);
2388       g=0.0275/p;
2389       avSigZ += g;      // (cm)
2390       z += g*gRandom->Gaus();
2391     } else {
2392       z += fSigmaZ*gRandom->Gaus();
2393     }
2394
2395     // smearing on TPC momentum
2396
2397     {                                                                             
2398       Float_t pmom,phi,theta,arg;
2399       
2400       pmom=TMath::Sqrt(px*px+py*py+pz*pz);
2401       phi=0.;
2402       if(TMath::Abs(px)>0. || TMath::Abs(py)>0.) phi=TMath::ATan2(py,px);
2403       if(phi<0.) phi=phi+2*TMath::Pi();
2404       arg=1.;
2405       if(pmom>0.) arg=pz/pmom;
2406       theta=0.;
2407       if(TMath::Abs(arg)<=1.) theta=TMath::ACos(arg);
2408       
2409       if(fSigmavsp) {
2410         if(pmom<=0) printf(" pmom=%f in g = TMath::Abs(TMath::Log(pmom)/TMath::Log(10)+0.5)/0.7\n",pmom);
2411         g = TMath::Abs(TMath::Log(pmom)/TMath::Log(10)+0.5)/0.7;
2412         g = 0.01*(g*g*g+1.5)*1.24;
2413         avSigP += g;
2414         pmom *= (1+g*gRandom->Gaus());
2415         
2416         if(p<10) {
2417           if(pmom<=0) printf(" pmom=%f in g = 1-TMath::Log(pmom)/TMath::Log(10)\n",pmom);
2418           g = 1-TMath::Log(pmom)/TMath::Log(10);
2419           g = 0.001*(g*g*g+0.3)*0.65;  // (radian)
2420         } else {
2421           g = 0.001*0.3*0.65;
2422         }
2423         avSigPHI += g;
2424         phi += g*gRandom->Gaus();
2425         avSigTHETA += g;
2426         theta += g*gRandom->Gaus();
2427         
2428       } else {
2429         pmom *= (1+fSigmap*gRandom->Gaus());
2430         phi += fSigmaPhi*gRandom->Gaus();
2431         theta += fSigmaTheta*gRandom->Gaus();
2432       }
2433       gxW=px;
2434       gyW=py;
2435       gzW=pz;
2436       
2437       px=pmom*TMath::Sin(theta)*TMath::Cos(phi);
2438       py=pmom*TMath::Sin(theta)*TMath::Sin(phi);
2439       pz=pmom*TMath::Cos(theta);
2440
2441       
2442       if(x*px+y*py<=0) {
2443         x=Wx;
2444         y=Wy;
2445         z=Wz;
2446         px=gxW;
2447         py=gyW;
2448         pz=gzW;
2449       }// if(x*px+y*py<=0)
2450     }
2451     
2452     p = TMath::Sqrt(px*px+py*py+pz*pz);
2453     
2454     particleCharge=charge[PDGtoGeantCode(pdgCode)-1];
2455     mass=massArray[PDGtoGeantCode(pdgCode)-1];
2456     mass=massArray[8-1];       //we take pion mass for all tracks
2457     //             mass=massArray[14-1];       //here we take proton mass for all tracks
2458     energy=TMath::Sqrt(p*p+mass*mass);
2459     chargeField=particleCharge*fField;
2460     
2461     g=fRadLenTPC/( (x*px+y*py)/(rho*p) );
2462     
2463     if(g<=0) printf(" error, g<=0: g=%f, itra=%i, x,y,px,py=%f, %f, %f, %f\n",g,itra,x,y,px,py);
2464     
2465     theta0=13.6*0.001*TMath::Sqrt(g)*(1.+0.038*TMath::Log(g))*energy/(p*p);
2466  
2467     
2468     // start Loop on test tracks
2469     sumOfTheta=0.;
2470     for(Int_t i=0;i<4;i++) {
2471       weightTestTracksOutTof[i]=0.;
2472     }
2473     
2474     itest=0;
2475     for(Int_t i=0;i<fMaxTestTracks;i++) {
2476       ntest[i]=0;
2477       testPixel[i]=0;
2478       testLength[i]=0.;
2479       testRho[i]=0.;
2480       testZ[i]=0.;
2481       testWeight[i]=0.;
2482     }
2483     
2484     iPlateFirst=0;
2485     TestTracks=0;
2486     iTestTrack=0;
2487     iTestGmax=0;
2488     
2489     length0=0;
2490     
2491     for (indexOfTestTrack=0; indexOfTestTrack<fMaxTestTracks; indexOfTestTrack++) {
2492
2493       iTestTrack++;
2494       gpx=px;
2495       gpy=py;
2496       gpz=pz;
2497       gp=p;
2498       if(indexOfTestTrack) {
2499         gtheta=theta0;
2500         EpMulScatt(gpx,gpy,gpz,gp,gtheta);
2501         
2502       } else {
2503         gtheta=0;
2504       }
2505       
2506       weight=TMath::Exp(-gtheta*gtheta/(2*theta0*theta0));
2507       sumOfTheta += gtheta;
2508       
2509       //    ==========================================================
2510       // Calculate crossing of the track in magnetic field with cylidrical surface
2511       // of radius RTOFINNER
2512       //   chargeField = qB, where q is a charge of a particle in units of e,
2513       //                     B is magnetic field in tesla
2514       //   see 3.3.1.1. in the book "Data analysis techniques for
2515       //   high-energy physics experiments", edited by M.Regler
2516       //   in Russian: "Metody analiza dannykh v fizicheskom eksperimente"
2517       //   Moskva, "Mir", 1993. ctr.306
2518       
2519       // Initial constants
2520       rotationFactor=1.;
2521       if(chargeField<0.) rotationFactor=-1.;
2522       rotationFactor=-rotationFactor;
2523       gpt=gpx;
2524       phi0=gpy;
2525       cylcor(gpt,phi0);
2526       phi0 -= rotationFactor*TMath::Pi()*0.5;
2527       //               phi0 -= h*PID2;
2528       coslam=gpt/gp;
2529       sinlam=gpz/gp;
2530       //      helixRadius=100.*gpt/TMath::Abs(0.299792458*chargeField);
2531       helixRadius=100.*gpt/TMath::Abs(AliTOFConstants::fgkSpeedOfLight*chargeField);
2532       xHelixCenter=x-helixRadius*TMath::Cos(phi0);
2533       yHelixCenter=y-helixRadius*TMath::Sin(phi0);
2534       zHelixCenter=z;
2535       helixFactor=rotationFactor*coslam/helixRadius;
2536       
2537       //   Solves the equation f(s)=r(s)-RTOFINNER=0 by the Newton's method:
2538       //   snew=s-f/f'
2539       istep=0;
2540       s=AliTOFConstants::fgkrmin-TMath::Sqrt(x*x+y*y);;
2541       do {
2542         istep++;
2543         xRespectToHelixCenter=helixRadius*TMath::Cos(phi0+s*helixFactor);
2544         yRespectToHelixCenter=helixRadius*TMath::Sin(phi0+s*helixFactor);
2545         gx=xHelixCenter+xRespectToHelixCenter;
2546         gy=yHelixCenter+yRespectToHelixCenter;
2547         gr=TMath::Sqrt(gx*gx+gy*gy);
2548         deltaRadius=gr-AliTOFConstants::fgkrmin;
2549         xp=-helixFactor*yRespectToHelixCenter;
2550         yp= helixFactor*xRespectToHelixCenter;
2551         fp=(gx*xp+gy*yp)/gr;
2552         ds=deltaRadius/fp;
2553         s -= ds;
2554         if(istep==20) {
2555           istep=0;
2556           break;
2557         }
2558       } while (TMath::Abs(ds)>0.01);
2559       
2560       
2561       if(istep==0) goto end;
2562       
2563       //   Steps along the circle till a pad
2564       wPixel=0;
2565       wLength=0.;
2566       iplate=0;
2567       iPadAlongX=0;
2568       grho=0.;
2569       ds=fStep;
2570       gxLast=xHelixCenter+helixRadius*TMath::Cos(phi0+s*helixFactor);
2571       gyLast=yHelixCenter+helixRadius*TMath::Sin(phi0+s*helixFactor);
2572       gzLast=zHelixCenter+s*sinlam;
2573
2574       
2575       do {
2576         istep++;
2577         s += ds;
2578         gx=xHelixCenter+helixRadius*TMath::Cos(phi0+s*helixFactor);
2579         gy=yHelixCenter+helixRadius*TMath::Sin(phi0+s*helixFactor);
2580         gz=zHelixCenter+s*sinlam;
2581         rho=TMath::Sqrt(gx*gx+gy*gy);
2582         
2583         IsInsideThePad(fg3,gx,gy,gz,npixel,zPad,xPad);
2584         
2585         iplate += npixel[1];
2586         iPadAlongX += npixel[4];
2587         
2588         if(indexOfTestTrack==0 && iplate && iPlateFirst==0) {
2589           iPlateFirst=1;
2590           length0=s;
2591
2592           radius=s*3*theta0;
2593           area=TMath::Pi()*radius*radius;
2594           normR=TMath::Sqrt(gx*gx+gy*gy);
2595           normS=TMath::Sqrt((gx-gxLast)*(gx-gxLast)+
2596                      (gy-gyLast)*(gy-gyLast)+
2597                      (gz-gzLast)*(gz-gzLast));
2598
2599           cosAngl=(gx*(gx-gxLast)+gy*(gy-gyLast))/(normR*normS);
2600           if(cosAngl<0) printf(" cosAngl<0: gx=%f,gy=%f,  gxLast=%f,gyLast=%f,gzLast=%f\n",gx,gy,gxLast,gyLast,gzLast);
2601
2602           area /= cosAngl;
2603           TestTracks=(Int_t) (2*area/(AliTOFConstants::fgkXPad * AliTOFConstants::fgkZPad));
2604
2605           if(TestTracks<12) TestTracks=12;
2606
2607           // Angles of entering into the TOF plate
2608
2609           Int_t iZ=0;
2610           if(TMath::Abs(gz)>300) {
2611             iZ=4;
2612           } else if(TMath::Abs(gz)>200) {
2613             iZ=3;
2614           } else if(TMath::Abs(gz)>100) {
2615             iZ=2;
2616           } else if(TMath::Abs(gz)>0) {
2617             iZ=1;
2618           }
2619           
2620           
2621         } // end of if(indexOfTestTrack==0 && iplate && iPlateFirst==0)
2622
2623         
2624         if(npixel[4]>0) {
2625
2626           iwork1=npixel[0];
2627           iwork2=npixel[1];
2628           iwork3=npixel[2];
2629           //               iwork4=npixel[3];
2630           iwork4=(npixel[3]-1)*AliTOFConstants::fgkNpadX+npixel[4];
2631
2632           Int_t ifirstindex=AliTOFConstants::fgkNSectors*(npixel[1]-1)+npixel[0];
2633           iMapValue=mapPixels[ifirstindex-1][iwork3-1][iwork4-1];
2634           if(iMapValue==0) {
2635             ipixel++;
2636             if(ipixel>fMaxPixels) {
2637               cout << "ipixel=" << ipixel << " > MAXPIXELS=" << fMaxPixels << endl;
2638               break;
2639             }
2640             mapPixels[ifirstindex-1][iwork3-1][iwork4-1]=ipixel;
2641             pixelArray[ipixel-1].SetGeom(iwork1,iwork2,iwork3,iwork4);
2642             iMapValue=ipixel;
2643           }
2644           
2645           wPixel=iMapValue;
2646           wLength=tpclength+s;
2647           wRho=rho;
2648           wZ=gz;
2649           
2650           ihit=kTOFhitFirst[ipart];
2651           
2652           if(ihit) {
2653             if(indexOfTestTrack==0) {
2654               {
2655                 idelR++;
2656                 delR=TMath::Sqrt((gx-hitArray[ihit-1].X())*(gx-hitArray[ihit-1].X())+
2657                           (gy-hitArray[ihit-1].Y())*(gy-hitArray[ihit-1].Y())+
2658                           (gz-hitArray[ihit-1].Z())*(gz-hitArray[ihit-1].Z()));
2659
2660               }
2661               
2662               if(delR>hitArray[ihit-1].GetRmin()) iRmin++;
2663               gxExtrap=gx;
2664               gyExtrap=gy;
2665               gzExtrap=gz;
2666             } else {
2667               delR=TMath::Sqrt((gx-gxExtrap)*(gx-gxExtrap)+
2668                         (gy-gyExtrap)*(gy-gyExtrap)+
2669                         (gz-gzExtrap)*(gz-gzExtrap));
2670             }
2671           }  //end of if(ihit)
2672           
2673           break;
2674           
2675         }  //end of npixel[4]
2676         
2677         if(rho<grho) {
2678           istep=0;
2679           break;
2680         }
2681         grho=rho;
2682         
2683         gxLast=gx;
2684         gyLast=gy;
2685         gzLast=gz;
2686         
2687       } while(rho<AliTOFConstants::fgkrmax); //end of do 
2688
2689       
2690       if(istep>0) {
2691         if(iplate) {
2692           if(iPadAlongX==0) {
2693             istep=-3;            // holes in TOF
2694           }
2695         } else {
2696           if(TMath::Abs(gz)<AliTOFConstants::fgkMaxhZtof) {
2697             //                   if(TMath::Abs(gz)<MAXZTOF2) {
2698             istep=-2;            // PHOS and RICH holes or holes in between TOF plates
2699           } else {
2700             istep=-1;            // out of TOF on z-size
2701           }
2702         }
2703       }
2704       
2705       if(iPadAlongX>0) {
2706         if(itest==0) {
2707           itest=1;
2708           ntest[itest-1]=1;
2709           testPixel[itest-1]=wPixel;
2710           testLength[itest-1]=wLength;
2711           testRho[itest-1]=wRho;
2712           testZ[itest-1]=wZ;
2713           testWeight[itest-1]=weight;
2714         } else {
2715           Int_t k;
2716           for(Int_t i=0;i<itest;i++) {
2717             k=0;
2718             if(testPixel[i]==wPixel) {
2719               k=1;
2720               ntest[i]++;
2721               testLength[i] += wLength;
2722               testRho[i] += wRho;
2723               testZ[i] += wZ;
2724               testWeight[i] += weight;
2725               break;
2726             }
2727           }  //end for i
2728           if(k==0) {
2729             itest++;
2730             ntest[itest-1]=1;
2731             testPixel[itest-1]=wPixel;
2732             testLength[itest-1]=wLength;
2733             testRho[itest-1]=wRho;
2734             testZ[itest-1]=wZ;
2735             testWeight[itest-1]=weight;
2736           }
2737         }
2738       }
2739       
2740     end: ;
2741       //   Statistics
2742       if(fMatchingStyle==1) {
2743         if(istep>-4 && istep<1) weightTestTracksOutTof[-istep] ++;
2744       } else {
2745         if(istep>-4 && istep<1) weightTestTracksOutTof[-istep] += weight;
2746       }
2747       
2748       if(fMatchingStyle==2) {
2749         if(indexOfTestTrack==0 && istep==0) break;
2750         if(indexOfTestTrack+1==TestTracks) break;
2751       }
2752       
2753     }  //end of indexOfTestTrack
2754
2755     snr += (Float_t) (indexOfTestTrack+1);
2756     
2757     //   Search for the "hole" with the largest weigth
2758     largestWeightTracksOutTof=0.;
2759     sumWeightTracksOutTof=0.;
2760     for(Int_t i=0;i<4;i++) {
2761       w=weightTestTracksOutTof[i];
2762       sumWeightTracksOutTof += w;
2763       if(w>largestWeightTracksOutTof) {
2764         largestWeightTracksOutTof=w;
2765         iMaxTestTracksOutTof=i;
2766       }
2767     }
2768     
2769     itestc=itest;
2770     if(itest>0) {
2771       for(Int_t i=0;i<itest;i++) {
2772         testLength[i] /= ntest[i];
2773         testRho[i] /= ntest[i];
2774         testZ[i] /= ntest[i];
2775       }
2776       //   Search for the pixel with the largest weigth
2777       wmax=0.;
2778       wstate=0;
2779       sw=0;
2780       sGeomWeigth=0;
2781       for(Int_t i=0;i<itest;i++) {
2782         istate=pixelArray[testPixel[i]-1].GetState();
2783         fstate=0;
2784         if(istate>0) {
2785           fstate=1;
2786           wstate++;
2787         }
2788         if(fMatchingStyle==1) {
2789           sGeomWeigth += ntest[i];
2790           w=(fpadefficiency*fstate+(1.-fpadefficiency)*(1-fstate))*ntest[i];
2791           if(pixelArray[testPixel[i]-1].GetTrackMatched()>0) w *= 0.1;
2792         } else {
2793           sGeomWeigth += testWeight[i];
2794           w=(fpadefficiency*fstate+(1.-fpadefficiency)*(1-fstate))*testWeight[i];
2795           if(pixelArray[testPixel[i]-1].GetTrackMatched()>0) w *= 0.1;
2796         }
2797         
2798         // weighting according to the Pulse Height (we use the square of weight)
2799         // if (fChargeFactorForMatching) w *= (pixelArray[testPixel[i]-1].GetCharge())*(pixelArray[testPixel[i]-1].GetCharge());
2800         if (fChargeFactorForMatching && fstate==1) w *= (pixelArray[testPixel[i]-1].GetCharge())*(pixelArray[testPixel[i]-1].GetCharge());
2801
2802         if(w>wmax) {
2803           wmax=w;
2804           imax=i;
2805         }
2806         sw += w;
2807       }
2808       wPixel=testPixel[imax];
2809       wLength=testLength[imax];
2810       istate=pixelArray[wPixel-1].GetState();
2811       
2812       //Choose the TOF dead space
2813       //               if(istate==0 && largestWeightTracksOutTof>wmax) {
2814       //               if(istate==0 && largestWeightTracksOutTof>=sw) {
2815       if(istate==0 && sumWeightTracksOutTof>sGeomWeigth) {
2816         itestc=itest;
2817         itest=0;
2818       }
2819     }
2820     
2821     if(itest>0) {
2822       
2823       //   Set for MyTrack: Pixel
2824       trackArray[itrack-1].SetPixel(wPixel);
2825       
2826       istate=pixelArray[wPixel-1].GetState();
2827       
2828       if(istate) {
2829         
2830         //   Set for MyTrack: Pixel, Length, TOF, MassTOF
2831         //fp
2832         //time=pixelArray[wPixel-1].GetTime();
2833         time=pixelArray[wPixel-1].GetRealTime();
2834         trackArray[itrack-1].SetLength(wLength);
2835         trackArray[itrack-1].SetTof(time);
2836         
2837         inverseOfParticleSpeed=time/wLength;
2838         //w=900.*inverseOfParticleSpeed*inverseOfParticleSpeed-1.;
2839         w=(100.*AliTOFConstants::fgkSpeedOfLight)*(100.*AliTOFConstants::fgkSpeedOfLight)*inverseOfParticleSpeed*inverseOfParticleSpeed-1.;
2840         w2=pvtx*pvtx;
2841         Float_t squareMass=w2*w;
2842         mass=TMath::Sqrt(TMath::Abs(squareMass));
2843         if(w<0.) mass=-mass;
2844         
2845         trackArray[itrack-1].SetMassTOF(mass);
2846         
2847         //   Set for MyTrack: Matching
2848         match=4;
2849         //                 if(ipart==pixelArray[wPixel-1].GetTrack()) match=3;
2850         if( (ipart==pixelArray[wPixel-1].GetTrack()) && hitArray[pixelArray[wPixel-1].GetHit()-1].GetNoise()==0)match=3;
2851         imatched=pixelArray[wPixel-1].GetTrackMatched();
2852         //   Set for TOFPixel the number of matched track
2853         pixelArray[wPixel-1].SetTrackMatched(itrack);
2854         
2855         if(imatched>0) {
2856           matchw=trackArray[imatched-1].GetMatching();
2857           if(match==3 && matchw==4) t34++;
2858           if(match==3 && matchw==2) t32++;
2859           if(match==4 && matchw==4) t44++;
2860           if(match==4 && matchw==3) t43++;
2861           if(match==4 && matchw==2) t42++;
2862           if(iTOFpixel[ipart]==0 || iTOFpixel[trackArray[imatched-1].GetTrack()]==0) {
2863             m20++;
2864           } else if(iTOFpixel[ipart]==iTOFpixel[trackArray[imatched-1].GetTrack()]) {
2865             m22++;
2866           } else {
2867             m23++;
2868             wPix=iTOFpixel[ipart];
2869             if(PRINT && iPrintM1==10 && iPrintM2<10) {
2870               if(iPrintM2==0) {
2871                 printf("*** test print for tracks matched with the pixel for with we had matched track\n");
2872               }
2873               iPrintM2++;
2874               printf(" m=2: ipart=%i, pdgCode=%i, p=%f, theta0=%f, %i Pixel(LP=%i,SP=%i,P=%i) \n", 
2875                      ipart,pdgCode,p,theta0,wPix,
2876                      pixelArray[wPix-1].GetSector(),pixelArray[wPix-1].GetPlate(),pixelArray[wPix-1].GetPixel());
2877               printf("      mat=%i, %i Pixel(LP=%i,SP=%i,P=%i), Test(n=%i,i=%i,w=%f,z=%f), wst=%i \n",
2878                      match,wPixel,
2879                      pixelArray[wPixel-1].GetSector(),pixelArray[wPixel-1].GetPlate(),pixelArray[wPixel-1].GetPixel(),
2880                      itest,imax,wmax,testZ[imax],wstate);
2881               Int_t fstat,istat;
2882               for(Int_t i=0;i<itest;i++) {
2883                 wPix=testPixel[i];
2884                 istat=pixelArray[wPix-1].GetState();
2885                 fstat=0;
2886                 if(istat>0) fstat=1;
2887                 w=(fpadefficiency*fstat+(1.-fpadefficiency)*(1-fstat))*ntest[i];
2888                 if(istat>0)
2889                   printf("                     %i: %i Pixel(LP=%i,SP=%i,P=%i), istat=%i, ntest=%i, w=%f\n",i+1,
2890                          wPix,pixelArray[wPix-1].GetSector(),pixelArray[wPix-1].GetPlate(),pixelArray[wPix-1].GetPixel(),
2891                          istat,ntest[i],w);
2892               }
2893               printf("      mat=%i, %i Pixel \n",matchw,trackArray[imatched-1].GetPad());
2894             }
2895           }
2896           if(wstate>1) m2state++;
2897           smat[matchw+4]--;
2898           match=2;
2899           trackArray[imatched-1].SetMatching(match);
2900           smat[match+4]++;
2901           
2902         }  // if(imatched>0)
2903         
2904       } else {  //else if(istate)
2905         
2906         match=1;
2907         if(iTOFpixel[ipart]==0) m10++;
2908         if(PRINT && iPrintM1<10) {
2909           Int_t wPix;
2910           wPix=iTOFpixel[ipart];
2911           if(wPix) {
2912             if(iPrintM1==0) {
2913               printf("*** test print for tracks fired a pixel but matched with non-fired pixel\n");
2914             }
2915             iPrintM1++;
2916             printf(" m=1: itra=%i,ipart=%i, pdgCode=%i, p=%f, theta0=%f, %i Pixel(LP=%i,SP=%i,P=%i) \n", 
2917                    itra,ipart,pdgCode,p,theta0,wPix,
2918                    pixelArray[wPix-1].GetSector(),pixelArray[wPix-1].GetPlate(),pixelArray[wPix-1].GetPixel());
2919             printf("      mat=%i, %i Pixel(LP=%i,SP=%i,P=%i), Test(n=%i,i=%i,w=%f,z=%f), wst=%i \n",
2920                    match,wPixel,
2921                    pixelArray[wPixel-1].GetSector(),pixelArray[wPixel-1].GetPlate(),pixelArray[wPixel-1].GetPixel(),
2922                    itest,imax,wmax,testZ[imax],wstate);
2923             
2924           }
2925         }  //end if(PRINT && iPrintM1<10)
2926         
2927       }  //end if(istate)
2928       
2929     } else {
2930       match=-1-iMaxTestTracksOutTof;
2931       
2932     }  //end itest
2933     
2934     trackArray[itrack-1].SetMatching(match);
2935     //             if(iTestGmax==1) hMTT->Fill(match);
2936     smat[match+4]++;
2937
2938     sumOfTheta /= iTestTrack;
2939     
2940     itest=itestc;
2941     
2942     //Test
2943     if(PRINT) {
2944       if(iTOFpixel[ipart] && match!=3) {
2945         particle = (TParticle*)gAlice->Particle(ipart);  //for V3.05
2946
2947         printf("          ipixel=%i (Sector=%i, Plate=%i, Strip=%i, Pixel=%i), fired by %i track\n",iTOFpixel[ipart],pixelArray[iTOFpixel[ipart]-1].GetSector(),pixelArray[iTOFpixel[ipart]-1].GetPlate(),pixelArray[iTOFpixel[ipart]-1].GetStrip(),pixelArray[iTOFpixel[ipart]-1].GetPixel(),pixelArray[iTOFpixel[ipart]-1].GetTrack());   
2948         printf("     indexOfTestTrack=%i itest=%i weightTestTracksOutTof[4]=%f weightTestTracksOutTof[2]=%f weightTestTracksOutTof[1]=%f weightTestTracksOutTof[0]=%f\n",indexOfTestTrack,itest,weightTestTracksOutTof[3],weightTestTracksOutTof[2],weightTestTracksOutTof[1],weightTestTracksOutTof[0]);
2949         if(itest) {
2950
2951           printf("     take ipixel=%i (Sector=%i, Plate=%i, Strip=%i, Pixel=%i), (fired by %i track), match=%i\n",wPixel,pixelArray[wPixel-1].GetSector(),pixelArray[wPixel-1].GetPlate(),pixelArray[wPixel-1].GetStrip(),pixelArray[wPixel-1].GetPixel(),pixelArray[wPixel-1].GetTrack(),match);   
2952         }
2953       }
2954     }
2955     if(PRINT && itra<10 ) {
2956
2957       if(itest) {
2958         cout << "      number of pixels with test tracks=" << itest << endl;
2959         for(Int_t i=0;i<itest;i++) {
2960           cout << "      " << i+1 << "  tr.=" << ntest[i] << "  w=" << testWeight[i] << "  pix.= " << testPixel[i] << " (" << 
2961             pixelArray[testPixel[i]-1].GetSector() << " " << " " << pixelArray[testPixel[i]-1].GetPlate() << " " << 
2962             pixelArray[testPixel[i]-1].GetPixel() << " )" << "  l= " << testLength[i] << " sig=" << 
2963             theta0*(testLength[i]-tpclength) << "  rho= " << testRho[i] << "  z= " << testZ[i] << endl;
2964         }
2965         cout << "      pixel=" << wPixel << "  state=" << istate << "  l=" << wLength << "  TOF=" << time << "  m=" << mass << "  match=" << match <<  endl;
2966         if(istate>0) cout << "      fired by track " << pixelArray[wPixel-1].GetTrack() << endl;
2967       }
2968     }
2969   }  //end of track
2970   
2971
2972   if(itr) {
2973     printf(" %f probe tracks per 1 real track\n",snr/itr);   
2974     itrack=itr;
2975   }
2976   
2977   
2978   cout << ipixel << " - total number of TOF pixels after matching" << endl;
2979   w=iRmin;
2980   if(idelR!=0) {
2981     w /= idelR;
2982     printf(" %i tracks with delR, %f of them have delR>Rmin \n",idelR,w);
2983   }
2984   w=iRmin1;
2985   if(idelR1!=0) {
2986     w /= idelR1;
2987     printf(" %i tracks with delR1 (|z|<175), %f of them have delR>Rmin \n",idelR1,w);
2988   }
2989   w=iRmin2;
2990   if(idelR2!=0) {
2991     w /= idelR2;
2992     printf(" %i tracks with delR2 (|z|>175), %f of them have delR>Rmin \n",idelR2,w);
2993   }
2994   
2995   cout << " ********************  End of matching **********" << endl;
2996   delete [] ntest;
2997   delete [] testPixel;
2998   delete [] testLength;
2999   delete [] testRho;
3000   delete [] testZ;
3001   delete [] testWeight;
3002 }
3003
3004 //____________________________________________________________________________
3005 void AliTOFReconstructioner::FillNtuple(Int_t ntracks, AliTOFTrack* trackArray, AliTOFRecHit* hitArray, AliTOFPad* pixelArray, Int_t* iTOFpixel, Int_t* iparticle, Float_t* toftime, Int_t& ipixelLastEntry, Int_t itrack){
3006   
3007   // itrack : total number of TPC selected tracks
3008   // for the caller is ntotTPCtracks
3009   
3010   cout << " ********************  Start of searching non-matched fired pixels **********" << endl;
3011   const Int_t charge[48]={ 0, 1,-1, 0, 1,-1, 0, 1,-1, 0,
3012                            1,-1, 0, 1,-1, 0, 0, 0, 1, 0,
3013                            -1, 0,-1,-1, 0, 0,-1, 0, 1, 0,
3014                            1, 1, 0, 0, 1,-1, 0, 0, 1,-1,
3015                            1, 1,-1, 0, 1, 1, 2, 0};
3016
3017   Int_t macthm1=0;
3018   Int_t macthm2=0;
3019   Int_t macthm3=0;
3020   Int_t macthm4=0;
3021   Int_t macth0=0;
3022   Int_t macth1=0;
3023   Int_t macth2=0;
3024   Int_t macth3=0;
3025   Int_t macth4=0;
3026   
3027   
3028   Float_t smat[9],smat0[9],smat1[9];
3029   for(Int_t i=0;i<9;i++) {
3030     smat[i]=0.;
3031     smat0[i]=0.;
3032     smat1[i]=0.;
3033   }
3034   
3035   Int_t nFiredPixelsNotMatchedWithTracks=0;
3036   Int_t istate;
3037   for (Int_t i=0; i<ipixelLastEntry; i++) {
3038     istate=pixelArray[i].GetState();
3039     if(istate==0) break;
3040     if(pixelArray[i].GetTrackMatched()==-1) nFiredPixelsNotMatchedWithTracks++;
3041   }
3042   printf("  %i fired pixels have not matched tracks\n",nFiredPixelsNotMatchedWithTracks);
3043   cout << " ********************  End of searching non-matched fired pixels **********" << endl;
3044   
3045   Int_t nTPCHitMissing=0;
3046   for(Int_t i=0; i<ipixelLastEntry; i++) {
3047     if(pixelArray[i].GetHit()>0) {
3048       if(hitArray[pixelArray[i].GetHit()-1].GetNoise()==0) {
3049         if(iparticle[pixelArray[i].GetTrack()]==0) nTPCHitMissing++;
3050       }
3051     }
3052   }
3053   printf("  %i pixels fired by track hit without a hit on the last layer of TPC\n",nTPCHitMissing);
3054   
3055   
3056   Int_t icharge=0;   // total number of charged particles
3057   Int_t iprim=0;     // number of primaries
3058   Int_t ipions=0;    // number of primary pions
3059   Int_t ikaons=0;    // number of primary kaons
3060   Int_t iprotons=0;  // number of primary protons
3061   Int_t ielectrons=0;// number of primary electrons
3062   Int_t imuons=0;    // number of primary muons
3063   Float_t particleTypeArray[6][5][2];
3064   
3065   for (Int_t index1=0;index1<6;index1++) {
3066     for (Int_t index2=0;index2<5;index2++) {
3067       for (Int_t index3=0;index3<2;index3++) {
3068         particleTypeArray[index1][index2][index3]=0.;
3069       }
3070     }
3071   }
3072   
3073   Int_t nTOFhitsWithNoTPCTracks=0; // to be moved later when used
3074   
3075   /*
3076   TObjArray *Particles = gAlice->Particles();
3077   Int_t numberOfParticles=Particles->GetEntries();
3078   cout << "numberOfParticles " << numberOfParticles << endl;
3079   // fpdbg
3080   if(numberOfParticles>fMaxAllTracks) numberOfParticles=fMaxAllTracks;
3081   */
3082
3083   for (Int_t i=0; i<ntracks; i++) { // starting loop on all primaries charged particles for current event)
3084
3085     /*
3086     cout << "particle " << i << endl;
3087     cout << "total " << numberOfParticles << endl;
3088     */
3089     TParticle *part = (TParticle *) gAlice->Particle(i);
3090     if(charge[PDGtoGeantCode(part->GetPdgCode())-1]) {
3091       icharge++;
3092       /*
3093       cout << "charged particles " << icharge << endl;
3094       */
3095       Int_t particleType=0;
3096       Int_t absPdgCode = TMath::Abs(part->GetPdgCode());
3097       switch (absPdgCode) {
3098       case 211:
3099         particleType=3;
3100         break ;
3101       case 321:
3102         particleType=2;
3103         break ;
3104       case 2212:
3105         particleType=1;
3106         break ;
3107       case 11:
3108         particleType=4;
3109         break ;
3110       case 13:
3111         particleType=5;
3112         break ;
3113       }
3114       
3115       if(part->GetFirstMother() < 0) {
3116         iprim++;
3117         switch (particleType) {
3118         case 1:
3119           iprotons++;
3120           break ;
3121         case 2:
3122           ikaons++;
3123           break ;
3124         case 3:
3125           ipions++;
3126           break ;
3127         case 4:
3128           ielectrons++;
3129           break ;
3130         case 5:
3131           imuons++;
3132           break ;
3133         }
3134       }
3135       
3136       Int_t match=0;
3137       Float_t wLength=-1.;
3138       Float_t time=-1.;
3139       Float_t mass=-1.;
3140       
3141       Int_t itr=iparticle[i]; // get the track number for the current charged particle
3142       
3143       if(iTOFpixel[i]>0 && itr==0) nTOFhitsWithNoTPCTracks++;
3144       
3145       if(itr) {
3146         match=trackArray[itr-1].GetMatching();
3147         //cout << "match " << match << endl;
3148         wLength=trackArray[itr-1].GetLength();
3149         //cout << "wLength " << wLength << endl;
3150         time=trackArray[itr-1].GetTof();
3151         mass=trackArray[itr-1].GetMassTOF();
3152         //cout << "mext " << mass << endl;
3153         //        if(PRINT && (i>789 && i<800) ) cout << i << " track:  l=" << wLength << "  TOF=" << time << "  m=" << mass << "  match=" << match <<  endl; 
3154         if(iTOFpixel[i]==0) {
3155           smat0[match+4]++;
3156           wLength=-wLength;
3157         }
3158       }
3159       Int_t ikparen=part->GetFirstMother();
3160       Int_t imam;
3161       if(ikparen<0) {
3162         imam=0;
3163       } else {
3164         imam=part->GetPdgCode();
3165       }
3166       
3167       Int_t evnumber=gAlice->GetEvNumber();
3168       if(match==-1) macthm1++;
3169       if(match==-2) macthm2++;
3170       if(match==-3) macthm3++;
3171       if(match==-4) macthm4++;
3172       if(match==0) macth0++;
3173       if(match==1) macth1++;
3174       if(match==2) macth2++;
3175       if(match==3) macth3++;
3176       if(match==4) macth4++;
3177       foutputntuple->Fill(evnumber,part->GetPdgCode(),imam,part->Vx(),part->Vy(),part->Vz(),part->Px(),part->Py(),part->Pz(),toftime[i],wLength,match,time,mass);
3178       
3179       
3180       
3181       // -----------------------------------------------------------
3182       // Filling 2 dimensional Histograms true time vs matched time
3183       // Filling 1 dimensional Histogram true time - matched time
3184       //
3185       // time              = time associated to the matched pad [ns]
3186       //                     it could be the average time of the cluster fired
3187       //
3188       // toftime[i]        = real time (including pulse height delays) [s]
3189       //
3190       //
3191       // if (time>=0) {
3192       // if (imam==0) TimeTrueMatched->Fill(time, toftime[i]*1E+09);
3193       // if (imam==0) DeltaTrueTimeMatched->Fill(time-toftime[i]*1E+09);
3194       // }
3195       //
3196       //---------------------------------------------------------------
3197       
3198       if(match==-4 || match>0) {
3199         Int_t matchW;
3200         matchW=match;
3201         if(match==-4) matchW=1;
3202         if(particleType) {
3203           particleTypeArray[particleType-1][matchW-1][1]++;
3204           particleTypeArray[5][matchW-1][1]++;
3205           particleTypeArray[particleType-1][4][1]++;
3206           particleTypeArray[5][4][1]++;
3207           if(part->GetFirstMother() < 0) {
3208             particleTypeArray[particleType-1][matchW-1][0]++;
3209             particleTypeArray[5][matchW-1][0]++;
3210             particleTypeArray[particleType-1][4][0]++;
3211             particleTypeArray[5][4][0]++;
3212             
3213             // fill histos for QA
3214             //if(particleType==3 && matchW==3) hPiWithTrueTime->Fill(sqrt((part->Px())*(part->Px())+(part->Py())*(part->Py())+(part->Pz())*(part->Pz())));
3215             //if(particleType==2 && matchW==3) hKWithTrueTime->Fill(sqrt((part->Px())*(part->Px())+(part->Py())*(part->Py())+(part->Pz())*(part->Pz())));
3216             //if(particleType==1 && matchW==3) hPWithTrueTime->Fill(sqrt((part->Px())*(part->Px())+(part->Py())*(part->Py())+(part->Pz())*(part->Pz())));
3217             //
3218             
3219           } // close if(part->GetFirstMother() < 0)
3220         } // close if(particleType)
3221       } // close if(match==-4 || match>0)
3222     } // close if(charge[PDGtoGeantCode(part->GetPdgCode())-1])
3223   } // close for (Int_t i=0; i<ntracks; i++) {
3224
3225   cout <<  " macthm1 " << macthm1 << endl;
3226   cout <<  " macthm2 " << macthm2 << endl;
3227   cout <<  " macthm3 " << macthm3 << endl;
3228   cout <<  " macthm4 " << macthm4 << endl;
3229   cout <<  " macth0 " << macth0 << endl;
3230   cout <<  " macth1 " << macth1 << endl;
3231   cout <<  " macth2 " << macth2 << endl;
3232   cout <<  " macth3 " << macth3 << endl;
3233   cout <<  " macth4 " << macth4 << endl;
3234   
3235
3236   printf(" %i TOF hits have not TPC track\n",nTOFhitsWithNoTPCTracks);
3237   Int_t imatch=0;
3238   for(Int_t i=0;i<9;i++) {
3239     if(itrack) cout << "   " << smat[i]*100./itrack << " % of them (="<<smat[i]<<") have match=" << i-4 << "  " << smat0[i] << " have not TOF hits" << endl;
3240     if(i==0 || i>4) imatch += (Int_t) (smat[i]);
3241     
3242     //     cout << "   " << smat[i]*100./itrack << " % of them (="<<smat[i]<<") have match=" << i-4 << "  " << smat0[i] << " have not TOF hits" << "  " << smat1[i] << " have (r.p)<0 for first hit" << endl;
3243   }
3244   
3245   if(fdbg){
3246     /*
3247     cout << " nparticles = " << numberOfParticles << "  charged = " << icharge << "  prim.=" << iprim << endl;
3248     */
3249     cout << " nparticles = " << ntracks << "  charged = " << icharge << "  prim.=" << iprim << endl;
3250     cout << ipions << " - primary pions" << endl;
3251     cout << ikaons << " - primary kaons" << endl;
3252     cout << iprotons << " - primary protons" << endl;
3253     cout << ielectrons << " - primary electrons" << endl;
3254     cout << imuons << " - primary muons reached TPC" << endl;
3255     cout << " ********** " << imatch << " TPC tracks are matched with TOF pixels (incl.match=-4) **********" << endl;
3256   }
3257   
3258   /*
3259     Float_t PrimaryInBarrel[6],Acceptance[6];
3260     PrimaryInBarrel[0]=ipions;
3261     PrimaryInBarrel[1]=ikaons;
3262     PrimaryInBarrel[2]=iprotons;
3263     PrimaryInBarrel[3]=ielectrons;
3264     PrimaryInBarrel[4]=imuons;
3265     PrimaryInBarrel[5]=ipions+ikaons+iprotons+ielectrons+imuons;
3266     
3267     //   cout << "   TPC acceptance for the primary species: 1-p, 2-K, 3-pi, 4-e, 5-mu, 6-all" << endl; 
3268     for(Int_t i=0; i<6; i++) {
3269      Acceptance[i]=0.;
3270      if(PrimaryInBarrel[i]) Acceptance[i]=100.*PrimaryReachedTPC[i]/PrimaryInBarrel[i];
3271      //hTPCacceptance[i]->Fill(Acceptance[i]);
3272      //     printf(" species: %i    %f\n",i+1,Acceptance[i]);     
3273      }
3274      
3275      //   cout << "   TOF acceptance for the primary species: 1-p, 2-K, 3-pi, 4-e, 5-mu, 6-all" << endl; 
3276      for(Int_t i=0; i<6; i++) {
3277      Acceptance[i]=0.;
3278      if(PrimaryInBarrel[i]) Acceptance[i]=100.*PrimaryReachedTOF[i]/PrimaryInBarrel[i];
3279      //hTOFacceptance[i]->Fill(Acceptance[i]);
3280      //     printf(" species: %i    %f\n",i+1,Acceptance[i]);     
3281      }
3282      
3283    for (Int_t index1=0;index1<6;index1++) {
3284    for (Int_t index2=0;index2<4;index2++) {
3285    for (Int_t index3=0;index3<2;index3++) {
3286    if(particleTypeArray[index1][4][index3]) particleTypeArray[index1][index2][index3]=
3287                                                     100.*particleTypeArray[index1][index2][index3]/particleTypeArray[index1][4][index3]; 
3288                                                     }
3289      }
3290      }
3291      
3292    cout << "species: 1-p, 2-K, 3-pi, 4-e, 5-mu, 6-all" << endl; 
3293    cout << " matched pixels(%): 1-unfired 2-double 3-true 4-wrong 5-total number of tracks" << endl;
3294    
3295    cout << "  primary tracks:" << endl; 
3296    for (Int_t i=0;i<6;i++) {
3297      cout << i+1 << "  " << particleTypeArray[i][0][0] << "  " << particleTypeArray[i][1][0] << "  " << particleTypeArray[i][2][0] << "  " << particleTypeArray[i][3][0] << "  " << particleTypeArray[i][4][0] << endl; 
3298      }
3299      
3300      //   cout<<"      contam.for all prim.(%)="<<100*particleTypeArray[5][3][0]/(particleTypeArray[5][3][0]+particleTypeArray[5][2][0])<<endl;
3301
3302      cout << "  all tracks:" << endl; 
3303      for (Int_t i=0;i<6;i++) {
3304      cout << i+1 << "  " << particleTypeArray[i][0][1] << "  " << particleTypeArray[i][1][1] << "  " << particleTypeArray[i][2][1] << "  " << particleTypeArray[i][3][1] << "  " << particleTypeArray[i][4][1] << endl; 
3305    } 
3306    
3307    //   cout<<"      contam.for all (%)="<<100*particleTypeArray[5][3][1]/(particleTypeArray[5][3][1]+particleTypeArray[5][2][1])<<endl;
3308    //  printf(" t34=%i, t32=%i, t44=%i, t43=%i, t42=%i\n",t34,t32,t44,t43,t42);
3309    //  printf(" m10=%f, m20=%f, m22=%f, m23=%f, m2state=%i\n",m10,m20,m22,m23,m2state);
3310   */
3311 }