Temporary removing AliEMCALv3 and AliEMCALHitv1 from the compilation until they are...
[u/mrichter/AliRoot.git] / JETAN / AliJetKineReader.cxx
CommitLineData
99e5fe42 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
39ClassImp(AliJetKineReader)
40
41AliJetKineReader::AliJetKineReader()
42{
43 // Constructor
44 fReaderHeader = 0x0;
45 fRunLoader = 0x0;
46 fMass = 0;
47 fPdgC = 0;
48}
49
50//____________________________________________________________________________
51
52AliJetKineReader::~AliJetKineReader()
53{
54 // Destructor
55 delete fReaderHeader;
56}
57
58//____________________________________________________________________________
59
60void 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
84void AliJetKineReader::FillMomentumArray(Int_t event)
85{
86//
87// Fill momentum array for event
88//
89 Int_t goodTrack = 0;
90 // Clear array
91 ClearArray();
dd1aa8ee 92 // Get event from runloader
99e5fe42 93 fRunLoader->GetEvent(event);
dd1aa8ee 94 // Get the stack
99e5fe42 95 AliStack* stack = fRunLoader->Stack();
dd1aa8ee 96 // Number of primaries
99e5fe42 97 Int_t nt = stack->GetNprimary();
dd1aa8ee 98 // Get cuts set by user and header
99e5fe42 99 Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
100 fAliHeader = fRunLoader->GetHeader();
101
102 // Loop over particles
48bb52d8 103 Int_t* flag = new Int_t[nt];
99e5fe42 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
dd1aa8ee 110 // Skip non-final state particles, neutrinos and particles with pt < pt_min
111
99e5fe42 112 if (
113 (status != 1)
dd1aa8ee 114 || (pdg == 12 || pdg == 14 || pdg == 16)
99e5fe42 115 || (pt < ptMin)
116 ) continue;
117
dd1aa8ee 118
48bb52d8 119 Float_t p = part->P();
dd1aa8ee 120 Float_t p0 = p;
48bb52d8 121 Float_t eta = part->Eta();
122 Float_t phi = part->Phi();
123
dd1aa8ee 124 // Fast simulation of TPC if requested
99e5fe42 125 if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
48bb52d8 126 // Charged particles only
99e5fe42 127 Float_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
dd1aa8ee 128 if (charge == 0) continue;
48bb52d8 129 // Simulate efficiency
dd1aa8ee 130 if (!Efficiency(p0, eta, phi)) continue;
48bb52d8 131 // Simulate resolution
dd1aa8ee 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
99e5fe42 149
150 fMass = part->GetCalcMass();
151 fPdgC = part->GetPdgCode();
152 // Fill momentum array
dd1aa8ee 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);
99e5fe42 161
162 new ((*fMomentumArray)[goodTrack])
dd1aa8ee 163 TLorentzVector(px, py, pz, e);
99e5fe42 164 goodTrack++;
165 }
48bb52d8 166 fSignalFlag.Set(goodTrack,flag);
dd1aa8ee 167 printf("\nNumber of good tracks in event %d= %d \n",event,goodTrack);
99e5fe42 168}
169
170
48bb52d8 171Float_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
dd1aa8ee 195 Float_t sigma = p * TMath::Sqrt(a * a + b * b * p * p)*0.01;
48bb52d8 196 pSmeared = p + gRandom->Gaus(0., sigma);
197 return pSmeared;
198}
199Bool_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