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