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