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