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