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