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