JETAN/AliJetKineReader.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 //
  17 // Jet kinematics reader
  18 // MC reader for jet analysis
  19 // Author: Andreas Morsch
  20 // andreas.morsch@cern.ch
  21 //
  22
  23 // From root ...
  24 #include <TClonesArray.h>
  25 #include <TPDGCode.h>
  26 #include <TParticle.h>
  27 #include <TParticlePDG.h>
  28 #include <TVector3.h>
  29 #include <TLorentzVector.h>
  30 #include <TSystem.h>
  31 #include <TDatabasePDG.h>
  32 #include <TRandom.h>
  33 // From AliRoot ...
  34 #include "AliJetKineReaderHeader.h"
  35 #include "AliJetKineReader.h"
  36 #include "AliRunLoader.h"
  37 #include "AliStack.h"
  38
  39 ClassImp(AliJetKineReader)
  40
  41 AliJetKineReader::AliJetKineReader()
  42 {
  43   // Constructor
  44   fReaderHeader = 0x0;
  45   fRunLoader    = 0x0;
  46   fMass         = 0;
  47   fPdgC         = 0;
  48 }
  49
  50 //____________________________________________________________________________
  51
  52 AliJetKineReader::~AliJetKineReader()
  53 {
  54   // Destructor
  55   delete fReaderHeader;
  56 }
  57
  58 //____________________________________________________________________________
  59
  60 void AliJetKineReader::OpenInputFiles()
  61 {
  62     // Opens the input file using the run loader
  63     const char* dirName = fReaderHeader->GetDirectory();
  64     char path[256];
  65     sprintf(path, "%s/galice.root",dirName);
  66     fRunLoader = AliRunLoader::Open(path);
  67     fRunLoader->LoadKinematics();
  68     fRunLoader->LoadHeader();
  69
  70     Int_t nMax = fRunLoader->GetNumberOfEvents();
  71     printf("\nTotal number of events = %d", nMax);
  72
  73   // set number of events in header
  74     if (fReaderHeader->GetLastEvent() == -1)
  75         fReaderHeader->SetLastEvent(nMax);
  76     else {
  77         Int_t nUsr = fReaderHeader->GetLastEvent();
  78         fReaderHeader->SetLastEvent(TMath::Min(nMax, nUsr));
  79     }
  80 }
  81
  82 //____________________________________________________________________________
  83
  84 void AliJetKineReader::FillMomentumArray(Int_t event)
  85 {
  86 //
  87 // Fill momentum array for event
  88 //
  89     Int_t goodTrack = 0;
  90     // Clear array
  91     ClearArray();
  92     // Get event from runloader
  93     fRunLoader->GetEvent(event);
  94     // Get the stack
  95     AliStack* stack = fRunLoader->Stack();
  96     // Number of primaries
  97     Int_t nt = stack->GetNprimary();
  98     // Get cuts set by user and header
  99     Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
 100     fAliHeader = fRunLoader->GetHeader();
 101
 102     // Loop over particles
 103     Int_t* flag  = new Int_t[nt];
 104     for (Int_t it = 0; it < nt; it++) {
 105         TParticle *part = stack->Particle(it);
 106         Int_t   status  = part->GetStatusCode();
 107         Int_t   pdg     = TMath::Abs(part->GetPdgCode());
 108         Float_t pt      = part->Pt();
 109
 110         // Skip non-final state particles, neutrinos and particles with pt < pt_min
 111
 112         if (
 113             (status != 1)
 114             || (pdg == 12 || pdg == 14 || pdg == 16)
 115             || (pt < ptMin)
 116             ) continue;
 117
 118
 119         Float_t p       = part->P();
 120         Float_t p0      = p;
 121         Float_t eta     = part->Eta();
 122         Float_t phi     = part->Phi();
 123
 124         // Fast simulation of TPC if requested
 125         if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
 126             // Charged particles only
 127             Float_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
 128             if (charge == 0)               continue;
 129             // Simulate efficiency
 130             if (!Efficiency(p0, eta, phi)) continue;
 131             // Simulate resolution
 132             p = SmearMomentum(4, p0);
 133         } // Fast TPC
 134
 135         // Fast simulation of EMCAL if requested
 136         if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
 137             Float_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
 138             // Charged particles
 139             if (charge != 0) {
 140                 // Simulate efficiency
 141                 if (!Efficiency(p0, eta, phi)) continue;
 142                 // Simulate resolution
 143                 p = SmearMomentum(4, p0);
 144             } // charged
 145             // Neutral particles
 146             // Exclude K0L, n, nbar
 147             if (pdg == kNeutron || pdg == kK0Long) continue;
 148         } // Fast EMCAL
 149
 150         fMass = part->GetCalcMass();
 151         fPdgC = part->GetPdgCode();
 152         // Fill momentum array
 153         Float_t r  = p/p0;
 154 //              printf("smeared %13.3f %13.3f %13.3f\n", p0, p, r);
 155
 156         Float_t px = r * part->Px();
 157         Float_t py = r * part->Py();
 158         Float_t pz = r * part->Pz();
 159         Float_t m  = part->GetMass();
 160         Float_t e  = TMath::Sqrt(px * px + py * py + pz * pz + m * m);
 161
 162         new ((*fMomentumArray)[goodTrack])
 163             TLorentzVector(px, py, pz, e);
 164         goodTrack++;
 165   }
 166     fSignalFlag.Set(goodTrack,flag);
 167     printf("\nNumber of good tracks in event %d= %d \n",event,goodTrack);
 168 }
 169
 170
 171 Float_t AliJetKineReader::SmearMomentum(Int_t ind, Float_t p)
 172 {
 173 //
 174 //  The relative momentum error, i.e.
 175 //  (delta p)/p = sqrt (a**2 + (b*p)**2) * 10**-2,
 176 //  where typically a = 0.75 and b = 0.16 - 0.24 depending on multiplicity
 177 //  (the lower value is for dn/d(eta) about 2000, and the higher one for 8000)
 178 //
 179 //  If we include information from TRD b will be by a factor 2/3 smaller.
 180 //
 181 //  ind = 1: high multiplicity
 182 //  ind = 2: low  multiplicity
 183 //  ind = 3: high multiplicity + TRD
 184 //  ind = 4: low  multiplicity + TRD
 185
 186     Float_t pSmeared;
 187     Float_t a = 0.75;
 188     Float_t b = 0.12;
 189
 190     if (ind == 1) b = 0.12;
 191     if (ind == 2) b = 0.08;
 192     if (ind == 3) b = 0.12;
 193     if (ind == 4) b = 0.08;
 194
 195     Float_t sigma =  p * TMath::Sqrt(a * a + b * b * p * p)*0.01;
 196     pSmeared = p + gRandom->Gaus(0., sigma);
 197     return pSmeared;
 198 }
 199 Bool_t AliJetKineReader::Efficiency(Float_t p, Float_t /*eta*/, Float_t phi)
 200 {
 201 //
 202 // Fast simulation of geometrical acceptance and tracking efficiency
 203 //
 204 //  Tracking
 205     Float_t eff = 0.99;
 206     if (p < 0.5) eff -= (0.5-p)*0.2/0.3;
 207 // Geometry
 208     if (p > 0.5) {
 209         phi *= 180. / TMath::Pi();
 210         // Sector number 0 - 17
 211         Int_t isec  = Int_t(phi / 20.);
 212         // Sector centre
 213         Float_t phi0 = isec * 20. + 10.;
 214         Float_t phir = TMath::Abs(phi-phi0);
 215         // 2 deg of dead space
 216         if (phir > 9.) eff = 0.;
 217     }
 218     if (gRandom->Rndm() > eff) {
 219         return kFALSE;
 220     } else {
 221         return kTRUE;
 222     }
 223 }
 224