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