/************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ /* $Id$ */ //_________________________________________________________________________ // Manager class for TOF reconstruction. // // //-- Authors: Bologna-ITEP-Salerno Group // // Description: Manager class for TOF reconstruction (derived from TTask) // Summary of the main methods: // - extraction of the TPC (assumed to be) reconstructed tracks // comment: it has to me moved as soon as possible into a separate // class AliTOFTrackReader (K. Safarik suggestion) // - geometrical propagation of the above tracks till TOF detector // - matching of the tracks with the TOF signals // // Remark: the GEANT3.21 geometry is used during the geometrical propagation // of the tracks in order to know the current volume reached by the track. // ////////////////////////////////////////////////////////////////////////////// #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "AliConst.h" #include "AliDetector.h" #include "AliHeader.h" #include "AliLoader.h" #include "AliRun.h" #include "AliRun.h" #include "AliRunLoader.h" #include "AliTOF.h" #include "AliTOFConstants.h" #include "AliTOFHitMap.h" #include "AliTOFPad.h" #include "AliTOFRecHit.h" #include "AliTOFReconstructioner.h" #include "AliTOFSDigit.h" #include "AliTOFTrack.h" #include "AliTOFhit.h" #include "AliTOFv1.h" #include "AliTOFv2.h" #include "AliTOFv2FHoles.h" #include "AliTOFv3.h" #include "AliTOFv4.h" #include "AliTOFv4T0.h" #include "AliMC.h" // #include "../TPC/AliTPC.h" // AliTPChit class or somewhere // this line has to be commented till TPC will provide fPx fPy fPz and fL in ClassImp(AliTOFReconstructioner) //____________________________________________________________________________ AliTOFReconstructioner::AliTOFReconstructioner():TTask("AliTOFReconstructioner","") { // default ctor fNevents = 0 ; foutputfile = 0; foutputntuple= 0; fZnoise = 0; ftail = 0; } //____________________________________________________________________________ AliTOFReconstructioner::AliTOFReconstructioner(char* headerFile, Option_t* opt, char *RecFile ):TTask("AliTOFReconstructioner","") { // // ctor // fNevents = 0 ; // Number of events to reconstruct, 0 means all evens in current file foutputfile = 0; foutputntuple= 0; fZnoise = 0; ftail = 0; Init(opt); // create output file if (RecFile){ foutputfile= new TFile(RecFile,"RECREATE","root file for matching"); } else { char outFileName[100]; strcpy(outFileName,"match"); strcat(outFileName,headerFile); foutputfile= new TFile(outFileName,"RECREATE","root file for matching"); } // initialize the ALIROOT geometry gAlice->Init(); gAlice->Print(); CreateNTuple(); // add Task to //root/Tasks folder TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; roottasks->Add(this) ; } //____________________________________________________________________________ void AliTOFReconstructioner::Init(Option_t* opt) { // Initialize the AliTOFReconstructioner setting parameters for // reconstruction. // Option values: Pb-Pb for Pb-Pb events // pp for pp events // set common parameters fdbg=1; fNevents = 1; fFirstEvent = 1; fLastEvent = 1; fTimeResolution =0.120; fpadefficiency =0.99 ; fEdgeEffect = 2 ; fEdgeTails = 0 ; fHparameter = 0.4 ; fH2parameter = 0.15; fKparameter = 0.5 ; fK2parameter = 0.35; fEffCenter = fpadefficiency; fEffBoundary = 0.65; fEff2Boundary = 0.90; fEff3Boundary = 0.08; fResCenter = 50. ; fResBoundary = 70. ; fResSlope = 40. ; fTimeWalkCenter = 0. ; fTimeWalkBoundary=0. ; fTimeWalkSlope = 0. ; fTimeDelayFlag = 1 ; fPulseHeightSlope=2.0 ; fTimeDelaySlope =0.060; // was fMinimumCharge = TMath::Exp(fPulseHeightSlope*fKparameter/2.); fMinimumCharge = TMath::Exp(-fPulseHeightSlope*fHparameter); fChargeSmearing=0.0 ; fLogChargeSmearing=0.13; fTimeSmearing =0.022; fAverageTimeFlag=0 ; fChargeFactorForMatching=1; fTrackingEfficiency=1.0; // 100% TPC tracking efficiency assumed fSigmavsp = 1. ; fSigmaZ = 0. ; fSigmarphi= 0. ; fSigmap = 0. ; fSigmaPhi = 0. ; fSigmaTheta=0. ; fField = 0.2 ; // fRadLenTPC : 0.2 includes TRD / 0.03 TPC only fRadLenTPC=0.06 ; // last value fCorrectionTRD=0. ; fLastTPCRow=111 ; fRadiusvtxBound=50. ; // expressed in [cm] fStep = 0.1 ; // expressed in [cm] step during propagation of the // track inside TOF volumes fMatchingStyle=2 ; /* previous values default fMaxPixels=70000 ; fMaxAllTracks=70000 ; fMaxTracks=15000 ; */ fMaxPixels=165000 ; fMaxAllTracks=500000 ; fMaxTracks=15000 ; fMaxTOFHits=35000 ; fPBound =0.0 ; // bending effect: P_t=0.3*z*B*R , z particle charge fNoiseSlope=20. ; // set parameters as specified in opt //pp case if(strstr(opt,"pp")){ fMaxTestTracks=500 ; fNoise = 26. ; fNoiseMeanTof= 26.4 ; // to check } //Pb-Pb case if(strstr(opt,"Pb-Pb")){ fMaxTestTracks=20 ; fNoise = 9400. ; fNoiseMeanTof= 26.4 ; } } //____________________________________________________________________________ AliTOFReconstructioner::~AliTOFReconstructioner() { // // dtor // if (foutputfile) { delete foutputfile; foutputfile = 0; } if (foutputntuple) { delete foutputntuple; foutputntuple = 0; } if (fZnoise) { delete fZnoise; fZnoise = 0; } if (ftail) { delete ftail; ftail = 0; } } //____________________________________________________________________________ void AliTOFReconstructioner::CreateNTuple() { // // Create a Ntuple where information about reconstructed charged particles // (both primaries and secondaries) are stored // Variables: event ipart imam xvtx yvtx zvtx pxvtx pyvtx pzvtx time leng matc text mext // Meaning: // event - event number (0, 1, ...) // ipart - PDG code of particles // imam - PDG code for the parent // =0 for primary particle // xvtx - x-coordinate of the vertex (cm) // yvtx - y-coordinate of the vertex (cm) // zvtx - z-coordinate of the vertex (cm) // pxvtx - x-coordinate of the momentum in the vertex (GeV) // pyvtx - y-coordinate of the momentum in the vertex (GeV) // pzvtx - z-coordinate of the momentum in the vertex (GeV) // time - time of flight from TOF for given track (ps) - TOF time for the // first TOF hit of the track // leng - track length to the TOF pixel (cm), evaluate as a sum of the // track length from the track vertex to TPC and the average // length of the extrapolated track from TPC to TOF. // for the track without TOF hits leng=-abs(leng) // matc - index of the (TPC track) - (TOF pixel) matching // =0 for tracks which are not tracks for matching, i.e. // there is not hit on the TPC or Rvxt>200 cm // >0 for tracks with positive matching procedure: // =1 or 2 for non-identified tracks: // =1, if the corresponding pixel is not fired, // =2, if the corresponding pixel is also matched to the // other track, // =3 or 4 for identified tracks: // =3, if identified with true time, // =4, if identified with wrong time. // <0 for tracks with negative mathing procedure: // =-1, if track do not reach the pixel plate (curved in the // magnetic field), // =-2, if track is out of z-size of the TOF, // =-3, if track is or into the RICH hole, or into the PHOS hole, or in the space between the plates, // =-4, if track is into the dead space of the TOF. // text - time of fligth from the matching procedure = time of the // pixel corresponding to the track (ps) // =0 for the tracks with matc<=1 // mext - mass of the track from the matching procedure // =p*sqrt(900*(text/leng)**2-1), if 900*(text/leng)**2-1>=0 // =-p*sqrt(abs(900*(text/leng)**2-1)), if 900*(text/leng)**2-1<0 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 } //__________________________________________________________________ Double_t TimeWithTailR(Double_t* x, Double_t* par) { // sigma - par[0], alpha - par[1], part - par[2] // at xpart*sigma - TMath::Exp(-x/alpha) Float_t xx =x[0]; Double_t f; if(xxStart("TOFReconstruction"); AliRunLoader *rl = AliRunLoader::Open(datafile); if (rl == 0x0) { Error("Exec","Can not open session for file %s",datafile); return; } // Get AliRun object from file or create it if not on file rl->LoadgAlice(); gAlice = rl->GetAliRun(); AliTOF* TOF = (AliTOF *) gAlice->GetDetector ("TOF"); AliDetector* TPC = gAlice->GetDetector("TPC"); if (!TOF) { Error("AliTOFReconstructioner","TOF not found"); delete rl; return; } if (!TPC) { Error("AliTOFReconstructioner","TPC Detector not found"); delete rl; return; } AliLoader* tpcloader = rl->GetLoader("TPCLoader"); if (tpcloader == 0x0) { Error("AliTOFReconstructioner","Can not get TPC Loader from Run Loader."); delete rl; return; } AliLoader* tofloader = rl->GetLoader("TOFLoader"); if (tofloader == 0x0) { Error("AliTOFReconstructioner","Can not get TOF Loader from Run Loader."); delete rl; return; } if (fEdgeTails) ftail = new TF1("tail",TimeWithTailR,-2,2,3); if (fNevents == 0) fNevents = rl->GetNumberOfEvents(); // You have to set the number of event with the ad hoc setter // see testrecon.C if (rl->GetHeader() == 0x0) rl->LoadHeader(); tofloader->LoadHits(); tpcloader->LoadHits(); for (Int_t ievent = 0; ievent < fNevents; ievent++) { // start loop on events rl->GetEvent(ievent); Int_t nparticles= rl->GetHeader()->GetNtrack(); if (nparticles <= 0) return; TClonesArray* tofhits=0; TClonesArray* tpchits=0; if (TOF) tofhits = TOF->Hits(); if (TPC) tpchits = TPC->Hits(); TTree *TH = tofloader->TreeH(); if (!TH) return; Int_t ntracks = (Int_t) (TH->GetEntries()); // primary tracks cout << "number of primary tracked tracks in current event " << ntracks << endl; // number of primary tracked tracks // array declaration and initialization // TOF arrays // Int_t mapPixels[AliTOFConstants::fgkNSectors*AliTOFConstants::fgkNPlates][AliTOFConstants::fgkNStripC][AliTOFConstants::fgkNpadZ*AliTOFConstants::fgkNpadX]; Int_t *** mapPixels = new Int_t**[AliTOFConstants::fgkNSectors*AliTOFConstants::fgkNPlates]; for (Int_t i=0; i total number of fired pads (at least one time) // TPC arrays AliTOFTrack* trackArray = new AliTOFTrack[fMaxTracks]; Int_t * iparticle = new Int_t[fMaxAllTracks]; InitArray(iparticle,fMaxAllTracks); Int_t * iTrackPt = new Int_t[fMaxTracks]; InitArray(iTrackPt, fMaxTracks); // array Float_t * ptTrack = new Float_t[fMaxTracks]; InitArray( ptTrack, fMaxTracks); // array for selected track pt Int_t ntotTPCtracks=0; // total number of selected TPC tracks // reading TOF hits if(TOF) ReadTOFHits(ntracks, TH, tofhits, mapPixels, kTOFhitFirst, pixelArray, iTOFpixel, toftime, hitArray,isHitOnFiredPad,ntotFiredPads); cout << "isHitOnFiredPad " << isHitOnFiredPad << " for event " << ievent << endl; // start debug for adding noise // adding noise Int_t nHitsNoNoise=isHitOnFiredPad; if(fNoise) AddNoiseFromOuter(option,mapPixels,pixelArray,hitArray,isHitOnFiredPad,ntotFiredPads); cout << "ntotFiredPads after adding noise " << ntotFiredPads << " for event " << ievent << endl; // set the hitArray distance to nearest hit SetMinDistance(hitArray,nHitsNoNoise); // these lines has to be commented till TPC will provide fPx fPy fPz // and fL in AliTPChit class // reading TPC hits /* if(TPC) ReadTPCHits(ntracks, TH, tpchits, iTrackPt, iparticle, ptTrack, trackArray,ntotTPCtracks); */ // geometrical matching if(TOF && TPC) Matching(trackArray,hitArray,mapPixels,pixelArray,kTOFhitFirst,ntotFiredPads,iTrackPt,iTOFpixel,ntotTPCtracks); // fill ntuple with reconstructed particles from current event FillNtuple(ntracks,trackArray,hitArray,pixelArray,iTOFpixel,iparticle,toftime,ntotFiredPads,ntotTPCtracks); // free used memory delete [] toftime; delete [] pixelArray; delete [] iTOFpixel; delete [] kTOFhitFirst; delete [] hitArray; delete [] trackArray; delete [] iparticle; delete [] iTrackPt; delete [] ptTrack; for (Int_t i=0; icd(); //foutputntuple->Write(0,TObject::kOverwrite); foutputntuple->Write(); foutputfile->Write(); foutputfile->Close(); gBenchmark->Stop("TOFReconstruction"); cout << "AliTOFReconstructioner:" << endl ; cout << " took " << gBenchmark->GetCpuTime("TOFReconstruction") << " seconds in order to make the reconstruction for " << fNevents << " events " << endl; cout << gBenchmark->GetCpuTime("TOFReconstruction")/fNevents << " seconds per event " << endl ; cout << endl ; } //__________________________________________________________________ void AliTOFReconstructioner::SetRecFile(char * file ) { // // Set the file name for reconstruction output // if(!fRecFile.IsNull()) cout << "Changing destination file for TOF reconstruction from " <<(char *)fRecFile.Data() << " to " << file << endl ; fRecFile=file ; } //__________________________________________________________________ void AliTOFReconstructioner::Print(Option_t* /*option*/)const { // // Print reconstruction output file name // cout << "------------------- "<< GetName() << " -------------" << endl ; if(fRecFile.IsNull()) cout << " Writing reconstructed particles to file galice.root "<< endl ; else cout << " Writing reconstructed particle to file " << (char*) fRecFile.Data() << endl ; } //__________________________________________________________________ void AliTOFReconstructioner::PrintParameters()const { // // Print parameters used for reconstruction // cout << " ------------------- "<< GetName() << " -------------" << endl ; cout << " Parameters used for TOF reconstruction " << endl ; // Printing the parameters cout << " Number of events: " << fNevents << endl; cout << " Recostruction from event "<< fFirstEvent << " to event "<< fLastEvent << endl; cout << " TOF geometry parameters " << endl; cout << " Min. radius of the TOF (cm) "<< AliTOFConstants::fgkrmin << endl; cout << " Max. radius of the TOF (cm) "<< AliTOFConstants::fgkrmax << endl; cout << " Number of TOF geom. levels "<< AliTOFConstants::fgkmaxtoftree<< endl; cout << " Number of TOF sectors "<< AliTOFConstants::fgkNSectors << endl; cout << " Number of TOF modules "<< AliTOFConstants::fgkNPlates << endl; cout << " Max. Number of strips in a module "<< AliTOFConstants::fgkNStripC << endl; cout << " Number of pads per strip "<< AliTOFConstants::fgkNpadX*AliTOFConstants::fgkNpadZ << endl; cout << " Number of strips in central module "<< AliTOFConstants::fgkNStripA << endl; cout << " Number of strips in intermediate modules "<< AliTOFConstants::fgkNStripB << endl; cout << " Number of strips in outer modules "<< AliTOFConstants::fgkNStripC << endl; cout << " Number of MRPC in x strip direction "<< AliTOFConstants::fgkNpadX<< endl; cout << " Size of MRPC (cm) along X "<< AliTOFConstants::fgkXPad<< endl; cout << " Number of MRPC in z strip direction "<< AliTOFConstants::fgkNpadZ<Gmedia(xTOF, numed); gcvolu=g3->Gcvolu(); nLevel=gcvolu->nlevel; if(fdbg) { for (Int_t i=0; inames[i]),4); cout<<"Vol. n."<number[i]<=2) { // sector type name: B071(1,2,...,10),B074(1,2,3,4,5-PHOS),B075(1,2,3-RICH) strncpy(name,(char*) (&gcvolu->names[2]),4); // volume copy: 1,2,...,10 for B071, 1,2,3,4,5 for B074, 1,2,3 for B075 copyNumber=gcvolu->number[2]; if(!strcmp(name,"B071")) { if (copyNumber>=6 && copyNumber<=8) { sector=copyNumber+10; } else if (copyNumber>=1 && copyNumber<=5){ sector=copyNumber+7; } else { sector=copyNumber-8; } } else if(!strcmp(name,"B075")) { sector=copyNumber+12; } else if(!strcmp(name,"B074")) { if (copyNumber>=1 && copyNumber<=3){ sector=copyNumber+4; } else { sector=copyNumber-1; } } } if(sector) { nGeom[0]=sector; if(nLevel>=4) { // we'll use the module value in z-direction: // 1 2 3 4 5 // the module order in z-direction: FTOC,FTOB,FTOA,FTOB,FTOC // the module copy: 2 2 0 1 1 // module type name: FTOA, FTOB, FTOC strncpy(name,(char*) (&gcvolu->names[4]),4); // module copy: copyNumber=gcvolu->number[4]; if(!strcmp(name,"FTOC")) { if (copyNumber==2) { module=1; } else { module=5; } } else if(!strcmp(name,"FTOB")) { if (copyNumber==2) { module=2; } else { module=4; } } else if(!strcmp(name,"FTOA")) { module=3; } } } if(module) { nGeom[1]=module; if(nLevel>=6) { // strip type name: FSTR strncpy(name,(char*) (&gcvolu->names[6]),4); // strip copy: copyNumber=gcvolu->number[6]; if(!strcmp(name,"FSTR")) strip=copyNumber; } } if(strip) { nGeom[2]=strip; if(nLevel>=8) { // padz type name: FSEZ strncpy(name,(char*) (&gcvolu->names[8]),4); // padz copy: copyNumber=gcvolu->number[8]; if(!strcmp(name,"FSEZ")) padz=copyNumber; } } if(padz) { nGeom[3]=padz; if(nLevel>=9) { // padx type name: FSEX strncpy(name,(char*) (&gcvolu->names[9]),4); // padx copy: copyNumber=gcvolu->number[9]; if(!strcmp(name,"FSEX")) padx=copyNumber; } } if(padx) { nGeom[4]=padx; zPad=gcvolu->glx[2]; // check here xPad=gcvolu->glx[0]; // check here } // printf(" nGeom[0,1,2,3,4]=%i,%i,%i,%i,%i\n",nGeom[0],nGeom[1],nGeom[2],nGeom[3],nGeom[4]); } //__________________________________________________________________ void AliTOFReconstructioner::EpMulScatt(Float_t& px, Float_t& py, Float_t& pz, Float_t& p, Float_t& theta) { // Momentum p - before mult.scat. // Momentum p2 - after mult.scat. // THE0 - r.m.s. of deviation angle in plane // (see RPP'96: Phys.Rev.D54 (1996) 134) Float_t pt,thex,they,tantx,tanty,p2px,p2py,p2pz,costhe,sinthe,cospsi,sinpsi,p2x,p2y,p2z,p2,g; pt=TMath::Sqrt(px*px+py*py); // angles for p in the ' frame with Z'along p if(fMatchingStyle==1) { thex=theta*gRandom->Gaus(); they=theta*gRandom->Gaus(); } else { thex=3*(-theta+2*theta*gRandom->Rndm()); they=3*(-theta+2*theta*gRandom->Rndm()); } tantx=TMath::Tan(thex); tanty=TMath::Tan(they); // p2p - p2 in the ' frame p2pz=p/TMath::Sqrt(1.+tantx*tantx+tanty*tanty); p2py=p2pz*tanty; p2px=p2pz*tantx; // choose X'so that PHI=0 (see Il'in, Pozdnyak Analiticheskaya geometriya, 1968, c.88 // for Euler angles PSI, THETA (PHI=0) costhe=pz/p; sinthe=pt/p; cospsi=-py/pt; sinpsi=px/pt; // g=p2py*costhe-p2pz*sinthe; p2x=p2px*cospsi-g*sinpsi; p2y=p2px*sinpsi+g*cospsi; p2z=p2py*sinthe+p2pz*costhe; p2=TMath::Sqrt(p2x*p2x+p2y*p2y+p2z*p2z); // Test angle g=(px*p2x+py*p2y+pz*p2z)/(p*p2); if(g>1) g=1; theta=TMath::ACos(g); px=p2x; py=p2y; pz=p2z; p=p2; } // std border effect algorithm //__________________________________________________________________ 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) { // Input: z0, x0 - hit position in the strip system (0,0 - center of the strip), cm // geantTime - time generated by Geant, ns // Output: nActivatedPads - the number of pads activated by the hit (1 || 2 || 4) // nFiredPads - the number of pads fired (really activated) by the hit (nFiredPads <= nActivatedPads) // qInduced[iPad]- charge induced on pad, arb. units // this array is initialized at zero by the caller // tofAfterSimul[iPad] - time calculated with edge effect algorithm, ns // this array is initialized at zero by the caller // 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. // The weight is given by the qInduced[iPad]/qCenterPad // this variable is initialized at zero by the caller // nPlace[iPad] - the number of the pad place, iPad = 0, 1, 2, 3 // this variable is initialized at zero by the caller // // Description of used variables: // eff[iPad] - efficiency of the pad // res[iPad] - resolution of the pad, ns // timeWalk[iPad] - time walk of the pad, ns // timeDelay[iPad] - time delay for neighbouring pad to hited pad, ns // PadId[iPad] - Pad Identifier // E | F --> PadId[iPad] = 5 | 6 // A | B --> PadId[iPad] = 1 | 2 // C | D --> PadId[iPad] = 3 | 4 // nTail[iPad] - the tail number, = 1 for tailA, = 2 for tailB // qCenterPad - charge extimated for each pad, arb. units // weightsSum - sum of weights extimated for each pad fired, arb. units const Float_t kSigmaForTail[2] = {AliTOFConstants::fgkSigmaForTail1,AliTOFConstants::fgkSigmaForTail2}; //for tail Int_t iz = 0, ix = 0; Float_t dX = 0., dZ = 0., x = 0., z = 0.; Float_t h = fHparameter, h2 = fH2parameter, k = fKparameter, k2 = fK2parameter; Float_t effX = 0., effZ = 0., resX = 0., resZ = 0., timeWalkX = 0., timeWalkZ = 0.; Float_t logOfqInd = 0.; Float_t weightsSum = 0.; Int_t nTail[4] = {0,0,0,0}; Int_t padId[4] = {0,0,0,0}; Float_t eff[4] = {0.,0.,0.,0.}; Float_t res[4] = {0.,0.,0.,0.}; // Float_t qCenterPad = fMinimumCharge * fMinimumCharge; Float_t qCenterPad = 1.; Float_t timeWalk[4] = {0.,0.,0.,0.}; Float_t timeDelay[4] = {0.,0.,0.,0.}; nActivatedPads = 0; nFiredPads = 0; (z0 <= 0) ? iz = 0 : iz = 1; dZ = z0 + (0.5 * AliTOFConstants::fgkNpadZ - iz - 0.5) * AliTOFConstants::fgkZPad; // hit position in the pad frame, (0,0) - center of the pad z = 0.5 * AliTOFConstants::fgkZPad - TMath::Abs(dZ); // variable for eff., res. and timeWalk. functions iz++; // z row: 1, ..., AliTOFConstants::fgkNpadZ = 2 ix = (Int_t)((x0 + 0.5 * AliTOFConstants::fgkNpadX * AliTOFConstants::fgkXPad) / AliTOFConstants::fgkXPad); dX = x0 + (0.5 * AliTOFConstants::fgkNpadX - ix - 0.5) * AliTOFConstants::fgkXPad; // hit position in the pad frame, (0,0) - center of the pad x = 0.5 * AliTOFConstants::fgkXPad - TMath::Abs(dX); // variable for eff., res. and timeWalk. functions; ix++; // x row: 1, ..., AliTOFConstants::fgkNpadX = 48 ////// Pad A: nActivatedPads++; nPlace[nActivatedPads-1] = (iz - 1) * AliTOFConstants::fgkNpadX + ix; qInduced[nActivatedPads-1] = qCenterPad; padId[nActivatedPads-1] = 1; if (fEdgeEffect == 0) { eff[nActivatedPads-1] = fEffCenter; if (gRandom->Rndm() < eff[nActivatedPads-1]) { nFiredPads = 1; res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + fResCenter * fResCenter); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns; isFired[nActivatedPads-1] = kTRUE; tofTime[nActivatedPads-1] = gRandom->Gaus(geantTime + fTimeWalkCenter, res[0]); averageTime = tofTime[nActivatedPads-1]; } } else { if(z < h) { if(z < h2) { effZ = fEffBoundary + (fEff2Boundary - fEffBoundary) * z / h2; } else { effZ = fEff2Boundary + (fEffCenter - fEff2Boundary) * (z - h2) / (h - h2); } resZ = fResBoundary + (fResCenter - fResBoundary) * z / h; timeWalkZ = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * z / h; nTail[nActivatedPads-1] = 1; } else { effZ = fEffCenter; resZ = fResCenter; timeWalkZ = fTimeWalkCenter; } if(x < h) { if(x < h2) { effX = fEffBoundary + (fEff2Boundary - fEffBoundary) * x / h2; } else { effX = fEff2Boundary + (fEffCenter - fEff2Boundary) * (x - h2) / (h - h2); } resX = fResBoundary + (fResCenter - fResBoundary) * x / h; timeWalkX = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * x / h; nTail[nActivatedPads-1] = 1; } else { effX = fEffCenter; resX = fResCenter; timeWalkX = fTimeWalkCenter; } (effZ 0) { if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) { nActivatedPads++; nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX; eff[nActivatedPads-1] = effZ; res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ; // ns nTail[nActivatedPads-1] = 2; if (fTimeDelayFlag) { // qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.); // qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.); qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * z); logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing); timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing); } else { timeDelay[nActivatedPads-1] = 0.; } padId[nActivatedPads-1] = 2; } } ////// Pad C, D, E, F: if(x < k2) { effX = fEffBoundary - (fEffBoundary - fEff3Boundary) * (x / k2); } else { effX = fEff3Boundary * (k - x) / (k - k2); } resX = fResBoundary + fResSlope*x/k; timeWalkX = fTimeWalkBoundary + fTimeWalkSlope*x/k; if(x < k && x > 0) { // C: if(ix > 1 && dX < 0) { nActivatedPads++; nPlace[nActivatedPads-1] = nPlace[0] - 1; eff[nActivatedPads-1] = effX; res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns nTail[nActivatedPads-1] = 2; if (fTimeDelayFlag) { // qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.); // qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.); qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x); logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing); timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing); } else { timeDelay[nActivatedPads-1] = 0.; } padId[nActivatedPads-1] = 3; // D: if(z < k && z > 0) { if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) { nActivatedPads++; nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX - 1; eff[nActivatedPads-1] = effX * effZ; (resZGaus(-fPulseHeightSlope * z, fLogChargeSmearing); } else { // qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.); // qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.); qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x); logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing); } timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing); } else { timeDelay[nActivatedPads-1] = 0.; } padId[nActivatedPads-1] = 4; } } // end D } // end C // E: if(ix < AliTOFConstants::fgkNpadX && dX > 0) { nActivatedPads++; nPlace[nActivatedPads-1] = nPlace[0] + 1; eff[nActivatedPads-1] = effX; res[nActivatedPads-1] = 0.001 * (TMath::Sqrt(10400 + resX * resX)); // ns timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns nTail[nActivatedPads-1] = 2; if (fTimeDelayFlag) { // qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.); // qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.); qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x); logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing); timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing); } else { timeDelay[nActivatedPads-1] = 0.; } padId[nActivatedPads-1] = 5; // F: if(z < k && z > 0) { if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) { nActivatedPads++; nPlace[nActivatedPads - 1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX + 1; eff[nActivatedPads - 1] = effX * effZ; (resZGaus(-fPulseHeightSlope * z, fLogChargeSmearing); } else { // qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.); // qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.); qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x); logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing); } timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing); } else { timeDelay[nActivatedPads-1] = 0.; } padId[nActivatedPads-1] = 6; } } // end F } // end E } // end if(x < k) for (Int_t iPad = 0; iPad < nActivatedPads; iPad++) { if (res[iPad] < fTimeResolution) res[iPad] = fTimeResolution; if(gRandom->Rndm() < eff[iPad]) { isFired[iPad] = kTRUE; nFiredPads++; if(fEdgeTails) { if(nTail[iPad] == 0) { tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]); } else { ftail->SetParameters(res[iPad], 2. * res[iPad], kSigmaForTail[nTail[iPad]-1]); Double_t timeAB = ftail->GetRandom(); tofTime[iPad] = geantTime + timeWalk[iPad] + timeDelay[iPad] + timeAB; } } else { tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]); } if (fAverageTimeFlag) { averageTime += tofTime[iPad] * qInduced[iPad]; weightsSum += qInduced[iPad]; } else { averageTime += tofTime[iPad]; weightsSum += 1.; } } } if (weightsSum!=0) averageTime /= weightsSum; } // end else (fEdgeEffect != 0) } /* new algorithm (to be checked) //__________________________________________________________________ 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) { // Input: z0, x0 - hit position in the strip system (0,0 - center of the strip), cm // geantTime - time generated by Geant, ns // Output: nActivatedPads - the number of pads activated by the hit (1 || 2 || 4) // nFiredPads - the number of pads fired (really activated) by the hit (nFiredPads <= nActivatedPads) // qInduced[iPad]- charge induced on pad, arb. units // this array is initialized at zero by the caller // tofAfterSimul[iPad] - time calculated with edge effect algorithm, ns // this array is initialized at zero by the caller // 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. // The weight is given by the qInduced[iPad]/qCenterPad // this variable is initialized at zero by the caller // nPlace[iPad] - the number of the pad place, iPad = 0, 1, 2, 3 // this variable is initialized at zero by the caller // // Description of used variables: // eff[iPad] - efficiency of the pad // res[iPad] - resolution of the pad, ns // timeWalk[iPad] - time walk of the pad, ns // timeDelay[iPad] - time delay for neighbouring pad to hited pad, ns // PadId[iPad] - Pad Identifier // E | F --> PadId[iPad] = 5 | 6 // A | B --> PadId[iPad] = 1 | 2 // C | D --> PadId[iPad] = 3 | 4 // nTail[iPad] - the tail number, = 1 for tailA, = 2 for tailB // qCenterPad - charge extimated for each pad, arb. units // weightsSum - sum of weights extimated for each pad fired, arb. units const Float_t kSigmaForTail[2] = {AliTOFConstants::fgkSigmaForTail1,AliTOFConstants::fgkSigmaForTail2}; //for tail Int_t iz = 0, ix = 0; Float_t dX = 0., dZ = 0., x = 0., z = 0.; Float_t h = fHparameter, h2 = fH2parameter, k = fKparameter, k2 = fK2parameter; Float_t effX = 0., effZ = 0., resX = 0., resZ = 0., timeWalkX = 0., timeWalkZ = 0.; Float_t logOfqInd = 0.; Float_t weightsSum = 0.; Int_t nTail[4] = {0,0,0,0}; Int_t padId[4] = {0,0,0,0}; Float_t eff[4] = {0.,0.,0.,0.}; Float_t res[4] = {0.,0.,0.,0.}; Float_t qCenterPad = fMinimumCharge * fMinimumCharge; Float_t timeWalk[4] = {0.,0.,0.,0.}; Float_t timeDelay[4] = {0.,0.,0.,0.}; nActivatedPads = 0; nFiredPads = 0; (z0 <= 0) ? iz = 0 : iz = 1; dZ = z0 + (0.5 * AliTOFConstants::fgkNpadZ - iz - 0.5) * AliTOFConstants::fgkZPad; // hit position in the pad frame, (0,0) - center of the pad z = 0.5 * AliTOFConstants::fgkZPad - TMath::Abs(dZ); // variable for eff., res. and timeWalk. functions iz++; // z row: 1, ..., AliTOFConstants::fgkNpadZ = 2 ix = (Int_t)((x0 + 0.5 * AliTOFConstants::fgkNpadX * AliTOFConstants::fgkXPad) / AliTOFConstants::fgkXPad); dX = x0 + (0.5 * AliTOFConstants::fgkNpadX - ix - 0.5) * AliTOFConstants::fgkXPad; // hit position in the pad frame, (0,0) - center of the pad x = 0.5 * AliTOFConstants::fgkXPad - TMath::Abs(dX); // variable for eff., res. and timeWalk. functions; ix++; // x row: 1, ..., AliTOFConstants::fgkNpadX = 48 ////// Pad A: nActivatedPads++; nPlace[nActivatedPads-1] = (iz - 1) * AliTOFConstants::fgkNpadX + ix; qInduced[nActivatedPads-1] = qCenterPad; padId[nActivatedPads-1] = 1; if (fEdgeEffect == 0) { eff[nActivatedPads-1] = fEffCenter; if (gRandom->Rndm() < eff[nActivatedPads-1]) { nFiredPads = 1; res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + fResCenter * fResCenter); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns; isFired[nActivatedPads-1] = kTRUE; tofTime[nActivatedPads-1] = gRandom->Gaus(geantTime + fTimeWalkCenter, res[0]); averageTime = tofTime[nActivatedPads-1]; } } else { if(z < h) { if(z < h2) { effZ = fEffBoundary + (fEff2Boundary - fEffBoundary) * z / h2; } else { effZ = fEff2Boundary + (fEffCenter - fEff2Boundary) * (z - h2) / (h - h2); } resZ = fResBoundary + (fResCenter - fResBoundary) * z / h; timeWalkZ = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * z / h; nTail[nActivatedPads-1] = 1; } else { effZ = fEffCenter; resZ = fResCenter; timeWalkZ = fTimeWalkCenter; } if(x < h) { if(x < h2) { effX = fEffBoundary + (fEff2Boundary - fEffBoundary) * x / h2; } else { effX = fEff2Boundary + (fEffCenter - fEff2Boundary) * (x - h2) / (h - h2); } resX = fResBoundary + (fResCenter - fResBoundary) * x / h; timeWalkX = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * x / h; nTail[nActivatedPads-1] = 1; } else { effX = fEffCenter; resX = fResCenter; timeWalkX = fTimeWalkCenter; } (effZ 0) { if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) { nActivatedPads++; nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX; eff[nActivatedPads-1] = effZ; res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ; // ns nTail[nActivatedPads-1] = 2; if (fTimeDelayFlag) { qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.); qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.); logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing); timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing); } else { timeDelay[nActivatedPads-1] = 0.; } padId[nActivatedPads-1] = 2; } } ////// Pad C, D, E, F: if(x < k2) { effX = fEffBoundary - (fEffBoundary - fEff3Boundary) * (x / k2); } else { effX = fEff3Boundary * (k - x) / (k - k2); } resX = fResBoundary + fResSlope*x/k; timeWalkX = fTimeWalkBoundary + fTimeWalkSlope*x/k; if(x < k && x > 0) { // C: if(ix > 1 && dX < 0) { nActivatedPads++; nPlace[nActivatedPads-1] = nPlace[0] - 1; eff[nActivatedPads-1] = effX; res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns nTail[nActivatedPads-1] = 2; if (fTimeDelayFlag) { qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.); qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.); logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing); timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing); } else { timeDelay[nActivatedPads-1] = 0.; } padId[nActivatedPads-1] = 3; // D: if(z < k && z > 0) { if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) { nActivatedPads++; nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX - 1; eff[nActivatedPads-1] = effX * effZ; (resZGaus(-fPulseHeightSlope * z, fLogChargeSmearing); } else { qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.); qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.); logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing); } timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing); } else { timeDelay[nActivatedPads-1] = 0.; } padId[nActivatedPads-1] = 4; } } // end D } // end C // E: if(ix < AliTOFConstants::fgkNpadX && dX > 0) { nActivatedPads++; nPlace[nActivatedPads-1] = nPlace[0] + 1; eff[nActivatedPads-1] = effX; res[nActivatedPads-1] = 0.001 * (TMath::Sqrt(10400 + resX * resX)); // ns timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns nTail[nActivatedPads-1] = 2; if (fTimeDelayFlag) { qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.); qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.); logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing); timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing); } else { timeDelay[nActivatedPads-1] = 0.; } padId[nActivatedPads-1] = 5; // F: if(z < k && z > 0) { if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) { nActivatedPads++; nPlace[nActivatedPads - 1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX + 1; eff[nActivatedPads - 1] = effX * effZ; (resZGaus(-fPulseHeightSlope * z, fLogChargeSmearing); } else { qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.); qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.); logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing); } timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing); } else { timeDelay[nActivatedPads-1] = 0.; } padId[nActivatedPads-1] = 6; } } // end F } // end E } // end if(x < k) for (Int_t iPad = 0; iPad < nActivatedPads; iPad++) { if (res[iPad] < fTimeResolution) res[iPad] = fTimeResolution; if(gRandom->Rndm() < eff[iPad]) { isFired[iPad] = kTRUE; nFiredPads++; if(fEdgeTails) { if(nTail[iPad] == 0) { tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]); } else { ftail->SetParameters(res[iPad], 2. * res[iPad], kSigmaForTail[nTail[iPad]-1]); Double_t timeAB = ftail->GetRandom(); tofTime[iPad] = geantTime + timeWalk[iPad] + timeDelay[iPad] + timeAB; } } else { tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]); } if (fAverageTimeFlag) { averageTime += tofTime[iPad] * qInduced[iPad]; weightsSum += qInduced[iPad]; } else { averageTime += tofTime[iPad]; weightsSum += 1.; } } } if (weightsSum!=0) averageTime /= weightsSum; } // end else (fEdgeEffect != 0) //cout << "timedelay " << timeDelay[0] << endl; //cout << "timedelay " << timeDelay[1] << endl; //cout << "timedelay " << timeDelay[2] << endl; //cout << "timedelay " << timeDelay[3] << endl; } */ //__________________________________________________________________ Int_t AliTOFReconstructioner::PDGtoGeantCode(Int_t pdgcode) { // // Gives the GEANT code from KF code of LUND JETSET // Int_t geantCode=0; // default value switch (pdgcode) { case 22: geantCode=1; // GAMMA break ; case -11: geantCode=2; // E+ break ; case 11: geantCode=3; // E- break ; case 12: geantCode=4; // NUE break ; case 14: geantCode=4; // NUMU break ; case -13: geantCode=5; // MU+ break ; case 13: geantCode=6; // MU- break ; case 111: geantCode=7; // PI0 break ; case 211: geantCode=8; // PI+ break ; case -211: geantCode=9; // PI- break ; case 130: geantCode=10; // K_L0 break ; case 321: geantCode=11; // K+ break ; case -321: geantCode=12; // K- break ; case 2112: geantCode=13; // N0 break ; case 2212: geantCode=14; // P+ break ; case -2212: geantCode=15; // P~- break ; case 310: geantCode=16; // K_S0 break ; case 221: geantCode=17; // ETA break ; case 3122: geantCode=18; // LAMBDA0 break ; case 3222: geantCode=19; // SIGMA+ break ; case 3212: geantCode=20; // SIGMA0 break ; case 3112: geantCode=21; // SIGMA- break ; case 3322: geantCode=22; // XI0 break ; case 3312: geantCode=23; // XI- break ; case 3334: geantCode=24; // OMEGA- break ; case -2112: geantCode=25; // N~0 break ; case -3122: geantCode=26; // LAMBDA~0 break ; case -3112: geantCode=27; // SIGMA~+ break ; case -3212: geantCode=28; // SIGMA~0 break ; case -3222: geantCode=29; // SIGMA~- break ; case -3322: geantCode=30; // XI~0 break ; case -3312: geantCode=31; // XI~+ break ; case -3334: geantCode=32; // OMEGA~+ break ; case 223: geantCode=33; // OMEGA(782) break ; case 333: geantCode=34; // PHI(1020) break ; case 411: geantCode=35; // D+ break ; case -411: geantCode=36; // D- break ; case 421: geantCode=37; // D0 break ; case -421: geantCode=38; // D~0 break ; case 431: geantCode=39; // D_S+ break ; case -431: geantCode=40; // D_S~- break ; case 4122: geantCode=41; // LAMBDA_C+ break ; case 213: geantCode=42; // RHP(770)+ break ; case -213: geantCode=43; // RHO(770)- break ; case 113: geantCode=44; // RHO(770)0 break ; default: geantCode=45; break; } return geantCode; } //__________________________________________________________________ Bool_t AliTOFReconstructioner::operator==( AliTOFReconstructioner const & tofrec)const { // Equal operator. // Reconstructioners are equal if their parameters are equal // split the member variables in analogous categories // time resolution and edge effect parameters Bool_t dummy0=(fTimeResolution==tofrec.fTimeResolution)&& (fpadefficiency==tofrec.fpadefficiency)&& (fEdgeEffect==tofrec.fEdgeEffect)&& (fEdgeTails==tofrec.fEdgeTails)&& (fHparameter==tofrec.fHparameter)&& (fH2parameter==tofrec.fH2parameter)&& (fKparameter==tofrec.fKparameter)&& (fK2parameter==tofrec.fK2parameter); // pad efficiency parameters Bool_t dummy1=(fEffCenter==tofrec.fEffCenter)&& (fEffBoundary==tofrec.fEffBoundary)&& (fEff2Boundary==tofrec.fEff2Boundary)&& (fEff3Boundary==tofrec.fEff3Boundary)&& (fResCenter==tofrec.fResCenter)&& (fResBoundary==tofrec.fResBoundary)&& (fResSlope==tofrec.fResSlope); // time walk parameters Bool_t dummy2=(fTimeWalkCenter==tofrec.fTimeWalkCenter)&& (fTimeWalkBoundary==tofrec.fTimeWalkBoundary)&& (fTimeWalkSlope==tofrec.fTimeWalkSlope)&& (fTimeDelayFlag==tofrec.fTimeDelayFlag)&& (fPulseHeightSlope==tofrec.fPulseHeightSlope)&& (fTimeDelaySlope==tofrec.fTimeDelaySlope); // ADC-TDC correlation parameters Bool_t dummy3=(fMinimumCharge==tofrec.fMinimumCharge)&& (fChargeSmearing==tofrec.fChargeSmearing )&& (fLogChargeSmearing==tofrec.fLogChargeSmearing )&& (fTimeSmearing==tofrec.fTimeSmearing )&& (fAverageTimeFlag==tofrec.fAverageTimeFlag)&& (fChargeFactorForMatching==tofrec.fChargeFactorForMatching)&& (fMatchingStyle==tofrec.fMatchingStyle); Bool_t dummy4=(fTrackingEfficiency==tofrec.fTrackingEfficiency)&& (fSigmavsp==tofrec.fSigmavsp)&& (fSigmaZ==tofrec.fSigmaZ)&& (fSigmarphi==tofrec.fSigmarphi)&& (fSigmap==tofrec.fSigmap)&& (fSigmaPhi==tofrec.fSigmaPhi)&& (fSigmaTheta==tofrec.fSigmaTheta)&& (fNoise==tofrec.fNoise)&& (fNoiseSlope==tofrec.fNoiseSlope)&& (fField==tofrec.fField)&& (fRadLenTPC==tofrec.fRadLenTPC)&& (fCorrectionTRD==tofrec.fCorrectionTRD)&& (fLastTPCRow==tofrec.fLastTPCRow)&& (fRadiusvtxBound==tofrec.fRadiusvtxBound)&& (fMaxTestTracks==tofrec.fMaxTestTracks)&& (fStep==tofrec.fStep)&& (fMaxPixels==tofrec.fMaxPixels)&& (fMaxAllTracks==tofrec.fMaxAllTracks)&& (fMaxTracks==tofrec.fMaxTracks)&& (fMaxTOFHits==tofrec.fMaxTOFHits)&& (fPBound==tofrec.fPBound); if( dummy0 && dummy1 && dummy2 && dummy3 && dummy4) return kTRUE ; else return kFALSE ; } //____________________________________________________________________________ void AliTOFReconstructioner::UseHitsFrom(const char * filename) { SetTitle(filename) ; } //____________________________________________________________________________ void AliTOFReconstructioner::InitArray(Float_t array[], Int_t nlocations) { // // Initialize the array of Float_t // for (Int_t i = 0; i < nlocations; i++) { array[i]=0.; } // end loop } //____________________________________________________________________________ void AliTOFReconstructioner::InitArray(Int_t array[], Int_t nlocations) { // // Initialize the array of Int_t // for (Int_t i = 0; i < nlocations; i++) { array[i]=0; } // end loop } //____________________________________________________________________________ void AliTOFReconstructioner::ReadTOFHits(Int_t ntracks, TTree* treehits, TClonesArray* tofhits, Int_t ***MapPixels, Int_t* kTOFhitFirst, AliTOFPad* pixelArray , Int_t* iTOFpixel, Float_t* toftime, AliTOFRecHit* hitArray, Int_t& isHitOnFiredPad, Int_t& ipixel) { // // Read TOF hits for the current event and fill arrays // // Start loop on primary tracks in the hits containers // // Noise meaning in ReadTOFHits: we use the word 'noise' in the // following cases // - signals produced by secondary particles // - signals produced by the next hits (out of the first) of a given track // (both primary and secondary) // - signals produced by edge effect TParticle *particle; Int_t nHitOutofTofVolumes; // number of hits out of TOF GEANT volumes (it happens in very // few cases) Int_t * npixel = new Int_t[AliTOFConstants::fgkmaxtoftree]; // array used by TOFRecon for check on TOF geometry Int_t npions=0; // number of pions for the current event Int_t nkaons=0; // number of kaons for the current event Int_t nprotons=0; // number of protons for the current event Int_t nelectrons=0;// number of electrons for the current event Int_t nmuons=0; // number of muons for the current event Float_t tofpos[3]; // TOF hit position and GEANT time Float_t zPad,xPad; Int_t nbytes = 0; Int_t ipart, nhits=0, nHitsFromPrimaries=0; Int_t ntotalTOFhits=0; // total number of TOF hits for the current event Int_t ipartLast=-1; // last track identifier Int_t iFirstHit; // flag to check if the current hit is the first hit on TOF for the // current track Int_t iNoiseHit=0; // flag used to tag noise hits (the noise meaning is reported in the // header of the ReadTOFHits method) Int_t nhitWithoutNoise;// number of hits not due to noise Int_t inoise=0,inoise2=0; Int_t nMultipleSignOnSamePad=0; // number of cases where a pad is fired more than one time Int_t nPixEdge=0; // additional pads fired due to edge effect in ReadTOFHits (local var) // array used for counting different types of primary particles Int_t particleTypeGEANT[50]={0,4,4,0,5,5,0,3,3,0, 2,2,0,1,1,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0}; Int_t particleType,particleInTOFtype[6][3]; for (Int_t i=0;i<6;i++) { for (Int_t j=0;j<3;j++) { particleInTOFtype[i][j]=0; } } // speed-up the code treehits->SetBranchStatus("*",0); // switch off all branches treehits->SetBranchStatus("TOF*",1); // switch on only TOF for (Int_t track=0; trackResetHits(); nbytes += treehits->GetEvent(track); nhits = tofhits->GetEntriesFast(); ntotalTOFhits+=nhits; // Start loop on hits connected to the current primary tracked particle // (including hits produced by secondary particles generaterd by the // current ptimary tracked particle) for (Int_t hit=0;hitUncheckedAt(hit); ipart = tofHit->GetTrack(); if(ipart>=fMaxAllTracks) break; Float_t geantTime= tofHit->GetTof(); // it is given in [s] particle = (TParticle*)gAlice->GetMCApp()->Particle(ipart); Int_t pdgCode=particle->GetPdgCode(); // Only high momentum tracks (see fPBound value) // momentum components at vertex Float_t pxvtx = particle->Px(); Float_t pyvtx = particle->Py(); Float_t pzvtx = particle->Pz(); Float_t pvtx = TMath::Sqrt(pxvtx*pxvtx+pyvtx*pyvtx+pzvtx*pzvtx); if(pvtx>fPBound) { if(particle->GetFirstMother() < 0) nHitsFromPrimaries++; // count primaries // x and y coordinates of the particle production vertex Float_t vx = particle->Vx(); Float_t vy = particle->Vy(); Float_t vr = TMath::Sqrt(vx*vx+vy*vy); // cylindrical radius of the particle production vertex Float_t x = tofHit->X(); tofpos[0]=x; Float_t y = tofHit->Y(); tofpos[1]=y; Float_t z = tofHit->Z(); tofpos[2]=z; /* var used for QA Float_t tofradius = TMath::Sqrt(x*x+y*y); // radius cilindrical coordinate of the TOF hit */ // momentum components (cosine) when striking the TOF Float_t pxtof = tofHit->GetPx(); Float_t pytof = tofHit->GetPy(); Float_t pztof = tofHit->GetPz(); // scalar product indicating the direction of the particle when striking the TOF /* var used for QA // (>0 for outgoing particles) Float_t isGoingOut = (x*pxtof+y*pytof+z*pztof)/TMath::Sqrt(x*x+y*y+z*z); */ Float_t momtof = tofHit->GetMom(); // now momentum components when striking the TOF pxtof *= momtof; pytof *= momtof; pztof *= momtof; particleType=particleTypeGEANT[PDGtoGeantCode(pdgCode)-1]; if(particleType) { particleInTOFtype[5][2]++; particleInTOFtype[particleType-1][2]++; } iFirstHit=0; // without noise hits if(ipart!=ipartLast) { iFirstHit=1; toftime[ipart]=geantTime; //time [s] // tofMom[ipart]=momtof; ipartLast=ipart; if(particle->GetFirstMother() < 0) { Int_t abspdgCode=TMath::Abs(pdgCode); switch (abspdgCode) { case 211: npions++; break ; case 321: nkaons++; break ; case 2212: nprotons++; break ; case 11: nelectrons++; break ; case 13: nmuons++; break ; } } if(vr>fRadiusvtxBound) { if(particleType) { particleInTOFtype[5][1]++; particleInTOFtype[particleType-1][1]++; } inoise++; inoise2++; } else { if(particleType) { particleInTOFtype[5][0]++; particleInTOFtype[particleType-1][0]++; } } } else { inoise++; if(particleType) { particleInTOFtype[5][1]++; particleInTOFtype[particleType-1][1]++; } } //end if(ipart!=ipartLast) IsInsideThePad(gMC,x,y,z,npixel,zPad,xPad); Int_t sec = tofHit->GetSector(); Int_t pla = tofHit->GetPlate(); Int_t str = tofHit->GetStrip(); if(sec!=npixel[0] || pla!=npixel[1] || str!=npixel[2]){// check on volume cout << "sector" << sec << " npixel[0] " << npixel[0] << endl; cout << "plate " << pla << " npixel[1] " << npixel[1] << endl; cout << "strip " << str << " npixel[2] " << npixel[2] << endl; } // close check on volume Int_t padz = tofHit->GetPadz(); Int_t padx = tofHit->GetPadx(); Float_t Zpad = tofHit->GetDz(); Float_t Xpad = tofHit->GetDx(); if (npixel[4]==0){ IsInsideThePad(gMC,x,y,z,npixel,zPad,xPad); if (npixel[4]==0){ nHitOutofTofVolumes++; } } else { Float_t zStrip=AliTOFConstants::fgkZPad*(padz-0.5-0.5*AliTOFConstants::fgkNpadZ)+Zpad; if(padz!=npixel[3]) printf(" : Zpad=%f, padz=%i, npixel[3]=%i, zStrip=%f\n",Zpad,padz,npixel[3],zStrip); Float_t xStrip=AliTOFConstants::fgkXPad*(padx-0.5-0.5*AliTOFConstants::fgkNpadX)+Xpad; Int_t nPlace[4]={0,0,0,0}; nPlace[0]=(padz-1)*AliTOFConstants::fgkNpadX+padx; Int_t nActivatedPads=0; Int_t nFiredPads=0; Bool_t isFired[4]={kFALSE,kFALSE,kFALSE,kFALSE}; Float_t tofAfterSimul[4]={0.,0.,0.,0.}; Float_t qInduced[4]={0.,0.,0.,0.}; Float_t averageTime=0.; BorderEffect(zStrip,xStrip,geantTime*1.0e+09,nActivatedPads,nFiredPads,isFired,nPlace,qInduced,tofAfterSimul,averageTime); // simulate edge effect if(nFiredPads) { for(Int_t indexOfPad=0; indexOfPadfRadiusvtxBound || iFirstHit==0) iNoiseHit=1; hitArray[isHitOnFiredPad-1].SetNoise(iNoiseHit); if(iFirstHit) kTOFhitFirst[ipart]=isHitOnFiredPad; }// close - the hit belongs to a fired pad Int_t iMapFirstIndex=AliTOFConstants::fgkNSectors*(npixel[1]-1)+npixel[0]-1; Int_t iMapValue=MapPixels[iMapFirstIndex][npixel[2]-1][nPlace[indexOfPad]-1]; if(iMapValue==0) { ipixel++; if(indexOfPad) { iNoiseHit=1; nPixEdge++; } else { iTOFpixel[ipart]=ipixel; } if(ipixel>fMaxPixels){ // check on the total number of activated pads cout << "ipixel=" << ipixel << " > fMaxPixels=" << fMaxPixels << endl; return; } // close check on the number of activated pads MapPixels[iMapFirstIndex][npixel[2]-1][nPlace[indexOfPad]-1]=ipixel; pixelArray[ipixel-1].SetGeom(npixel[0],npixel[1],npixel[2],nPlace[indexOfPad]); pixelArray[ipixel-1].SetTrack(ipart); if(iNoiseHit) { pixelArray[ipixel-1].AddState(1); } else { if(tofAfterSimul[indexOfPad]<0) cout << "Time of Flight after detector simulation is negative" << endl; pixelArray[ipixel-1].AddState(10); } pixelArray[ipixel-1].SetTofChargeHit(tofAfterSimul[indexOfPad],qInduced[indexOfPad],geantTime*1.0e+09,isHitOnFiredPad); } else { //else if(iMapValue==0) if(indexOfPad==0) iTOFpixel[ipart]=iMapValue; nMultipleSignOnSamePad++; if(tofAfterSimul[indexOfPad] < pixelArray[iMapValue-1].GetRealTime() ) { pixelArray[iMapValue-1].SetTrack(ipart); // if(indexOfPad==0) pixelArray[iMapValue-1].SetTrack(ipart); if(indexOfPad) iNoiseHit=1; if(iNoiseHit) { pixelArray[iMapValue-1].AddState(1); } else { pixelArray[iMapValue-1].AddState(10); } pixelArray[iMapValue-1].SetRealTime(tofAfterSimul[indexOfPad]); pixelArray[iMapValue-1].SetGeantTime(geantTime*1.0e+09); pixelArray[iMapValue-1].SetHit(isHitOnFiredPad); } // close if(tofAfterSimul[indexOfPad] < pixelArray[iMapValue-1].GetTime() ) } //end of Pixel filling } // close if(isFired[indexOfPad]) } //end loop on activated pads indexOfPad } // close if(nFiredPads) } //end of hit with npixel[3]!=0 } //high momentum tracks } //end on TOF hits } //end on primary tracks if(fdbg) { cout << ntotalTOFhits << " - total number of TOF hits " << nHitsFromPrimaries << " - primary " << endl; cout << inoise << " - noise hits, " << inoise2<< " - first crossing of a track with Rvtx>" << fRadiusvtxBound << endl; // cout << inoise << " - noise hits (" << 100.*inoise/ihit << " %), " << inoise2 //<< " - first crossing of a track with Rvtx>" << RVTXBOUND << endl; nhitWithoutNoise=isHitOnFiredPad; cout << ipixel << " fired pixels (" << nMultipleSignOnSamePad << " multiple fired pads, " << endl; //j << " fired by noise, " << j1 << " noise+track)" << endl; printf(" %i additional pads are fired due to edge effect\n",nPixEdge); cout << npions << " primary pions reached TOF" << endl; cout << nkaons << " primary kaons reached TOF" << endl; cout << nprotons << " primary protons reached TOF" << endl; cout << nelectrons<<" primary electrons reached TOF" << endl; cout << nmuons << " primary muons reached TOF" << endl; cout << "number of TOF hits for different species: 1-p, 2-K, 3-pi, 4-e, 5-mu, 6-all" << endl; cout << " first number - track hits, second - noise ones, third - all" << endl; for (Int_t i=0;i<6;i++) cout << i+1 << " " << particleInTOFtype[i][0] << " " << particleInTOFtype[i][1] << " " << particleInTOFtype[i][2] << endl; Int_t primaryReachedTOF[6]; primaryReachedTOF[0]=npions; primaryReachedTOF[1]=nkaons; primaryReachedTOF[2]=nprotons; primaryReachedTOF[3]=nelectrons; primaryReachedTOF[4]=nmuons; primaryReachedTOF[5]=npions+nkaons+nprotons+nelectrons+nmuons; cout << " Reading TOF hits done" << endl; } delete [] npixel; } //____________________________________________________________________________ void AliTOFReconstructioner::AddNoiseFromOuter(Option_t *option, Int_t ***MapPixels, AliTOFPad* pixelArray , AliTOFRecHit* hitArray, Int_t& isHitOnFiredPad, Int_t& ipixel) { // // Add noise hits from outer regions (forward and backward) according // to parameterized fZNoise distribution (to be used with events // generated in the barrel region) Float_t * zLen = new Float_t[AliTOFConstants::fgkNPlates+1]; Float_t * zStrips = new Float_t[AliTOFConstants::fgkNPlates]; zStrips[0]=(Float_t) (AliTOFConstants::fgkNStripC); zStrips[1]=(Float_t) (AliTOFConstants::fgkNStripB); zStrips[2]=(Float_t) (AliTOFConstants::fgkNStripA); zStrips[3]=(Float_t) (AliTOFConstants::fgkNStripB); zStrips[4]=(Float_t) (AliTOFConstants::fgkNStripC); zLen[5]=AliTOFConstants::fgkzlenA*0.5+AliTOFConstants::fgkzlenB+AliTOFConstants::fgkzlenC; zLen[4]=zLen[5]-AliTOFConstants::fgkzlenC; zLen[3]=zLen[4]-AliTOFConstants::fgkzlenB; zLen[2]=zLen[3]-AliTOFConstants::fgkzlenA; zLen[1]=zLen[2]-AliTOFConstants::fgkzlenB; zLen[0]=zLen[1]-AliTOFConstants::fgkzlenC; Int_t isector; // random sector number Int_t iplate; // random plate number Int_t istrip; // random strip number in the plate Int_t ipadAlongX; // random pad number along x direction Int_t ipadAlongZ; // random pad number along z direction Int_t ipad; Int_t nPixEdge=0; // additional pads fired due to edge effect when adding noise from outer // regions // x -> time of flight given in ns TF1 *noiseTof = new TF1("noiseTof","exp(-x/20)",0,100); if(strstr(option,"pp")){ fZnoise = new TF1("fZnoise","257.8-0.178*x-0.000457*x*x",-AliTOFConstants::fgkMaxhZtof,AliTOFConstants::fgkMaxhZtof); } if(strstr(option,"Pb-Pb")){ fZnoise = new TF1("fZnoise","182.2-0.09179*x-0.0001931*x*x",-AliTOFConstants::fgkMaxhZtof,AliTOFConstants::fgkMaxhZtof); } if(fNoise) { if(fdbg) cout << " Start adding additional noise hits from outer regions" << endl; for(Int_t i=0;iRndm())+1; //the sector number // non-flat z-distribution of additional hits Float_t zNoise=fZnoise->GetRandom(); // holes for PHOS and HMPID if(((AliTOF *) gAlice->GetDetector("TOF"))->IsVersion()==2) { // to be checked the holes case if(isector>12 && isector<16) { // sectors 13,14,15 - RICH do { iplate=(Int_t) (AliTOFConstants::fgkNPlates*gRandom->Rndm())+1; } while (iplate==2 || iplate==3 || iplate==4); // } else if(isector>11 && isector<17) { // sectors 12,13,14,15,16 - PHOS } else if(isector>2 && isector<8) { // sectors 3,4,5,6,7 - PHOS do { iplate=(Int_t) (AliTOFConstants::fgkNPlates*gRandom->Rndm())+1; } while (iplate==3); } else { iplate=(Int_t) (AliTOFConstants::fgkNPlates*gRandom->Rndm())+1; } } else { iplate=0; do { iplate++; } while(zNoise>zLen[iplate]); } // end of holes if(iplate<1 || iplate>5) { printf(" iplate<1 or iplate>5, iplate=%i\n",iplate); return; } Float_t nStripes=0; if(iplate>1) { for (Int_t i=0;iRndm())+1; ipadAlongZ = (Int_t)(AliTOFConstants::fgkNpadZ*gRandom->Rndm())+1; ipad=(Int_t)(ipadAlongZ-1)*AliTOFConstants::fgkNpadX+ipadAlongX; //the pad number Float_t xStrip=(ipadAlongX-1)*AliTOFConstants::fgkXPad+AliTOFConstants::fgkXPad*gRandom->Rndm()-0.5*AliTOFConstants::fgkNpadX*AliTOFConstants::fgkXPad;//x-coor.in the strip frame Float_t zStrip=(ipadAlongZ-1)*AliTOFConstants::fgkZPad+AliTOFConstants::fgkZPad*gRandom->Rndm()-0.5*AliTOFConstants::fgkNpadZ*AliTOFConstants::fgkZPad;//z-coor.in the strip frame Int_t nPlace[4]={0,0,0,0}; nPlace[0]=ipad; Int_t nActivatedPads=0; Int_t nFiredPads=0; Bool_t isFired[4]={kFALSE,kFALSE,kFALSE,kFALSE}; Float_t tofAfterSimul[4]={0.,0.,0.,0.}; Float_t qInduced[4]={0.,0.,0.,0.}; Float_t averageTime=0.; Float_t toffornoise=10.+noiseTof->GetRandom(); // 10 ns offset + parameterization [ns] BorderEffect(zStrip,xStrip,toffornoise,nActivatedPads,nFiredPads,isFired,nPlace,qInduced,tofAfterSimul,averageTime); // simulate edge effect if(nFiredPads) { for(Int_t indexOfPad=0; indexOfPad10) nAll++; } } if(fdbg) { cout << " after adding " << fNoise << " noise hits: " << ipixel << " fired pixels (" << nNoiseSignals << " fired by noise, " << nAll << " noise+track)" << endl; printf(" %i additional pixels are fired by noise due to edge effect\n",nPixEdge); cout << " End of adding additional noise hits from outer regions" << endl; } Float_t occupancy; // numberOfPads for AliTOFV4 (Full coverage) // - to be upgraded checking the used TOF version - Float_t numberOfPads=AliTOFConstants::fgkPadXSector*AliTOFConstants::fgkNSectors; occupancy=100.*ipixel/numberOfPads; // percentage of fired pads printf(" Overall TOF occupancy (percentage of fired pads after adding noise) = %f\n",occupancy); delete [] zLen; delete [] zStrips; } //____________________________________________________________________________ void AliTOFReconstructioner::SetMinDistance(AliTOFRecHit* hitArray, Int_t ilastEntry) { // // Set the distance to the nearest hit for hitArray // ilastEntry is the index of the last entry of hitArray // starting the setting for the distance to the nearest TOFhit (cm) for(Int_t i=0; iSetBranchStatus("*",0); // switch off all branches treehits->SetBranchStatus("TPC*",1); // switch on only TPC for (Int_t track=0; trackResetHits(); nbytes += treehits->GetEvent(track); nhits = tpchits->GetEntriesFast(); for (Int_t hit=0;hitUncheckedAt(hit); Int_t row = tpcHit->fPadRow; ipart = tpcHit->GetTrack(); if(ipart>=fMaxAllTracks) break; particle = (TParticle*)gAlice->Particle(ipart); Int_t pdgCode=particle->GetPdgCode(); // only high momentum tracks // momentum components at production vertex Float_t pxvtx = particle->Px(); Float_t pyvtx = particle->Py(); Float_t pzvtx = particle->Pz(); Float_t pvtx = TMath::Sqrt(pxvtx*pxvtx+pyvtx*pyvtx+pzvtx*pzvtx); if(pvtx>fPBound && row == fLastTPCRow) { Float_t vx = particle->Vx(); Float_t vy = particle->Vy(); Float_t vr = TMath::Sqrt(vx*vx+vy*vy); Float_t x = tpcHit->X(); Float_t y = tpcHit->Y(); Float_t z = tpcHit->Z(); pos[0]=x; pos[1]=y; pos[2]=z; Float_t pxtpc = tpcHit->fPx; Float_t pytpc = tpcHit->fPy; Float_t pztpc = tpcHit->fPz; mom[0]=pxtpc; mom[1]=pytpc; mom[2]=pztpc; Float_t momtpc = TMath::Sqrt(pxtpc*pxtpc+pytpc*pytpc+pztpc*pztpc); if(x*pxtpc+y*pytpc>0) { // only tracks going out of TPC Float_t isoutgoing = x*pxtpc+y*pytpc+z*pztpc; isoutgoing /= (momtpc*TMath::Sqrt(x*x+y*y+z*z)); tpclen = tpcHit->fL; if(ipart!=idummy) { if(particle->GetFirstMother() < 0) { Int_t abspdgCode=TMath::Abs(pdgCode); switch (abspdgCode) { case 211: npions++; break ; case 321: nkaons++; break ; case 2212: nprotons++; break ; case 11: nelectrons++; break ; case 13: nmuons++; break ; } } // close if(particle->GetFirstMother() < 0) } // close if(ipart!=idummy) if(gRandom->Rndm()GetFirstMother() < 0) iprim++; if(itrack>fMaxTracks) { cout << "itrack=" << itrack << " > MAXTRACKS=" << fMaxTracks << endl; return; } // close if(itrack>fMaxTracks) iparticle[ipart]=itrack; trackArray[itrack-1].SetTrack(ipart,pvtx,pdgCode,tpclen,pos,mom,trdpos); pt=TMath::Sqrt(pxtpc*pxtpc+pytpc*pytpc); // pt: transverse momentum at TPC // Filling iTrackPt[MAXTRACKS] by itrack ordering on Pt if(itrack==1) { iTrackPt[itrack-1]=itrack; ptTrack[itrack-1]=pt; } else { for (Int_t i=0; iptTrack[i]) { for(Int_t j=i; jfRadiusvtxBound) nTpcHitsLastRow++; if(ipart==idummy) nTpcDoubleHitsLastRow++; idummy=ipart; } // close if(x*px+y*py>0) } // close if(pvtx>fPBound && row == fLastTPCRow) } //end of hits } // close loop on tracks if(fdbg) { cout << ntotalTPChits << " TPC hits in the last TPC row " << fLastTPCRow << endl; cout << " " << nTpcHitsLastRow << " - hits with Rvtx>fRadiusvtxBound=" << fRadiusvtxBound << endl; cout << " " << nTpcDoubleHitsLastRow << " double TPC hits" << endl; cout << itrack << " - extracted TPC tracks " << iprim << " - primary" << endl; cout << npions << " primary pions reached TPC" << endl; cout << nkaons << " primary kaons reached TPC" << endl; cout << nprotons << " primary protons reached TPC" << endl; cout << nelectrons<< " primary electrons reached TPC" << endl; cout << nmuons << " primary muons reached TPC" << endl; } // if(fdbg) Int_t primaryInTPC[6]={0,0,0,0,0,0}; primaryInTPC[0]=npions; primaryInTPC[1]=nkaons; primaryInTPC[2]=nprotons; primaryInTPC[3]=nelectrons; primaryInTPC[4]=nmuons; primaryInTPC[5]=npions+nkaons+nprotons+nelectrons+nmuons; if(fdbg) { printf(" contents of iTrackPt[MAXTRACKS],PtTrack[MAXTRACKS]\n"); for (Int_t i=0; i=0; i--) { printf(" %i : iTrackPt=%i, PtTrack=%f\n",i+1,iTrackPt[i],ptTrack[i]); } }// if(fdbg) } */ //____________________________________________________________________________ void cylcor(Float_t& x, Float_t& y) { Float_t rho,phi; rho=TMath::Sqrt(x*x+y*y); phi=0.; if(TMath::Abs(x)>0. || TMath::Abs(y)>0.) phi=TMath::ATan2(y,x); if(phi<0.) phi=phi+2.*TMath::Pi(); x=rho; y=phi; } //____________________________________________________________________________ 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) { Int_t TestTracks,iTestTrack,itest,wPixel=0,itestc; Int_t * ntest = new Int_t[fMaxTestTracks]; Int_t * testPixel = new Int_t[fMaxTestTracks]; Float_t wLength=0.,wRho=0.,wZ=0.; Float_t * testLength = new Float_t[fMaxTestTracks]; Float_t * testRho = new Float_t[fMaxTestTracks]; Float_t * testZ = new Float_t[fMaxTestTracks]; Float_t weight; Float_t * testWeight = new Float_t[fMaxTestTracks]; Float_t rotationFactor,phi0,coslam,sinlam,helixRadius,xHelixCenter,yHelixCenter,zHelixCenter,helixFactor; Int_t npixel[5],iMapValue,iwork1,iwork2,iwork3,iwork4,ihit=0; Int_t charge[48]={ 0, 1,-1, 0, 1,-1, 0, 1,-1, 0, 1,-1, 0, 1,-1, 0, 0, 0, 1, 0, -1, 0,-1,-1, 0, 0,-1, 0, 1, 0, 1, 1, 0, 0, 1,-1, 0, 0, 1,-1, 1, 1,-1, 0, 1, 1, 2, 0}; Float_t theta0,gpx,gpy,gpz,gp,gpt,gtheta,gx,gy,gz,gr,gxLast,gyLast,gzLast,chargeField; Float_t sumOfTheta=0.,weightTestTracksOutTof[4]; Float_t s,ds,xRespectToHelixCenter,yRespectToHelixCenter,deltaRadius,fp,xp,yp,grho; Float_t mass,energy,g; Int_t itrack=0,itr,particleCharge,istep,iplate=0,iPadAlongX=0; Int_t itra,t34=0,t32=0,t44=0,t43=0,t42=0; Int_t wstate=0,m2state=0,wPix; Int_t idelR=0,idelR1=0,idelR2=0,iRmin=0,iRmin1=0,iRmin2=0; Float_t massArray[50] = {0.0,0.00051,0.00051,0.0,0.1057,0.1057,0.135,0.1396,0.1396,0.4977, 0.4936,0.4936,0.9396,0.9383,0.9383,0.4977,0.5488,1.1156,1.1894,1.1926,1.1926, 1.3149,1.3213,1.6724,0.9396,1.1156,1.1894,1.1926,1.1974,1.3149, 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; Float_t delR; Float_t radius,area,normR,normS,cosAngl; Int_t iPlateFirst,iTestGmax=0; Int_t fstate,iPrintM1=0,iPrintM2=0; Float_t gxExtrap=0.,gyExtrap=0.,gzExtrap=0.; Float_t avSigZ=0,avSigRPHI=0,avSigP=0,avSigPHI=0,avSigTHETA=0; Float_t gxW,gyW,gzW; Float_t length0; Float_t snr=0; Int_t indexOfTestTrack; Float_t zPad,xPad; Int_t istate=0,imax=0,match,iMaxTestTracksOutTof=0,matchw; Float_t w,wmax=0.,inverseOfParticleSpeed,w2,smat[9],largestWeightTracksOutTof,sw; Float_t sumWeightTracksOutTof,sGeomWeigth; Int_t imatched; Int_t m10=0,m20=0,m22=0,m23=0; Int_t PRINT=0; TParticle *particle; Float_t time=0.; itr=ntotTpcTracks; printf(" itr=%i\n",itr); for (itra=1; itra0. || TMath::Abs(y)>0.) phi=TMath::ATan2(y,x); if(phi<0.) phi=phi+2.*TMath::Pi(); /* var used for QA Float_t phiTPC=phi*kRaddeg; */ if(fSigmavsp) { if(p==0) printf(" p=%f in g=0.022/p\n",p); g=0.022/p; avSigRPHI += g; // (cm) if(rho==0) printf(" rho=%f in phi += g*gRandom->Gaus()/rho\n",rho); phi += g*gRandom->Gaus()/rho; } else { if(rho==0) printf(" rho=%f in phi += (SIGMARPHI*gRandom->Gaus()/rho\n",rho); phi += (fSigmarphi*gRandom->Gaus()/rho); } x=rho*TMath::Cos(phi); y=rho*TMath::Sin(phi); /* var used for QA Float_t zTPC=z; */ if(fSigmavsp) { if(p==0) printf(" p=%f in g=0.0275/p\n",p); g=0.0275/p; avSigZ += g; // (cm) z += g*gRandom->Gaus(); } else { z += fSigmaZ*gRandom->Gaus(); } // smearing on TPC momentum { Float_t pmom,phi,theta,arg; pmom=TMath::Sqrt(px*px+py*py+pz*pz); phi=0.; if(TMath::Abs(px)>0. || TMath::Abs(py)>0.) phi=TMath::ATan2(py,px); if(phi<0.) phi=phi+2*TMath::Pi(); arg=1.; if(pmom>0.) arg=pz/pmom; theta=0.; if(TMath::Abs(arg)<=1.) theta=TMath::ACos(arg); if(fSigmavsp) { if(pmom<=0) printf(" pmom=%f in g = TMath::Abs(TMath::Log(pmom)/TMath::Log(10)+0.5)/0.7\n",pmom); g = TMath::Abs(TMath::Log(pmom)/TMath::Log(10)+0.5)/0.7; g = 0.01*(g*g*g+1.5)*1.24; avSigP += g; pmom *= (1+g*gRandom->Gaus()); if(p<10) { if(pmom<=0) printf(" pmom=%f in g = 1-TMath::Log(pmom)/TMath::Log(10)\n",pmom); g = 1-TMath::Log(pmom)/TMath::Log(10); g = 0.001*(g*g*g+0.3)*0.65; // (radian) } else { g = 0.001*0.3*0.65; } avSigPHI += g; phi += g*gRandom->Gaus(); avSigTHETA += g; theta += g*gRandom->Gaus(); } else { pmom *= (1+fSigmap*gRandom->Gaus()); phi += fSigmaPhi*gRandom->Gaus(); theta += fSigmaTheta*gRandom->Gaus(); } gxW=px; gyW=py; gzW=pz; px=pmom*TMath::Sin(theta)*TMath::Cos(phi); py=pmom*TMath::Sin(theta)*TMath::Sin(phi); pz=pmom*TMath::Cos(theta); if(x*px+y*py<=0) { x=Wx; y=Wy; z=Wz; px=gxW; py=gyW; pz=gzW; }// if(x*px+y*py<=0) } p = TMath::Sqrt(px*px+py*py+pz*pz); particleCharge=charge[PDGtoGeantCode(pdgCode)-1]; mass=massArray[PDGtoGeantCode(pdgCode)-1]; mass=massArray[8-1]; //we take pion mass for all tracks // mass=massArray[14-1]; //here we take proton mass for all tracks energy=TMath::Sqrt(p*p+mass*mass); chargeField=particleCharge*fField; g=fRadLenTPC/( (x*px+y*py)/(rho*p) ); 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); theta0=13.6*0.001*TMath::Sqrt(g)*(1.+0.038*TMath::Log(g))*energy/(p*p); // start Loop on test tracks sumOfTheta=0.; for(Int_t i=0;i<4;i++) { weightTestTracksOutTof[i]=0.; } itest=0; for(Int_t i=0;i0.01); if(istep==0) goto end; // Steps along the circle till a pad wPixel=0; wLength=0.; iplate=0; iPadAlongX=0; grho=0.; ds=fStep; gxLast=xHelixCenter+helixRadius*TMath::Cos(phi0+s*helixFactor); gyLast=yHelixCenter+helixRadius*TMath::Sin(phi0+s*helixFactor); gzLast=zHelixCenter+s*sinlam; do { istep++; s += ds; gx=xHelixCenter+helixRadius*TMath::Cos(phi0+s*helixFactor); gy=yHelixCenter+helixRadius*TMath::Sin(phi0+s*helixFactor); gz=zHelixCenter+s*sinlam; rho=TMath::Sqrt(gx*gx+gy*gy); IsInsideThePad(gMC,gx,gy,gz,npixel,zPad,xPad); iplate += npixel[1]; iPadAlongX += npixel[4]; if(indexOfTestTrack==0 && iplate && iPlateFirst==0) { iPlateFirst=1; length0=s; radius=s*3*theta0; area=TMath::Pi()*radius*radius; normR=TMath::Sqrt(gx*gx+gy*gy); normS=TMath::Sqrt((gx-gxLast)*(gx-gxLast)+ (gy-gyLast)*(gy-gyLast)+ (gz-gzLast)*(gz-gzLast)); cosAngl=(gx*(gx-gxLast)+gy*(gy-gyLast))/(normR*normS); if(cosAngl<0) printf(" cosAngl<0: gx=%f,gy=%f, gxLast=%f,gyLast=%f,gzLast=%f\n",gx,gy,gxLast,gyLast,gzLast); area /= cosAngl; TestTracks=(Int_t) (2*area/(AliTOFConstants::fgkXPad * AliTOFConstants::fgkZPad)); if(TestTracks<12) TestTracks=12; // Angles of entering into the TOF plate Int_t iZ=0; if(TMath::Abs(gz)>300) { iZ=4; } else if(TMath::Abs(gz)>200) { iZ=3; } else if(TMath::Abs(gz)>100) { iZ=2; } else if(TMath::Abs(gz)>0) { iZ=1; } } // end of if(indexOfTestTrack==0 && iplate && iPlateFirst==0) if(npixel[4]>0) { iwork1=npixel[0]; iwork2=npixel[1]; iwork3=npixel[2]; // iwork4=npixel[3]; iwork4=(npixel[3]-1)*AliTOFConstants::fgkNpadX+npixel[4]; Int_t ifirstindex=AliTOFConstants::fgkNSectors*(npixel[1]-1)+npixel[0]; iMapValue=mapPixels[ifirstindex-1][iwork3-1][iwork4-1]; if(iMapValue==0) { ipixel++; if(ipixel>fMaxPixels) { cout << "ipixel=" << ipixel << " > MAXPIXELS=" << fMaxPixels << endl; break; } mapPixels[ifirstindex-1][iwork3-1][iwork4-1]=ipixel; pixelArray[ipixel-1].SetGeom(iwork1,iwork2,iwork3,iwork4); iMapValue=ipixel; } wPixel=iMapValue; wLength=tpclength+s; wRho=rho; wZ=gz; ihit=kTOFhitFirst[ipart]; if(ihit) { if(indexOfTestTrack==0) { { idelR++; delR=TMath::Sqrt((gx-hitArray[ihit-1].X())*(gx-hitArray[ihit-1].X())+ (gy-hitArray[ihit-1].Y())*(gy-hitArray[ihit-1].Y())+ (gz-hitArray[ihit-1].Z())*(gz-hitArray[ihit-1].Z())); } if(delR>hitArray[ihit-1].GetRmin()) iRmin++; gxExtrap=gx; gyExtrap=gy; gzExtrap=gz; } else { delR=TMath::Sqrt((gx-gxExtrap)*(gx-gxExtrap)+ (gy-gyExtrap)*(gy-gyExtrap)+ (gz-gzExtrap)*(gz-gzExtrap)); } } //end of if(ihit) break; } //end of npixel[4] if(rho0) { if(iplate) { if(iPadAlongX==0) { istep=-3; // holes in TOF } } else { if(TMath::Abs(gz)0) { if(itest==0) { itest=1; ntest[itest-1]=1; testPixel[itest-1]=wPixel; testLength[itest-1]=wLength; testRho[itest-1]=wRho; testZ[itest-1]=wZ; testWeight[itest-1]=weight; } else { Int_t k=0; for(Int_t i=0;i-4 && istep<1) weightTestTracksOutTof[-istep] ++; } else { if(istep>-4 && istep<1) weightTestTracksOutTof[-istep] += weight; } if(fMatchingStyle==2) { if(indexOfTestTrack==0 && istep==0) break; if(indexOfTestTrack+1==TestTracks) break; } } //end of indexOfTestTrack snr += (Float_t) (indexOfTestTrack+1); // Search for the "hole" with the largest weigth largestWeightTracksOutTof=0.; sumWeightTracksOutTof=0.; for(Int_t i=0;i<4;i++) { w=weightTestTracksOutTof[i]; sumWeightTracksOutTof += w; if(w>largestWeightTracksOutTof) { largestWeightTracksOutTof=w; iMaxTestTracksOutTof=i; } } itestc=itest; if(itest>0) { for(Int_t i=0;i0) { fstate=1; wstate++; } if(fMatchingStyle==1) { sGeomWeigth += ntest[i]; w=(fpadefficiency*fstate+(1.-fpadefficiency)*(1-fstate))*ntest[i]; if(pixelArray[testPixel[i]-1].GetTrackMatched()>0) w *= 0.1; } else { sGeomWeigth += testWeight[i]; w=(fpadefficiency*fstate+(1.-fpadefficiency)*(1-fstate))*testWeight[i]; if(pixelArray[testPixel[i]-1].GetTrackMatched()>0) w *= 0.1; } // weighting according to the Pulse Height (we use the square of weight) // if (fChargeFactorForMatching) w *= (pixelArray[testPixel[i]-1].GetCharge())*(pixelArray[testPixel[i]-1].GetCharge()); if (fChargeFactorForMatching && fstate==1) w *= (pixelArray[testPixel[i]-1].GetCharge())*(pixelArray[testPixel[i]-1].GetCharge()); if(w>wmax) { wmax=w; imax=i; } sw += w; } wPixel=testPixel[imax]; wLength=testLength[imax]; istate=pixelArray[wPixel-1].GetState(); //Choose the TOF dead space // if(istate==0 && largestWeightTracksOutTof>wmax) { // if(istate==0 && largestWeightTracksOutTof>=sw) { if(istate==0 && sumWeightTracksOutTof>sGeomWeigth) { itestc=itest; itest=0; } } if(itest>0) { // Set for MyTrack: Pixel trackArray[itrack-1].SetPixel(wPixel); istate=pixelArray[wPixel-1].GetState(); if(istate) { // Set for MyTrack: Pixel, Length, TOF, MassTOF //fp //time=pixelArray[wPixel-1].GetTime(); time=pixelArray[wPixel-1].GetRealTime(); trackArray[itrack-1].SetLength(wLength); trackArray[itrack-1].SetTof(time); inverseOfParticleSpeed=time/wLength; //w=900.*inverseOfParticleSpeed*inverseOfParticleSpeed-1.; w=(100.*AliTOFConstants::fgkSpeedOfLight)*(100.*AliTOFConstants::fgkSpeedOfLight)*inverseOfParticleSpeed*inverseOfParticleSpeed-1.; w2=pvtx*pvtx; Float_t squareMass=w2*w; mass=TMath::Sqrt(TMath::Abs(squareMass)); if(w<0.) mass=-mass; trackArray[itrack-1].SetMassTOF(mass); // Set for MyTrack: Matching match=4; // if(ipart==pixelArray[wPixel-1].GetTrack()) match=3; if( (ipart==pixelArray[wPixel-1].GetTrack()) && hitArray[pixelArray[wPixel-1].GetHit()-1].GetNoise()==0)match=3; imatched=pixelArray[wPixel-1].GetTrackMatched(); // Set for TOFPixel the number of matched track pixelArray[wPixel-1].SetTrackMatched(itrack); if(imatched>0) { matchw=trackArray[imatched-1].GetMatching(); if(match==3 && matchw==4) t34++; if(match==3 && matchw==2) t32++; if(match==4 && matchw==4) t44++; if(match==4 && matchw==3) t43++; if(match==4 && matchw==2) t42++; if(iTOFpixel[ipart]==0 || iTOFpixel[trackArray[imatched-1].GetTrack()]==0) { m20++; } else if(iTOFpixel[ipart]==iTOFpixel[trackArray[imatched-1].GetTrack()]) { m22++; } else { m23++; wPix=iTOFpixel[ipart]; if(PRINT && iPrintM1==10 && iPrintM2<10) { if(iPrintM2==0) { printf("*** test print for tracks matched with the pixel for with we had matched track\n"); } iPrintM2++; printf(" m=2: ipart=%i, pdgCode=%i, p=%f, theta0=%f, %i Pixel(LP=%i,SP=%i,P=%i) \n", ipart,pdgCode,p,theta0,wPix, pixelArray[wPix-1].GetSector(),pixelArray[wPix-1].GetPlate(),pixelArray[wPix-1].GetPixel()); printf(" mat=%i, %i Pixel(LP=%i,SP=%i,P=%i), Test(n=%i,i=%i,w=%f,z=%f), wst=%i \n", match,wPixel, pixelArray[wPixel-1].GetSector(),pixelArray[wPixel-1].GetPlate(),pixelArray[wPixel-1].GetPixel(), itest,imax,wmax,testZ[imax],wstate); Int_t fstat,istat; for(Int_t i=0;i0) fstat=1; w=(fpadefficiency*fstat+(1.-fpadefficiency)*(1-fstat))*ntest[i]; if(istat>0) printf(" %i: %i Pixel(LP=%i,SP=%i,P=%i), istat=%i, ntest=%i, w=%f\n",i+1, wPix,pixelArray[wPix-1].GetSector(),pixelArray[wPix-1].GetPlate(),pixelArray[wPix-1].GetPixel(), istat,ntest[i],w); } printf(" mat=%i, %i Pixel \n",matchw,trackArray[imatched-1].GetPad()); } } if(wstate>1) m2state++; smat[matchw+4]--; match=2; trackArray[imatched-1].SetMatching(match); smat[match+4]++; } // if(imatched>0) } else { //else if(istate) match=1; if(iTOFpixel[ipart]==0) m10++; if(PRINT && iPrintM1<10) { Int_t wPix; wPix=iTOFpixel[ipart]; if(wPix) { if(iPrintM1==0) { printf("*** test print for tracks fired a pixel but matched with non-fired pixel\n"); } iPrintM1++; printf(" m=1: itra=%i,ipart=%i, pdgCode=%i, p=%f, theta0=%f, %i Pixel(LP=%i,SP=%i,P=%i) \n", itra,ipart,pdgCode,p,theta0,wPix, pixelArray[wPix-1].GetSector(),pixelArray[wPix-1].GetPlate(),pixelArray[wPix-1].GetPixel()); printf(" mat=%i, %i Pixel(LP=%i,SP=%i,P=%i), Test(n=%i,i=%i,w=%f,z=%f), wst=%i \n", match,wPixel, pixelArray[wPixel-1].GetSector(),pixelArray[wPixel-1].GetPlate(),pixelArray[wPixel-1].GetPixel(), itest,imax,wmax,testZ[imax],wstate); } } //end if(PRINT && iPrintM1<10) } //end if(istate) } else { match=-1-iMaxTestTracksOutTof; } //end itest trackArray[itrack-1].SetMatching(match); // if(iTestGmax==1) hMTT->Fill(match); smat[match+4]++; sumOfTheta /= iTestTrack; itest=itestc; //Test if(PRINT) { if(iTOFpixel[ipart] && match!=3) { particle = (TParticle*)gAlice->GetMCApp()->Particle(ipart); //for V3.05 printf(" ipixel=%i (Sector=%i, Plate=%i, Strip=%i, Pixel=%i), fired by %i track\n",iTOFpixel[ipart], pixelArray[iTOFpixel[ipart]-1].GetSector(),pixelArray[iTOFpixel[ipart]-1].GetPlate(), pixelArray[iTOFpixel[ipart]-1].GetStrip(),pixelArray[iTOFpixel[ipart]-1].GetPixel(), pixelArray[iTOFpixel[ipart]-1].GetTrack()); printf(" indexOfTestTrack=%i itest=%i weightTestTracksOutTof[4]=%f weightTestTracksOutTof[2]=%f weightTestTracksOutTof[1]=%f weightTestTracksOutTof[0]=%f\n", indexOfTestTrack,itest,weightTestTracksOutTof[3],weightTestTracksOutTof[2], weightTestTracksOutTof[1],weightTestTracksOutTof[0]); if(itest) { printf(" take ipixel=%i (Sector=%i, Plate=%i, Strip=%i, Pixel=%i), (fired by %i track), match=%i\n", wPixel,pixelArray[wPixel-1].GetSector(),pixelArray[wPixel-1].GetPlate(),pixelArray[wPixel-1].GetStrip(), pixelArray[wPixel-1].GetPixel(),pixelArray[wPixel-1].GetTrack(),match); } } } if(PRINT && itra<10 ) { if(itest) { cout << " number of pixels with test tracks=" << itest << endl; for(Int_t i=0;iRmin \n",idelR,w); } w=iRmin1; if(idelR1!=0) { w /= idelR1; printf(" %i tracks with delR1 (|z|<175), %f of them have delR>Rmin \n",idelR1,w); } w=iRmin2; if(idelR2!=0) { w /= idelR2; printf(" %i tracks with delR2 (|z|>175), %f of them have delR>Rmin \n",idelR2,w); } cout << " ******************** End of matching **********" << endl; delete [] ntest; delete [] testPixel; delete [] testLength; delete [] testRho; delete [] testZ; delete [] testWeight; } //____________________________________________________________________________ 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){ // itrack : total number of TPC selected tracks // for the caller is ntotTPCtracks cout << " ******************** Start of searching non-matched fired pixels **********" << endl; const Int_t charge[48]={ 0, 1,-1, 0, 1,-1, 0, 1,-1, 0, 1,-1, 0, 1,-1, 0, 0, 0, 1, 0, -1, 0,-1,-1, 0, 0,-1, 0, 1, 0, 1, 1, 0, 0, 1,-1, 0, 0, 1,-1, 1, 1,-1, 0, 1, 1, 2, 0}; Int_t macthm1=0; Int_t macthm2=0; Int_t macthm3=0; Int_t macthm4=0; Int_t macth0=0; Int_t macth1=0; Int_t macth2=0; Int_t macth3=0; Int_t macth4=0; Float_t smat[9],smat0[9],smat1[9]; for(Int_t i=0;i<9;i++) { smat[i]=0.; smat0[i]=0.; smat1[i]=0.; } Int_t nFiredPixelsNotMatchedWithTracks=0; Int_t istate; for (Int_t i=0; i0) { if(hitArray[pixelArray[i].GetHit()-1].GetNoise()==0) { if(iparticle[pixelArray[i].GetTrack()]==0) nTPCHitMissing++; } } } printf(" %i pixels fired by track hit without a hit on the last layer of TPC\n",nTPCHitMissing); Int_t icharge=0; // total number of charged particles Int_t iprim=0; // number of primaries Int_t ipions=0; // number of primary pions Int_t ikaons=0; // number of primary kaons Int_t iprotons=0; // number of primary protons Int_t ielectrons=0;// number of primary electrons Int_t imuons=0; // number of primary muons Float_t particleTypeArray[6][5][2]; for (Int_t index1=0;index1<6;index1++) { for (Int_t index2=0;index2<5;index2++) { for (Int_t index3=0;index3<2;index3++) { particleTypeArray[index1][index2][index3]=0.; } } } Int_t nTOFhitsWithNoTPCTracks=0; // to be moved later when used /* TObjArray *Particles = gAlice->Particles(); Int_t numberOfParticles=Particles->GetEntries(); cout << "numberOfParticles " << numberOfParticles << endl; // fpdbg if(numberOfParticles>fMaxAllTracks) numberOfParticles=fMaxAllTracks; */ for (Int_t i=0; iGetMCApp()->Particle(i); if(charge[PDGtoGeantCode(part->GetPdgCode())-1]) { icharge++; /* cout << "charged particles " << icharge << endl; */ Int_t particleType=0; Int_t absPdgCode = TMath::Abs(part->GetPdgCode()); switch (absPdgCode) { case 211: particleType=3; break ; case 321: particleType=2; break ; case 2212: particleType=1; break ; case 11: particleType=4; break ; case 13: particleType=5; break ; } if(part->GetFirstMother() < 0) { iprim++; switch (particleType) { case 1: iprotons++; break ; case 2: ikaons++; break ; case 3: ipions++; break ; case 4: ielectrons++; break ; case 5: imuons++; break ; } } Int_t match=0; Float_t wLength=-1.; Float_t time=-1.; Float_t mass=-1.; Int_t itr=iparticle[i]; // get the track number for the current charged particle if(iTOFpixel[i]>0 && itr==0) nTOFhitsWithNoTPCTracks++; if(itr) { match=trackArray[itr-1].GetMatching(); //cout << "match " << match << endl; wLength=trackArray[itr-1].GetLength(); //cout << "wLength " << wLength << endl; time=trackArray[itr-1].GetTof(); mass=trackArray[itr-1].GetMassTOF(); //cout << "mext " << mass << endl; // if(PRINT && (i>789 && i<800) ) cout << i << " track: l=" << wLength << " TOF=" << time << " m=" << mass << " match=" << match << endl; if(iTOFpixel[i]==0) { smat0[match+4]++; wLength=-wLength; } } Int_t ikparen=part->GetFirstMother(); Int_t imam; if(ikparen<0) { imam=0; } else { imam=part->GetPdgCode(); } Int_t evnumber=gAlice->GetEvNumber(); if(match==-1) macthm1++; if(match==-2) macthm2++; if(match==-3) macthm3++; if(match==-4) macthm4++; if(match==0) macth0++; if(match==1) macth1++; if(match==2) macth2++; if(match==3) macth3++; if(match==4) macth4++; foutputntuple->Fill(evnumber,part->GetPdgCode(),imam,part->Vx(),part->Vy(),part->Vz(),part->Px(),part->Py(),part->Pz(),toftime[i],wLength,match,time,mass); // ----------------------------------------------------------- // Filling 2 dimensional Histograms true time vs matched time // Filling 1 dimensional Histogram true time - matched time // // time = time associated to the matched pad [ns] // it could be the average time of the cluster fired // // toftime[i] = real time (including pulse height delays) [s] // // // if (time>=0) { // if (imam==0) TimeTrueMatched->Fill(time, toftime[i]*1E+09); // if (imam==0) DeltaTrueTimeMatched->Fill(time-toftime[i]*1E+09); // } // //--------------------------------------------------------------- if(match==-4 || match>0) { Int_t matchW; matchW=match; if(match==-4) matchW=1; if(particleType) { particleTypeArray[particleType-1][matchW-1][1]++; particleTypeArray[5][matchW-1][1]++; particleTypeArray[particleType-1][4][1]++; particleTypeArray[5][4][1]++; if(part->GetFirstMother() < 0) { particleTypeArray[particleType-1][matchW-1][0]++; particleTypeArray[5][matchW-1][0]++; particleTypeArray[particleType-1][4][0]++; particleTypeArray[5][4][0]++; // fill histos for QA //if(particleType==3 && matchW==3) hPiWithTrueTime->Fill(sqrt((part->Px())*(part->Px())+(part->Py())*(part->Py())+(part->Pz())*(part->Pz()))); //if(particleType==2 && matchW==3) hKWithTrueTime->Fill(sqrt((part->Px())*(part->Px())+(part->Py())*(part->Py())+(part->Pz())*(part->Pz()))); //if(particleType==1 && matchW==3) hPWithTrueTime->Fill(sqrt((part->Px())*(part->Px())+(part->Py())*(part->Py())+(part->Pz())*(part->Pz()))); // } // close if(part->GetFirstMother() < 0) } // close if(particleType) } // close if(match==-4 || match>0) } // close if(charge[PDGtoGeantCode(part->GetPdgCode())-1]) } // close for (Int_t i=0; i4) imatch += (Int_t) (smat[i]); // cout << " " << smat[i]*100./itrack << " % of them (="<Fill(Acceptance[i]); // printf(" species: %i %f\n",i+1,Acceptance[i]); } // cout << " TOF acceptance for the primary species: 1-p, 2-K, 3-pi, 4-e, 5-mu, 6-all" << endl; for(Int_t i=0; i<6; i++) { Acceptance[i]=0.; if(PrimaryInBarrel[i]) Acceptance[i]=100.*PrimaryReachedTOF[i]/PrimaryInBarrel[i]; //hTOFacceptance[i]->Fill(Acceptance[i]); // printf(" species: %i %f\n",i+1,Acceptance[i]); } for (Int_t index1=0;index1<6;index1++) { for (Int_t index2=0;index2<4;index2++) { for (Int_t index3=0;index3<2;index3++) { if(particleTypeArray[index1][4][index3]) particleTypeArray[index1][index2][index3]= 100.*particleTypeArray[index1][index2][index3]/particleTypeArray[index1][4][index3]; } } } cout << "species: 1-p, 2-K, 3-pi, 4-e, 5-mu, 6-all" << endl; cout << " matched pixels(%): 1-unfired 2-double 3-true 4-wrong 5-total number of tracks" << endl; cout << " primary tracks:" << endl; for (Int_t i=0;i<6;i++) { cout << i+1 << " " << particleTypeArray[i][0][0] << " " << particleTypeArray[i][1][0] << " " << particleTypeArray[i][2][0] << " " << particleTypeArray[i][3][0] << " " << particleTypeArray[i][4][0] << endl; } // cout<<" contam.for all prim.(%)="<<100*particleTypeArray[5][3][0]/(particleTypeArray[5][3][0]+particleTypeArray[5][2][0])<