]>
Commit | Line | Data |
---|---|---|
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 | |
40 | ClassImp(AliJetKineReader) | |
41 | ||
1b7d5d7e | 42 | AliJetKineReader::AliJetKineReader(): |
b45b0c92 | 43 | AliJetReader(), |
44 | fRunLoader(0x0), | |
45 | fAliHeader(0x0) | |
99e5fe42 | 46 | { |
1b7d5d7e | 47 | // Default constructor |
99e5fe42 | 48 | } |
49 | ||
50 | //____________________________________________________________________________ | |
51 | ||
52 | AliJetKineReader::~AliJetKineReader() | |
53 | { | |
54 | // Destructor | |
0ffa8579 | 55 | delete fAliHeader; |
99e5fe42 | 56 | } |
57 | ||
58 | //____________________________________________________________________________ | |
59 | ||
60 | void 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 | 85 | Bool_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 | |
201 | fSignalFlag.Set(goodTrack,sflag); | |
202 | fCutFlag.Set(goodTrack,cflag); | |
203 | printf("\nIn event %d, number of good tracks %d \n", event, goodTrack); | |
0ffa8579 | 204 | return kTRUE; |
99e5fe42 | 205 | } |
206 | ||
207 | ||
48bb52d8 | 208 | Float_t AliJetKineReader::SmearMomentum(Int_t ind, Float_t p) |
209 | { | |
83a444b1 | 210 | // The relative momentum error, i.e. |
211 | // (delta p)/p = sqrt (a**2 + (b*p)**2) * 10**-2, | |
212 | // where typically a = 0.75 and b = 0.16 - 0.24 depending on multiplicity | |
213 | // (the lower value is for dn/d(eta) about 2000, and the higher one for 8000) | |
214 | // | |
215 | // If we include information from TRD, b will be a factor 2/3 smaller. | |
216 | // | |
217 | // ind = 1: high multiplicity | |
218 | // ind = 2: low multiplicity | |
219 | // ind = 3: high multiplicity + TRD | |
220 | // ind = 4: low multiplicity + TRD | |
221 | ||
222 | Float_t pSmeared; | |
223 | Float_t a = 0.75; | |
224 | Float_t b = 0.12; | |
225 | ||
226 | if (ind == 1) b = 0.12; | |
227 | if (ind == 2) b = 0.08; | |
228 | if (ind == 3) b = 0.12; | |
229 | if (ind == 4) b = 0.08; | |
230 | ||
231 | Float_t sigma = p * TMath::Sqrt(a * a + b * b * p * p)*0.01; | |
232 | pSmeared = p + gRandom->Gaus(0., sigma); | |
233 | return pSmeared; | |
48bb52d8 | 234 | } |
83a444b1 | 235 | |
236 | ||
48bb52d8 | 237 | Bool_t AliJetKineReader::Efficiency(Float_t p, Float_t /*eta*/, Float_t phi) |
238 | { | |
83a444b1 | 239 | // Fast simulation of geometrical acceptance and tracking efficiency |
240 | ||
241 | // Tracking | |
242 | Float_t eff = 0.99; | |
243 | if (p < 0.5) eff -= (0.5-p)*0.2/0.3; | |
244 | // Geometry | |
245 | if (p > 0.5) { | |
246 | phi *= 180. / TMath::Pi(); | |
247 | // Sector number 0 - 17 | |
248 | Int_t isec = Int_t(phi / 20.); | |
249 | // Sector centre | |
250 | Float_t phi0 = isec * 20. + 10.; | |
251 | Float_t phir = TMath::Abs(phi-phi0); | |
252 | // 2 deg of dead space | |
253 | if (phir > 9.) eff = 0.; | |
254 | } | |
255 | ||
256 | if (gRandom->Rndm() > eff) { | |
257 | return kFALSE; | |
258 | } else { | |
259 | return kTRUE; | |
260 | } | |
261 | ||
48bb52d8 | 262 | } |
263 | ||
0ffa8579 | 264 | Bool_t AliJetKineReader::GetGenJets(AliJet* genJets) |
265 | { | |
266 | // Get generated jets from mc header | |
267 | AliHeader* alih = GetAliHeader(); | |
268 | if (alih == 0) return kFALSE; | |
269 | AliGenEventHeader * genh = alih->GenEventHeader(); | |
270 | if (genh == 0) return kFALSE; | |
271 | Int_t nj =((AliGenPythiaEventHeader*)genh)->NTriggerJets(); | |
272 | Int_t* m = new Int_t[nj]; | |
273 | Int_t* k = new Int_t[nj]; | |
274 | for (Int_t i=0; i< nj; i++) { | |
275 | Float_t p[4]; | |
276 | ((AliGenPythiaEventHeader*)genh)->TriggerJet(i,p); | |
277 | genJets->AddJet(p[0],p[1],p[2],p[3]); | |
278 | m[i]=1; | |
279 | k[i]=i; | |
280 | } | |
281 | genJets->SetNinput(nj); | |
282 | genJets->SetMultiplicities(m); | |
283 | genJets->SetInJet(k); | |
284 | return kTRUE; | |
285 | } |