]>
Commit | Line | Data |
---|---|---|
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 | ||
16 | /* $Id$ */ | |
17 | ||
18 | //-------------------------------------------------- | |
19 | // Filling of CalTrkTrack objects in the MC reader task | |
20 | // | |
21 | // Author: magali.estienne@subatech.in2p3.fr | |
22 | // alexandre.shabetai@cern.ch | |
23 | //------------------------------------------------- | |
24 | ||
25 | // --- ROOT system --- | |
26 | #include <TDatabasePDG.h> | |
27 | #include <TRandom.h> | |
28 | #include <TPDGCode.h> | |
29 | ||
30 | // --- AliRoot system --- | |
31 | #include "AliJetFillCalTrkTrackMC.h" | |
32 | #include "AliAODMCParticle.h" | |
33 | #include "AliJetKineReaderHeader.h" | |
34 | #include "AliVEvent.h" | |
35 | #include "AliMCEvent.h" | |
36 | ||
37 | using std::cout; | |
38 | using std::endl; | |
39 | ClassImp(AliJetFillCalTrkTrackMC) | |
40 | ||
41 | ///////////////////////////////////////////////////////////////////////// | |
42 | ||
43 | AliJetFillCalTrkTrackMC::AliJetFillCalTrkTrackMC(): | |
44 | AliJetFillCalTrkEvent(), | |
45 | fHadCorr(0), | |
46 | fApplyMIPCorrection(kTRUE), | |
47 | fVEvt(0x0), | |
48 | fMCEvent(0x0) | |
49 | { | |
50 | // constructor | |
51 | } | |
52 | ||
53 | //----------------------------------------------------------------------- | |
54 | AliJetFillCalTrkTrackMC::AliJetFillCalTrkTrackMC(AliVEvent* evt): | |
55 | AliJetFillCalTrkEvent(), | |
56 | fHadCorr(0), | |
57 | fApplyMIPCorrection(kTRUE), | |
58 | fVEvt(evt), | |
59 | fMCEvent(0x0) | |
60 | { | |
61 | // constructor | |
62 | } | |
63 | ||
64 | //----------------------------------------------------------------------- | |
65 | AliJetFillCalTrkTrackMC::AliJetFillCalTrkTrackMC(const AliJetFillCalTrkTrackMC &det): | |
66 | AliJetFillCalTrkEvent(det), | |
67 | fHadCorr(det.fHadCorr), | |
68 | fApplyMIPCorrection(det.fApplyMIPCorrection), | |
69 | fVEvt(det.fVEvt), | |
70 | fMCEvent(det.fMCEvent) | |
71 | { | |
72 | // Copy constructor | |
73 | } | |
74 | ||
75 | //----------------------------------------------------------------------- | |
76 | AliJetFillCalTrkTrackMC& AliJetFillCalTrkTrackMC::operator=(const AliJetFillCalTrkTrackMC& other) | |
77 | { | |
78 | // Assignment | |
79 | if (this != &other) { | |
80 | fHadCorr = other.fHadCorr; | |
81 | fApplyMIPCorrection = other.fApplyMIPCorrection; | |
82 | fVEvt = other.fVEvt; | |
83 | fMCEvent = other.fMCEvent; | |
84 | } | |
85 | ||
86 | return (*this); | |
87 | ||
88 | } | |
89 | ||
90 | //----------------------------------------------------------------------- | |
91 | AliJetFillCalTrkTrackMC::~AliJetFillCalTrkTrackMC() | |
92 | { | |
93 | // destructor | |
94 | } | |
95 | ||
96 | //----------------------------------------------------------------------- | |
97 | void AliJetFillCalTrkTrackMC::Exec(Option_t const * /*option*/) | |
98 | { | |
99 | // Main method. | |
100 | // Fill AliJetFillCalTrkTrackMC the with the charged particle information | |
101 | ||
102 | fDebug = fReaderHeader->GetDebug(); | |
103 | fOpt = fReaderHeader->GetDetector(); | |
104 | TString datatype = fReaderHeader->GetDataType(); | |
105 | datatype.ToUpper(); | |
106 | ||
107 | Int_t type=0; | |
108 | ||
109 | if (datatype.Contains("AODMC2B")) { type=kTrackAODMCChargedAcceptance; } | |
110 | else if (datatype.Contains("AODMC2")) { type=kTrackAODMCCharged; } | |
111 | else if (datatype.Contains("AODMC")) { type=kTrackAODMCAll; } | |
112 | ||
113 | else if (!datatype.Contains("AOD") && datatype.Contains("MC2")) {type=kTrackKineCharged;} | |
114 | else if (!datatype.Contains("AOD") && datatype.Contains("MC")) {type=kTrackKineAll;} | |
115 | ||
116 | else { cout<< "Unknown Data type !" << endl; } | |
117 | ||
118 | // temporary storage of signal and pt cut flag | |
119 | Bool_t sflag = 0; | |
120 | Bool_t cflag = 0; | |
121 | ||
122 | // get cuts set by user | |
123 | Float_t ptMin = fReaderHeader->GetPtCut(); | |
124 | Float_t etaMin = fReaderHeader->GetFiducialEtaMin(); | |
125 | Float_t etaMax = fReaderHeader->GetFiducialEtaMax(); | |
126 | Float_t phiMin = fReaderHeader->GetFiducialPhiMin(); | |
127 | Float_t phiMax = fReaderHeader->GetFiducialPhiMax(); | |
128 | ||
129 | if(fDebug>2) Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type); | |
130 | ||
131 | Int_t goodTrackStd = 0; | |
132 | Int_t goodTrackNonStd = 0; | |
133 | ||
134 | if (type == kTrackKineAll || type == kTrackKineCharged){ | |
135 | ||
136 | FillKine(); | |
137 | } // Kine | |
138 | ||
139 | else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) { | |
140 | if(!fVEvt) return; | |
141 | TClonesArray *mcArray = dynamic_cast<TClonesArray*>(fVEvt->FindListObject(AliAODMCParticle::StdBranchName())); | |
142 | if(!mcArray) { | |
143 | Printf("%s:%d No MC particle branch found",(char*)__FILE__,__LINE__); | |
144 | return; | |
145 | } | |
146 | for(int it = 0;it < mcArray->GetEntriesFast();++it){ | |
147 | cflag=sflag=0; | |
148 | AliAODMCParticle *part = (AliAODMCParticle*)(mcArray->At(it)); | |
149 | if(!part->IsPhysicalPrimary())continue; | |
150 | Int_t pdg = TMath::Abs(part->GetPdgCode()); | |
151 | // exclude neutrinos anyway | |
152 | if((pdg == 12 || pdg == 14 || pdg == 16)) continue; | |
153 | cflag = ( part->Pt()>ptMin ) ? 1 : 0; | |
154 | sflag = ( TMath::Abs(part->GetLabel()) < 10000 ) ? 1 : 0; | |
155 | if(type == kTrackAODMCAll){ | |
156 | fCalTrkEvent->AddCalTrkTrack(part,cflag,sflag); | |
157 | goodTrackStd++; | |
158 | } | |
159 | else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){ | |
160 | if(part->Charge()==0)continue; | |
161 | if(kTrackAODMCCharged){ | |
162 | fCalTrkEvent->AddCalTrkTrack(part,cflag,sflag); | |
163 | goodTrackStd++; | |
164 | } | |
165 | else { | |
166 | if((part->Eta()>etaMax) || (part->Eta()<etaMin)) continue; | |
167 | Float_t phi = ( (part->Phi()) < 0) ? (part->Phi()) + 2. * TMath::Pi() : (part->Phi()); | |
168 | if((phi>phiMax) || (phi<phiMin)) continue; | |
169 | fCalTrkEvent->AddCalTrkTrack(part,cflag,sflag); | |
170 | goodTrackStd++; | |
171 | } | |
172 | } | |
173 | } | |
174 | if(fDebug>0) printf("Number of good tracks selected: %5d \n", goodTrackStd+goodTrackNonStd); | |
175 | } // AODMCparticle | |
176 | ||
177 | } | |
178 | ||
179 | //----------------------------------------------------------------------- | |
180 | Bool_t AliJetFillCalTrkTrackMC::FillKine() | |
181 | { | |
182 | // Fill event | |
183 | Int_t goodTrack = 0; | |
184 | ||
185 | // Get the stack | |
186 | if(!fMCEvent) {cout<<"could not get MCEvent!"<<endl; return kFALSE;} | |
187 | ||
188 | Int_t nt = fMCEvent->GetNumberOfTracks(); | |
189 | ||
190 | // Get cuts set by user and header | |
191 | Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut(); | |
192 | Float_t etaMin = fReaderHeader->GetFiducialEtaMin(); | |
193 | Float_t etaMax = fReaderHeader->GetFiducialEtaMax(); | |
194 | ||
195 | TLorentzVector p4; | |
196 | // Loop over particles | |
197 | for (Int_t it = 0; it < nt; it++) { | |
198 | if(!fMCEvent->IsPhysicalPrimary(it)) continue; | |
199 | AliVParticle* part = fMCEvent->GetTrack(it); | |
200 | Int_t pdg = TMath::Abs(part->PdgCode()); | |
201 | Float_t pt = part->Pt(); | |
202 | ||
203 | if( (pdg == 12 || pdg == 14 || pdg == 16)) continue; | |
204 | ||
205 | Float_t p = part->P(); | |
206 | Float_t p0 = p; | |
207 | Float_t eta = part->Eta(); | |
208 | Float_t phi = part->Phi(); | |
209 | Float_t charge = part->Charge(); | |
210 | ||
211 | if (((AliJetKineReaderHeader*)fReaderHeader)->ChargedOnly()) { | |
212 | if (charge == 0) continue; | |
213 | } // End charged only | |
214 | ||
215 | // Fast simulation of EMCAL if requested | |
216 | if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) { | |
217 | // Charged particles only | |
218 | if (charge != 0){ | |
219 | // Simulate efficiency | |
220 | if (!Efficiency(p0, eta, phi)) continue; | |
221 | // Simulate resolution | |
222 | p = SmearMomentum(4, p0); | |
223 | } // end "if" charged particles | |
224 | // Neutral particles (exclude K0L, n, nbar) | |
225 | if (pdg == kNeutron || pdg == kK0Long) continue; | |
226 | } // End fast EMCAL | |
227 | ||
228 | // Fast simulation of TPC if requested | |
229 | if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) { | |
230 | // Charged particles only | |
231 | if (charge == 0) continue; | |
232 | // Simulate efficiency | |
233 | if (!Efficiency(p0, eta, phi)) continue; | |
234 | // Simulate resolution | |
235 | p = SmearMomentum(4, p0); | |
236 | } // End fast TPC | |
237 | ||
238 | // Fill momentum array | |
239 | Float_t r = p/p0; | |
240 | Float_t px = r * part->Px(); | |
241 | Float_t py = r * part->Py(); | |
242 | Float_t pz = r * part->Pz(); | |
243 | Float_t m = part->M(); | |
244 | ||
245 | Float_t e = TMath::Sqrt(px * px + py * py + pz * pz + m * m); | |
246 | p4 = TLorentzVector(px, py, pz, e); | |
247 | if ((p4.Eta()>etaMax) || (p4.Eta()<etaMin)) continue; | |
248 | ||
249 | // Flag used to store the r factor | |
250 | Float_t ptReso = r; | |
251 | Bool_t cflag = kFALSE, sflag = kTRUE; | |
252 | if (pt > ptMin) cflag = kTRUE; // track surviving pt cut | |
253 | ||
254 | fCalTrkEvent->AddCalTrkTrackKine(part,cflag,sflag,ptReso); | |
255 | ||
256 | goodTrack++; | |
257 | ||
258 | } // track loop | |
259 | ||
260 | if(fDebug>0) printf(" Number of good tracks %d \n", goodTrack); | |
261 | ||
262 | return kTRUE; | |
263 | ||
264 | } | |
265 | ||
266 | //----------------------------------------------------------------------- | |
267 | Float_t AliJetFillCalTrkTrackMC::SmearMomentum(Int_t ind, Float_t p) | |
268 | { | |
269 | // The relative momentum error, i.e. | |
270 | // (delta p)/p = sqrt (a**2 + (b*p)**2) * 10**-2, | |
271 | // where typically a = 0.75 and b = 0.16 - 0.24 depending on multiplicity | |
272 | // (the lower value is for dn/d(eta) about 2000, and the higher one for 8000) | |
273 | // | |
274 | // If we include information from TRD, b will be a factor 2/3 smaller. | |
275 | // | |
276 | // ind = 1: high multiplicity | |
277 | // ind = 2: low multiplicity | |
278 | // ind = 3: high multiplicity + TRD | |
279 | // ind = 4: low multiplicity + TRD | |
280 | ||
281 | Float_t pSmeared; | |
282 | Float_t a = 0.75; | |
283 | Float_t b = 0.12; | |
284 | ||
285 | if (ind == 1) b = 0.12; | |
286 | if (ind == 2) b = 0.08; | |
287 | if (ind == 3) b = 0.12; | |
288 | if (ind == 4) b = 0.08; | |
289 | ||
290 | Float_t sigma = p * TMath::Sqrt(a * a + b * b * p * p)*0.01; | |
291 | pSmeared = p + gRandom->Gaus(0., sigma); | |
292 | return pSmeared; | |
293 | ||
294 | } | |
295 | ||
296 | //----------------------------------------------------------------------- | |
297 | Bool_t AliJetFillCalTrkTrackMC::Efficiency(Float_t p, Float_t /*eta*/, Float_t phi) | |
298 | { | |
299 | // Fast simulation of geometrical acceptance and tracking efficiency | |
300 | ||
301 | // Tracking | |
302 | Float_t eff = 0.99; | |
303 | if (p < 0.5) eff -= (0.5-p)*0.2/0.3; | |
304 | // Geometry | |
305 | if (p > 0.5) { | |
306 | phi *= 180. / TMath::Pi(); | |
307 | // Sector number 0 - 17 | |
308 | Int_t isec = Int_t(phi / 20.); | |
309 | // Sector centre | |
310 | Float_t phi0 = isec * 20. + 10.; | |
311 | Float_t phir = TMath::Abs(phi-phi0); | |
312 | // 2 deg of dead space | |
313 | if (phir > 9.) eff = 0.; | |
314 | } | |
315 | ||
316 | if (gRandom->Rndm() > eff) { | |
317 | return kFALSE; | |
318 | } else { | |
319 | return kTRUE; | |
320 | } | |
321 | ||
322 | } | |
323 | ||
324 |