JETAN/AliJetFillCalTrkTrackMC.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //--------------------------------------------------
  19 // Filling of CalTrkTrack objects in the MC reader task
  20 //
  21 // Author: magali.estienne@subatech.in2p3.fr
  22 //         alexandre.shabetai@cern.ch
  23 //-------------------------------------------------
  24
  25 // --- ROOT system ---
  26 #include <TDatabasePDG.h>
  27 #include <TRandom.h>
  28 #include <TPDGCode.h>
  29
  30 // --- AliRoot system ---
  31 #include "AliJetFillCalTrkTrackMC.h"
  32 #include "AliAODMCParticle.h"
  33 #include "AliJetKineReaderHeader.h"
  34 #include "AliVEvent.h"
  35 #include "AliMCEvent.h"
  36
  37 using std::cout;
  38 using std::endl;
  39 ClassImp(AliJetFillCalTrkTrackMC)
  40
  41 /////////////////////////////////////////////////////////////////////////
  42
  43 AliJetFillCalTrkTrackMC::AliJetFillCalTrkTrackMC():
  44   AliJetFillCalTrkEvent(),
  45   fHadCorr(0),
  46   fApplyMIPCorrection(kTRUE),
  47   fVEvt(0x0),
  48   fMCEvent(0x0)
  49 {
  50   // constructor
  51 }
  52
  53 //-----------------------------------------------------------------------
  54 AliJetFillCalTrkTrackMC::AliJetFillCalTrkTrackMC(AliVEvent* evt):
  55   AliJetFillCalTrkEvent(),
  56   fHadCorr(0),
  57   fApplyMIPCorrection(kTRUE),
  58   fVEvt(evt),
  59   fMCEvent(0x0)
  60 {
  61   // constructor
  62 }
  63
  64 //-----------------------------------------------------------------------
  65 AliJetFillCalTrkTrackMC::AliJetFillCalTrkTrackMC(const AliJetFillCalTrkTrackMC &det):
  66   AliJetFillCalTrkEvent(det),
  67   fHadCorr(det.fHadCorr),
  68   fApplyMIPCorrection(det.fApplyMIPCorrection),
  69   fVEvt(det.fVEvt),
  70   fMCEvent(det.fMCEvent)
  71 {
  72   // Copy constructor
  73 }
  74
  75 //-----------------------------------------------------------------------
  76 AliJetFillCalTrkTrackMC& AliJetFillCalTrkTrackMC::operator=(const AliJetFillCalTrkTrackMC& other)
  77 {
  78   // Assignment
  79   if (this != &other) {
  80    fHadCorr = other.fHadCorr;
  81    fApplyMIPCorrection = other.fApplyMIPCorrection;
  82    fVEvt = other.fVEvt;
  83    fMCEvent = other.fMCEvent;
  84   }
  85
  86   return (*this);
  87
  88 }
  89
  90 //-----------------------------------------------------------------------
  91 AliJetFillCalTrkTrackMC::~AliJetFillCalTrkTrackMC()
  92 {
  93   // destructor
  94 }
  95
  96 //-----------------------------------------------------------------------
  97 void AliJetFillCalTrkTrackMC::Exec(Option_t const * /*option*/)
  98 {
  99   // Main method.
 100   // Fill AliJetFillCalTrkTrackMC the with the charged particle information
 101
 102   fDebug = fReaderHeader->GetDebug();
 103   fOpt = fReaderHeader->GetDetector();
 104   TString datatype = fReaderHeader->GetDataType();
 105   datatype.ToUpper();
 106
 107   Int_t type=0;
 108
 109   if      (datatype.Contains("AODMC2B")) { type=kTrackAODMCChargedAcceptance; }
 110   else if (datatype.Contains("AODMC2"))  { type=kTrackAODMCCharged; }
 111   else if (datatype.Contains("AODMC"))   { type=kTrackAODMCAll; }
 112
 113   else if (!datatype.Contains("AOD") && datatype.Contains("MC2")) {type=kTrackKineCharged;}
 114   else if (!datatype.Contains("AOD") && datatype.Contains("MC"))  {type=kTrackKineAll;}
 115
 116   else { cout<< "Unknown Data type !" << endl; }
 117
 118   // temporary storage of signal and pt cut flag
 119   Bool_t sflag = 0;
 120   Bool_t cflag = 0;
 121
 122   // get cuts set by user
 123   Float_t ptMin  = fReaderHeader->GetPtCut();
 124   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
 125   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
 126   Float_t phiMin = fReaderHeader->GetFiducialPhiMin();
 127   Float_t phiMax = fReaderHeader->GetFiducialPhiMax();
 128
 129   if(fDebug>2) Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
 130
 131   Int_t goodTrackStd = 0;
 132   Int_t goodTrackNonStd = 0;
 133
 134   if (type == kTrackKineAll || type == kTrackKineCharged){
 135
 136     FillKine();
 137   } // Kine
 138
 139   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
 140     if(!fVEvt) return;
 141     TClonesArray *mcArray = dynamic_cast<TClonesArray*>(fVEvt->FindListObject(AliAODMCParticle::StdBranchName()));
 142     if(!mcArray) {
 143       Printf("%s:%d No MC particle branch found",(char*)__FILE__,__LINE__);
 144       return;
 145     }
 146     for(int it = 0;it < mcArray->GetEntriesFast();++it){
 147       cflag=sflag=0;
 148       AliAODMCParticle *part = (AliAODMCParticle*)(mcArray->At(it));
 149       if(!part->IsPhysicalPrimary())continue;
 150       Int_t   pdg     = TMath::Abs(part->GetPdgCode());
 151       // exclude neutrinos anyway
 152       if((pdg == 12 || pdg == 14 || pdg == 16)) continue;
 153       cflag = ( part->Pt()>ptMin ) ? 1 : 0;
 154       sflag = ( TMath::Abs(part->GetLabel()) < 10000 ) ? 1 : 0;
 155       if(type == kTrackAODMCAll){
 156         fCalTrkEvent->AddCalTrkTrack(part,cflag,sflag);
 157         goodTrackStd++;
 158       }
 159       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
 160         if(part->Charge()==0)continue;
 161         if(kTrackAODMCCharged){
 162           fCalTrkEvent->AddCalTrkTrack(part,cflag,sflag);
 163           goodTrackStd++;
 164         }
 165         else {
 166           if((part->Eta()>etaMax) || (part->Eta()<etaMin)) continue;
 167           Float_t phi = ( (part->Phi()) < 0) ? (part->Phi()) + 2. * TMath::Pi() : (part->Phi());
 168           if((phi>phiMax) || (phi<phiMin)) continue;
 169           fCalTrkEvent->AddCalTrkTrack(part,cflag,sflag);
 170           goodTrackStd++;
 171         }
 172       }
 173     }
 174    if(fDebug>0) printf("Number of good tracks selected: %5d \n", goodTrackStd+goodTrackNonStd);
 175   } // AODMCparticle
 176
 177 }
 178
 179 //-----------------------------------------------------------------------
 180 Bool_t AliJetFillCalTrkTrackMC::FillKine()
 181 {
 182   // Fill event
 183   Int_t goodTrack = 0;
 184
 185   // Get the stack
 186   if(!fMCEvent) {cout<<"could not get MCEvent!"<<endl; return kFALSE;}
 187
 188   Int_t nt = fMCEvent->GetNumberOfTracks();
 189
 190   // Get cuts set by user and header
 191   Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
 192   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
 193   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
 194
 195   TLorentzVector p4;
 196   // Loop over particles
 197   for (Int_t it = 0; it < nt; it++) {
 198     if(!fMCEvent->IsPhysicalPrimary(it)) continue;
 199     AliVParticle* part = fMCEvent->GetTrack(it);
 200     Int_t   pdg    = TMath::Abs(part->PdgCode());
 201     Float_t pt     = part->Pt();
 202
 203     if( (pdg == 12 || pdg == 14 || pdg == 16)) continue;
 204
 205     Float_t p      = part->P();
 206     Float_t p0     = p;
 207     Float_t eta    = part->Eta();
 208     Float_t phi    = part->Phi();
 209     Float_t charge = part->Charge();
 210
 211     if (((AliJetKineReaderHeader*)fReaderHeader)->ChargedOnly()) {
 212       if (charge == 0) continue;
 213     } // End charged only
 214
 215     // Fast simulation of EMCAL if requested
 216     if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
 217       // Charged particles only
 218       if (charge != 0){
 219         // Simulate efficiency
 220         if (!Efficiency(p0, eta, phi)) continue;
 221         // Simulate resolution
 222         p = SmearMomentum(4, p0);
 223       } // end "if" charged particles
 224       // Neutral particles (exclude K0L, n, nbar)
 225       if (pdg == kNeutron || pdg == kK0Long) continue;
 226     } // End fast EMCAL
 227
 228     // Fast simulation of TPC if requested
 229     if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
 230       // Charged particles only
 231       if (charge == 0)               continue;
 232       // Simulate efficiency
 233       if (!Efficiency(p0, eta, phi)) continue;
 234       // Simulate resolution
 235       p = SmearMomentum(4, p0);
 236     } // End fast TPC
 237
 238     // Fill momentum array
 239     Float_t r  = p/p0;
 240     Float_t px = r * part->Px();
 241     Float_t py = r * part->Py();
 242     Float_t pz = r * part->Pz();
 243     Float_t m =  part->M();
 244
 245     Float_t e  = TMath::Sqrt(px * px + py * py + pz * pz + m * m);
 246     p4 = TLorentzVector(px, py, pz, e);
 247     if ((p4.Eta()>etaMax) || (p4.Eta()<etaMin)) continue;
 248
 249     // Flag used to store the r factor
 250     Float_t ptReso = r;
 251     Bool_t cflag = kFALSE, sflag = kTRUE;
 252     if (pt > ptMin) cflag = kTRUE; // track surviving pt cut
 253
 254     fCalTrkEvent->AddCalTrkTrackKine(part,cflag,sflag,ptReso);
 255
 256     goodTrack++;
 257
 258   } // track loop
 259
 260   if(fDebug>0) printf(" Number of good tracks %d \n", goodTrack);
 261
 262   return kTRUE;
 263
 264 }
 265
 266 //-----------------------------------------------------------------------
 267 Float_t AliJetFillCalTrkTrackMC::SmearMomentum(Int_t ind, Float_t p)
 268 {
 269   //  The relative momentum error, i.e.
 270   //  (delta p)/p = sqrt (a**2 + (b*p)**2) * 10**-2,
 271   //  where typically a = 0.75 and b = 0.16 - 0.24 depending on multiplicity
 272   //  (the lower value is for dn/d(eta) about 2000, and the higher one for 8000)
 273   //
 274   //  If we include information from TRD, b will be a factor 2/3 smaller.
 275   //
 276   //  ind = 1: high multiplicity
 277   //  ind = 2: low  multiplicity
 278   //  ind = 3: high multiplicity + TRD
 279   //  ind = 4: low  multiplicity + TRD
 280
 281   Float_t pSmeared;
 282   Float_t a = 0.75;
 283   Float_t b = 0.12;
 284
 285   if (ind == 1) b = 0.12;
 286   if (ind == 2) b = 0.08;
 287   if (ind == 3) b = 0.12;
 288   if (ind == 4) b = 0.08;
 289
 290   Float_t sigma =  p * TMath::Sqrt(a * a + b * b * p * p)*0.01;
 291   pSmeared = p + gRandom->Gaus(0., sigma);
 292   return pSmeared;
 293
 294 }
 295
 296 //-----------------------------------------------------------------------
 297 Bool_t AliJetFillCalTrkTrackMC::Efficiency(Float_t p, Float_t /*eta*/, Float_t phi)
 298 {
 299   // Fast simulation of geometrical acceptance and tracking efficiency
 300
 301   //  Tracking
 302   Float_t eff = 0.99;
 303   if (p < 0.5) eff -= (0.5-p)*0.2/0.3;
 304   // Geometry
 305   if (p > 0.5) {
 306     phi *= 180. / TMath::Pi();
 307     // Sector number 0 - 17
 308     Int_t isec  = Int_t(phi / 20.);
 309     // Sector centre
 310     Float_t phi0 = isec * 20. + 10.;
 311     Float_t phir = TMath::Abs(phi-phi0);
 312     // 2 deg of dead space
 313     if (phir > 9.) eff = 0.;
 314   }
 315
 316   if (gRandom->Rndm() > eff) {
 317     return kFALSE;
 318   } else {
 319     return kTRUE;
 320   }
 321
 322 }
 323
 324