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