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