Coding violations corrected.
[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
83a444b1 16//-------------------------------------------------------------------------
99e5fe42 17// Jet kinematics reader
18// MC reader for jet analysis
83a444b1 19// Author: Andreas Morsch (andreas.morsch@cern.ch)
20//-------------------------------------------------------------------------
99e5fe42 21
22// From root ...
23#include <TClonesArray.h>
24#include <TPDGCode.h>
25#include <TParticle.h>
26#include <TParticlePDG.h>
99e5fe42 27#include <TLorentzVector.h>
28#include <TSystem.h>
29#include <TDatabasePDG.h>
30#include <TRandom.h>
31// From AliRoot ...
85575952 32#include "AliAODJet.h"
d625f15e 33#include "AliPDG.h"
99e5fe42 34#include "AliJetKineReaderHeader.h"
35#include "AliJetKineReader.h"
85575952 36#include "AliMCEventHandler.h"
7f7091b8 37#include "AliMCEvent.h"
99e5fe42 38#include "AliStack.h"
0ffa8579 39#include "AliHeader.h"
40#include "AliGenPythiaEventHeader.h"
85575952 41#include "AliLog.h"
99e5fe42 42
43ClassImp(AliJetKineReader)
44
1b7d5d7e 45AliJetKineReader::AliJetKineReader():
b45b0c92 46 AliJetReader(),
85575952 47 fAliHeader(0),
b0394c65 48 fMCEvent(0),
85575952 49 fGenJets(0)
99e5fe42 50{
1b7d5d7e 51 // Default constructor
99e5fe42 52}
53
54//____________________________________________________________________________
55
56AliJetKineReader::~AliJetKineReader()
57{
58 // Destructor
dee50a65 59 if (fAliHeader) {
60 delete fAliHeader;
61 fAliHeader = 0;
62 }
99e5fe42 63}
64
65//____________________________________________________________________________
66
67void AliJetKineReader::OpenInputFiles()
68{
83a444b1 69 // Opens the input file using the run loader
85575952 70 // Not used anymore
99e5fe42 71}
72
73//____________________________________________________________________________
74
ae24a5a1 75Bool_t AliJetKineReader::FillMomentumArray()
99e5fe42 76{
83a444b1 77 // Fill momentum array for event
85575952 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];
83a444b1 83
85575952 84 ClearArray();
85 // Get the stack
b0394c65 86 AliStack* stack = fMCEvent->Stack();
85575952 87 // Number of primaries
88 Int_t nt = stack->GetNprimary();
4faf4c3d 89
85575952 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();
b0394c65 94 fAliHeader = fMCEvent->Header();
4faf4c3d 95
85575952 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();
4faf4c3d 104
85575952 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();
eadf8168 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
85575952 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
4faf4c3d 135
136 // Fast simulation of EMCAL if requested
85575952 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
4faf4c3d 147 // Neutral particles (exclude K0L, n, nbar)
85575952 148 if (pdg == kNeutron || pdg == kK0Long) continue;
149 } // End fast EMCAL
150
4faf4c3d 151 // Fill momentum array
85575952 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();
d625f15e 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 }
85575952 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);
b0394c65 188 AliInfo(Form(" Number of good tracks %d \n", goodTrack));
85575952 189
190 delete[] sflag;
191 delete[] cflag;
192
193 return kTRUE;
99e5fe42 194}
195
196
48bb52d8 197Float_t AliJetKineReader::SmearMomentum(Int_t ind, Float_t p)
198{
83a444b1 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;
48bb52d8 223}
83a444b1 224
225
48bb52d8 226Bool_t AliJetKineReader::Efficiency(Float_t p, Float_t /*eta*/, Float_t phi)
227{
83a444b1 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
48bb52d8 251}
252
85575952 253TClonesArray* AliJetKineReader::GetGeneratedJets()
0ffa8579 254{
85575952 255 //
0ffa8579 256 // Get generated jets from mc header
85575952 257 //
258 if (!fGenJets) {
259 fGenJets = new TClonesArray("AliAODJets", 100);
260 } else {
261 fGenJets->Clear();
262 }
263
264
0ffa8579 265 AliHeader* alih = GetAliHeader();
85575952 266 if (alih == 0) return 0;
0ffa8579 267 AliGenEventHeader * genh = alih->GenEventHeader();
85575952 268 if (genh == 0) return 0;
269
0ffa8579 270 Int_t nj =((AliGenPythiaEventHeader*)genh)->NTriggerJets();
85575952 271 for (Int_t i = 0; i < nj; i++) {
0ffa8579 272 Float_t p[4];
273 ((AliGenPythiaEventHeader*)genh)->TriggerJet(i,p);
85575952 274 new ((*fGenJets)[i]) AliAODJet(p[0], p[1], p[2], p[3]);
0ffa8579 275 }
85575952 276
277 return fGenJets;
278}
279
c4688b3a 280void AliJetKineReader::SetInputEvent(const TObject* /*esd*/, const TObject* /*aod*/, const TObject* mc)
85575952 281{
282 // Connect the input event
b0394c65 283 fMCEvent = (AliMCEvent*) mc;
0ffa8579 284}