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