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