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