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