]>
Commit | Line | Data |
---|---|---|
762e8424 | 1 | // $Id$ |
2 | // | |
3 | // Jet modelling task. | |
4 | // | |
297edd60 | 5 | // Author: S.Aiola, C.Loizides |
762e8424 | 6 | |
ffe32451 | 7 | #include "AliJetModelBaseTask.h" |
8 | ||
762e8424 | 9 | #include <TClonesArray.h> |
787a3c4f | 10 | #include <TH1I.h> |
762e8424 | 11 | #include <TLorentzVector.h> |
12 | #include <TRandom3.h> | |
b16bb001 | 13 | #include <TList.h> |
43032ce2 | 14 | #include <TF2.h> |
762e8424 | 15 | |
b16bb001 | 16 | #include "AliVEvent.h" |
e44e8726 | 17 | #include "AliAODCaloCluster.h" |
b16bb001 | 18 | #include "AliESDCaloCluster.h" |
19 | #include "AliVCluster.h" | |
762e8424 | 20 | #include "AliEMCALDigit.h" |
ffe32451 | 21 | #include "AliEMCALRecPoint.h" |
b16bb001 | 22 | #include "AliESDCaloCells.h" |
23 | #include "AliAODCaloCells.h" | |
787a3c4f | 24 | #include "AliAODMCParticle.h" |
b16bb001 | 25 | #include "AliVCaloCells.h" |
ffe32451 | 26 | #include "AliPicoTrack.h" |
b16bb001 | 27 | #include "AliEMCALGeometry.h" |
28 | #include "AliLog.h" | |
5be3857d | 29 | #include "AliNamedArrayI.h" |
762e8424 | 30 | |
31 | ClassImp(AliJetModelBaseTask) | |
32 | ||
33 | //________________________________________________________________________ | |
34 | AliJetModelBaseTask::AliJetModelBaseTask() : | |
35 | AliAnalysisTaskSE("AliJetModelBaseTask"), | |
36 | fGeomName(), | |
37 | fTracksName(), | |
38 | fOutTracksName(), | |
39 | fCaloName(), | |
40 | fOutCaloName(), | |
b16bb001 | 41 | fCellsName(), |
42 | fOutCellsName(), | |
787a3c4f | 43 | fMCParticlesName(), |
44 | fOutMCParticlesName(), | |
4358e58a | 45 | fIsMC(kFALSE), |
762e8424 | 46 | fSuffix(), |
47 | fEtaMin(-1), | |
48 | fEtaMax(1), | |
49 | fPhiMin(0), | |
50 | fPhiMax(TMath::Pi() * 2), | |
51 | fPtMin(0), | |
52 | fPtMax(0), | |
53 | fCopyArray(kTRUE), | |
54 | fNClusters(0), | |
b16bb001 | 55 | fNCells(0), |
762e8424 | 56 | fNTracks(0), |
7f76e479 | 57 | fMarkMC(99999), |
e44e8726 | 58 | fPtSpectrum(0), |
43032ce2 | 59 | fPtPhiEvPlDistribution(0), |
27f77648 | 60 | fDensitySpectrum(0), |
50da5cf2 | 61 | fDifferentialV2(0), |
62 | fAddV2(kFALSE), | |
63 | fFlowFluctuations(kFALSE), | |
b16bb001 | 64 | fQAhistos(kFALSE), |
43032ce2 | 65 | fPsi(0), |
ffe32451 | 66 | fIsInit(0), |
762e8424 | 67 | fGeom(0), |
68 | fClusters(0), | |
69 | fOutClusters(0), | |
70 | fTracks(0), | |
b16bb001 | 71 | fOutTracks(0), |
72 | fCaloCells(0), | |
73 | fOutCaloCells(0), | |
74 | fAddedCells(0), | |
787a3c4f | 75 | fMCParticles(0), |
76 | fMCParticlesMap(0), | |
77 | fOutMCParticles(0), | |
78 | fOutMCParticlesMap(0), | |
79 | fMCLabelShift(0), | |
b16bb001 | 80 | fEsdMode(kFALSE), |
81 | fOutput(0) | |
762e8424 | 82 | { |
83 | // Default constructor. | |
8e48f7f5 | 84 | |
85 | fVertex[0] = 0; | |
86 | fVertex[1] = 0; | |
87 | fVertex[2] = 0; | |
762e8424 | 88 | } |
89 | ||
90 | //________________________________________________________________________ | |
b16bb001 | 91 | AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) : |
762e8424 | 92 | AliAnalysisTaskSE(name), |
93 | fGeomName(""), | |
94 | fTracksName("PicoTracks"), | |
95 | fOutTracksName("PicoTracksEmbedded"), | |
96 | fCaloName("CaloClustersCorr"), | |
97 | fOutCaloName("CaloClustersCorrEmbedded"), | |
b16bb001 | 98 | fCellsName(""), |
99 | fOutCellsName(""), | |
787a3c4f | 100 | fMCParticlesName(""), |
101 | fOutMCParticlesName(""), | |
4358e58a | 102 | fIsMC(kFALSE), |
762e8424 | 103 | fSuffix("Processed"), |
104 | fEtaMin(-1), | |
105 | fEtaMax(1), | |
106 | fPhiMin(0), | |
107 | fPhiMax(TMath::Pi() * 2), | |
108 | fPtMin(50), | |
109 | fPtMax(60), | |
110 | fCopyArray(kTRUE), | |
111 | fNClusters(0), | |
b16bb001 | 112 | fNCells(0), |
762e8424 | 113 | fNTracks(1), |
7f76e479 | 114 | fMarkMC(99999), |
e44e8726 | 115 | fPtSpectrum(0), |
43032ce2 | 116 | fPtPhiEvPlDistribution(0), |
27f77648 | 117 | fDensitySpectrum(0), |
50da5cf2 | 118 | fDifferentialV2(0), |
119 | fAddV2(kFALSE), | |
120 | fFlowFluctuations(kFALSE), | |
b16bb001 | 121 | fQAhistos(drawqa), |
43032ce2 | 122 | fPsi(0), |
ffe32451 | 123 | fIsInit(0), |
762e8424 | 124 | fGeom(0), |
125 | fClusters(0), | |
126 | fOutClusters(0), | |
127 | fTracks(0), | |
b16bb001 | 128 | fOutTracks(0), |
129 | fCaloCells(0), | |
130 | fOutCaloCells(0), | |
131 | fAddedCells(0), | |
787a3c4f | 132 | fMCParticles(0), |
133 | fMCParticlesMap(0), | |
134 | fOutMCParticles(0), | |
135 | fOutMCParticlesMap(0), | |
136 | fMCLabelShift(0), | |
b16bb001 | 137 | fEsdMode(kFALSE), |
138 | fOutput(0) | |
762e8424 | 139 | { |
140 | // Standard constructor. | |
b16bb001 | 141 | |
142 | if (fQAhistos) { | |
143 | DefineOutput(1, TList::Class()); | |
144 | } | |
8e48f7f5 | 145 | |
146 | fVertex[0] = 0; | |
147 | fVertex[1] = 0; | |
148 | fVertex[2] = 0; | |
762e8424 | 149 | } |
150 | ||
151 | //________________________________________________________________________ | |
152 | AliJetModelBaseTask::~AliJetModelBaseTask() | |
153 | { | |
154 | // Destructor | |
155 | } | |
156 | ||
b16bb001 | 157 | //________________________________________________________________________ |
158 | void AliJetModelBaseTask::UserCreateOutputObjects() | |
159 | { | |
160 | // Create user output. | |
161 | if (!fQAhistos) | |
162 | return; | |
163 | ||
164 | OpenFile(1); | |
165 | fOutput = new TList(); | |
166 | fOutput->SetOwner(); | |
167 | ||
168 | PostData(1, fOutput); | |
169 | } | |
170 | ||
762e8424 | 171 | //________________________________________________________________________ |
ffe32451 | 172 | void AliJetModelBaseTask::UserExec(Option_t *) |
173 | { | |
174 | // Execute per event. | |
175 | ||
d63c8c07 | 176 | if (!fIsInit) |
177 | fIsInit = ExecOnce(); | |
178 | ||
179 | if (!fIsInit) | |
180 | return; | |
181 | ||
6fd2f1a9 | 182 | fVertex[0] = 0; |
183 | fVertex[1] = 0; | |
184 | fVertex[2] = 0; | |
185 | ||
186 | const AliVVertex *vert = InputEvent()->GetPrimaryVertex(); | |
187 | if (vert) | |
188 | vert->GetXYZ(fVertex); | |
189 | ||
507f6e75 | 190 | if (fCopyArray) { |
191 | if (fOutTracks) | |
192 | fOutTracks->Delete(); | |
193 | if (fOutClusters) | |
194 | fOutClusters->Delete(); | |
195 | if (fOutMCParticles) | |
196 | fOutMCParticles->Delete(); | |
197 | } | |
27f77648 | 198 | |
199 | if (fDensitySpectrum) { | |
43032ce2 | 200 | fNTracks = TMath::Nint(fDensitySpectrum->GetRandom()); |
201 | fNCells = TMath::Nint(fDensitySpectrum->GetRandom()); | |
202 | fNClusters = TMath::Nint(fDensitySpectrum->GetRandom()); | |
27f77648 | 203 | } |
63fac07f | 204 | |
5be3857d | 205 | // Clear map |
206 | if (fOutMCParticlesMap) | |
4358e58a | 207 | fOutMCParticlesMap->Clear(); |
d63c8c07 | 208 | |
cd6431de | 209 | AliVCaloCells *tempCaloCells = 0; |
210 | ||
b16bb001 | 211 | if (fCaloCells) { |
212 | fAddedCells = 0; | |
cd6431de | 213 | if (!fCopyArray) { |
214 | tempCaloCells = fCaloCells; | |
215 | fCaloCells = static_cast<AliVCaloCells*>(tempCaloCells->Clone(Form("%s_old",fCaloCells->GetName()))); | |
216 | } | |
b16bb001 | 217 | } |
218 | ||
50da5cf2 | 219 | if (fPtPhiEvPlDistribution || fAddV2) |
43032ce2 | 220 | fPsi = gRandom->Rndm() * TMath::Pi(); |
221 | ||
ffe32451 | 222 | Run(); |
cd6431de | 223 | |
fde82e42 | 224 | if (fCaloCells && !fCopyArray) { |
cd6431de | 225 | delete fCaloCells; |
226 | fCaloCells = tempCaloCells; | |
227 | } | |
ffe32451 | 228 | } |
229 | ||
b16bb001 | 230 | //________________________________________________________________________ |
231 | Bool_t AliJetModelBaseTask::ExecOnce() | |
232 | { | |
233 | // Init task. | |
234 | ||
235 | delete gRandom; | |
236 | gRandom = new TRandom3(0); | |
237 | ||
238 | fEsdMode = InputEvent()->InheritsFrom("AliESDEvent"); | |
239 | ||
240 | if (fPtMax < fPtMin) { | |
241 | AliWarning (Form("PtMax (%f) < PtMin (%f), setting PtMax = PtMin = %f", fPtMax, fPtMin, fPtMin)); | |
242 | fPtMax = fPtMin; | |
243 | } | |
244 | ||
245 | if (fEtaMax < fEtaMin) { | |
246 | AliWarning (Form("EtaMax (%f) < EtaMin (%f), setting EtaMax = EtaMin = %f", fEtaMax, fEtaMin, fEtaMin)); | |
247 | fEtaMax = fEtaMin; | |
248 | } | |
249 | ||
250 | if (fPhiMax < fPhiMin) { | |
251 | AliWarning (Form("PhiMax (%f) < PhiMin (%f), setting PhiMax = PhiMin = %f", fPhiMax, fPhiMin, fPhiMin)); | |
252 | fPhiMax = fPhiMin; | |
253 | } | |
254 | ||
787a3c4f | 255 | if (!fCellsName.IsNull()) { |
b16bb001 | 256 | fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCellsName)); |
257 | if (!fCaloCells) { | |
787a3c4f | 258 | AliWarning(Form("%s: Couldn't retrieve calo cells with name %s!", GetName(), fCellsName.Data())); |
b16bb001 | 259 | } |
787a3c4f | 260 | else if (!fCaloCells->InheritsFrom("AliVCaloCells")) { |
b16bb001 | 261 | AliError(Form("%s: Collection %s does not contain a AliVCaloCells object!", GetName(), fCellsName.Data())); |
262 | fCaloCells = 0; | |
263 | return kFALSE; | |
264 | } | |
265 | ||
266 | if (!fOutCaloCells) { | |
267 | fOutCellsName = fCellsName; | |
787a3c4f | 268 | if (fCopyArray) |
b16bb001 | 269 | fOutCellsName += fSuffix; |
787a3c4f | 270 | if (fCopyArray || !fCaloCells) { |
b16bb001 | 271 | if (fEsdMode) |
272 | fOutCaloCells = new AliESDCaloCells(fOutCellsName,fOutCellsName); | |
273 | else | |
274 | fOutCaloCells = new AliAODCaloCells(fOutCellsName,fOutCellsName); | |
275 | ||
276 | if (InputEvent()->FindListObject(fOutCellsName)) { | |
277 | AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCellsName.Data())); | |
278 | return kFALSE; | |
279 | } | |
280 | else { | |
281 | InputEvent()->AddObject(fOutCaloCells); | |
282 | } | |
283 | } | |
284 | else { | |
285 | fOutCaloCells = fCaloCells; | |
286 | } | |
287 | } | |
288 | } | |
289 | ||
27f77648 | 290 | if (!fTracksName.IsNull()) { |
b16bb001 | 291 | fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName)); |
292 | if (!fTracks) { | |
787a3c4f | 293 | AliWarning(Form("%s: Couldn't retrieve tracks with name %s!", GetName(), fTracksName.Data())); |
b16bb001 | 294 | } |
787a3c4f | 295 | else if (!fTracks->GetClass()->GetBaseClass("AliPicoTrack")) { |
b16bb001 | 296 | AliError(Form("%s: Collection %s does not contain AliPicoTrack objects!", GetName(), fTracksName.Data())); |
297 | fTracks = 0; | |
298 | return kFALSE; | |
299 | } | |
300 | ||
301 | if (!fOutTracks) { | |
302 | fOutTracksName = fTracksName; | |
787a3c4f | 303 | if (fCopyArray) |
b16bb001 | 304 | fOutTracksName += fSuffix; |
787a3c4f | 305 | if (fCopyArray || !fTracks) { |
306 | fOutTracks = new TClonesArray("AliPicoTrack"); | |
b16bb001 | 307 | fOutTracks->SetName(fOutTracksName); |
308 | if (InputEvent()->FindListObject(fOutTracksName)) { | |
309 | AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutTracksName.Data())); | |
310 | return kFALSE; | |
311 | } | |
312 | else { | |
313 | InputEvent()->AddObject(fOutTracks); | |
314 | } | |
315 | } | |
316 | else { | |
317 | fOutTracks = fTracks; | |
318 | } | |
319 | } | |
320 | } | |
321 | ||
50da5cf2 | 322 | if(fAddV2 && (!fDifferentialV2)) { |
323 | AliWarning(Form("%s: Cannot add v2 without diffential v2!", GetName())); | |
324 | } | |
325 | ||
27f77648 | 326 | if (!fCaloName.IsNull()) { |
b16bb001 | 327 | fClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName)); |
328 | ||
329 | if (!fClusters) { | |
787a3c4f | 330 | AliWarning(Form("%s: Couldn't retrieve clusters with name %s!", GetName(), fCaloName.Data())); |
b16bb001 | 331 | } |
787a3c4f | 332 | else if (!fClusters->GetClass()->GetBaseClass("AliVCluster")) { |
b16bb001 | 333 | AliError(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloName.Data())); |
334 | fClusters = 0; | |
335 | return kFALSE; | |
336 | } | |
337 | ||
338 | if (!fOutClusters) { | |
339 | fOutCaloName = fCaloName; | |
787a3c4f | 340 | if (fCopyArray) |
b16bb001 | 341 | fOutCaloName += fSuffix; |
507f74bc | 342 | TString className; |
343 | if (fClusters) | |
344 | className = fClusters->GetClass()->GetName(); | |
345 | else if (fEsdMode) | |
346 | className = "AliESDCaloCluster"; | |
347 | else | |
348 | className = "AliAODCaloCluster"; | |
349 | ||
787a3c4f | 350 | if (fCopyArray || !fClusters) { |
507f74bc | 351 | fOutClusters = new TClonesArray(className.Data()); |
b16bb001 | 352 | fOutClusters->SetName(fOutCaloName); |
353 | if (InputEvent()->FindListObject(fOutCaloName)) { | |
354 | AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCaloName.Data())); | |
355 | return kFALSE; | |
356 | } | |
357 | else { | |
358 | InputEvent()->AddObject(fOutClusters); | |
359 | } | |
360 | } | |
361 | else { | |
362 | fOutClusters = fClusters; | |
363 | } | |
364 | } | |
365 | } | |
366 | ||
787a3c4f | 367 | if (!fMCParticlesName.IsNull()) { |
368 | fMCParticles = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCParticlesName)); | |
369 | if (!fMCParticles) { | |
370 | AliWarning(Form("%s: Couldn't retrieve MC particles with name %s!", GetName(), fMCParticlesName.Data())); | |
371 | } | |
372 | else { | |
373 | if (!fMCParticles->GetClass()->GetBaseClass("AliAODMCParticle")) { | |
374 | AliError(Form("%s: Collection %s does not contain AliAODMCParticle objects!", GetName(), fMCParticlesName.Data())); | |
375 | fMCParticles = 0; | |
376 | return kFALSE; | |
377 | } | |
378 | ||
5be3857d | 379 | fMCParticlesMap = dynamic_cast<AliNamedArrayI*>(InputEvent()->FindListObject(fMCParticlesName + "_Map")); |
787a3c4f | 380 | |
381 | if (!fMCParticlesMap) { | |
382 | AliWarning(Form("%s: Could not retrieve map for MC particles %s! Will assume MC labels consistent with indexes...", GetName(), fMCParticlesName.Data())); | |
a7477843 | 383 | fMCParticlesMap = new AliNamedArrayI(fMCParticlesName + "_Map", 99999); |
384 | for (Int_t i = 0; i < 99999; i++) { | |
5be3857d | 385 | fMCParticlesMap->AddAt(i,i); |
787a3c4f | 386 | } |
387 | } | |
388 | } | |
389 | ||
390 | if (!fOutMCParticles) { | |
391 | fOutMCParticlesName = fMCParticlesName; | |
392 | if (fCopyArray) | |
393 | fOutMCParticlesName += fSuffix; | |
394 | if (fCopyArray || !fMCParticles) { | |
395 | ||
396 | fOutMCParticles = new TClonesArray("AliAODMCParticle"); | |
397 | fOutMCParticles->SetName(fOutMCParticlesName); | |
398 | if (InputEvent()->FindListObject(fOutMCParticlesName)) { | |
399 | AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutMCParticlesName.Data())); | |
400 | return kFALSE; | |
401 | } | |
402 | else { | |
403 | InputEvent()->AddObject(fOutMCParticles); | |
404 | } | |
405 | ||
a7477843 | 406 | fOutMCParticlesMap = new AliNamedArrayI(fOutMCParticlesName + "_Map",99999); |
787a3c4f | 407 | if (InputEvent()->FindListObject(fOutMCParticlesName + "_Map")) { |
408 | AliFatal(Form("%s: Map %s_Map is already present in the event!", GetName(), fOutMCParticlesName.Data())); | |
409 | return kFALSE; | |
410 | } | |
411 | else { | |
412 | InputEvent()->AddObject(fOutMCParticlesMap); | |
413 | } | |
414 | } | |
415 | else { | |
416 | fOutMCParticles = fMCParticles; | |
417 | fOutMCParticlesMap = fMCParticlesMap; | |
418 | } | |
419 | } | |
420 | } | |
421 | ||
2103dc6a | 422 | if (!fGeom) { |
423 | if (fGeomName.Length() > 0) { | |
424 | fGeom = AliEMCALGeometry::GetInstance(fGeomName); | |
425 | if (!fGeom) { | |
426 | AliFatal(Form("Could not get geometry with name %s!", fGeomName.Data())); | |
427 | return kFALSE; | |
428 | } | |
429 | } else { | |
4358e58a | 430 | fGeom = AliEMCALGeometry::GetInstance(); |
431 | if (!fGeom) { | |
432 | AliFatal("Could not get default geometry!"); | |
433 | return kFALSE; | |
434 | } | |
b16bb001 | 435 | } |
cfc2ac24 | 436 | } |
b16bb001 | 437 | |
438 | return kTRUE; | |
439 | } | |
440 | ||
441 | //________________________________________________________________________ | |
442 | Int_t AliJetModelBaseTask::SetNumberOfOutCells(Int_t n) | |
443 | { | |
444 | if (fOutCaloCells->GetNumberOfCells() < n) { | |
445 | fOutCaloCells->DeleteContainer(); | |
446 | fOutCaloCells->CreateContainer(n); | |
447 | } | |
448 | else { | |
449 | fOutCaloCells->SetNumberOfCells(n); | |
450 | } | |
451 | ||
452 | fAddedCells = 0; | |
453 | ||
454 | return n; | |
455 | } | |
456 | ||
b16bb001 | 457 | //________________________________________________________________________ |
458 | Int_t AliJetModelBaseTask::AddCell(Double_t e, Double_t eta, Double_t phi) | |
459 | { | |
460 | // Add a cell to the event. | |
461 | ||
462 | Int_t absId = 0; | |
463 | if (eta < -100 || phi < 0) { | |
464 | GetRandomCell(eta, phi, absId); | |
465 | } | |
466 | else { | |
467 | fGeom->EtaPhiFromIndex(absId, eta, phi); | |
468 | } | |
469 | ||
470 | if (absId == -1) { | |
471 | AliWarning(Form("Unable to embed cell in eta = %f, phi = %f!" | |
472 | " Maybe the eta-phi range is not inside the EMCal acceptance (eta = [%f, %f], phi = [%f, %f])", | |
473 | eta, phi, fEtaMin, fEtaMax, fPhiMin, fPhiMax)); | |
474 | return 0; | |
475 | } | |
476 | ||
477 | if (e < 0) { | |
478 | Double_t pt = GetRandomPt(); | |
479 | TLorentzVector nPart; | |
480 | nPart.SetPtEtaPhiM(pt, eta, phi, 0); | |
481 | e = nPart.E(); | |
482 | } | |
483 | ||
484 | return AddCell(e, absId); | |
485 | } | |
486 | ||
487 | //________________________________________________________________________ | |
787a3c4f | 488 | Int_t AliJetModelBaseTask::AddCell(Double_t e, Int_t absId, Double_t time, Int_t label) |
b16bb001 | 489 | { |
490 | // Add a cell to the event. | |
491 | ||
a7477843 | 492 | if (label < 0) |
493 | label = 0; | |
494 | ||
495 | label += fMarkMC + fMCLabelShift; | |
787a3c4f | 496 | |
fde82e42 | 497 | Short_t pos = -1; |
498 | if (fCaloCells) | |
499 | pos = fCaloCells->GetCellPosition(absId); | |
500 | ||
4358e58a | 501 | Double_t efrac = 1; |
502 | Bool_t increaseOnSuccess = kFALSE; | |
503 | ||
504 | if (pos < 0) { | |
505 | increaseOnSuccess = kTRUE; | |
506 | pos = fAddedCells; | |
507 | } | |
508 | else { | |
509 | Short_t cellNumber = -1; | |
510 | Double_t old_e = 0; | |
511 | Double_t old_time = 0; | |
512 | Int_t old_label = 0; | |
513 | Double_t old_efrac = 0; | |
514 | fOutCaloCells->GetCell(pos, cellNumber, old_e, old_time, old_label, old_efrac); | |
515 | ||
fde82e42 | 516 | efrac = e / (old_e + e); |
517 | ||
518 | if (old_label > 0 && e < old_efrac * old_e) { | |
4358e58a | 519 | label = old_label; |
fde82e42 | 520 | efrac = old_efrac; |
4358e58a | 521 | time = old_time; |
522 | } | |
523 | ||
4358e58a | 524 | e += old_e; |
525 | } | |
526 | ||
527 | Bool_t r = fOutCaloCells->SetCell(pos, absId, e, time, label, efrac); | |
b16bb001 | 528 | |
529 | if (r) { | |
4358e58a | 530 | if (increaseOnSuccess) |
531 | fAddedCells++; | |
b16bb001 | 532 | return fAddedCells; |
533 | } | |
4358e58a | 534 | else |
b16bb001 | 535 | return 0; |
b16bb001 | 536 | } |
537 | ||
507f74bc | 538 | //________________________________________________________________________ |
539 | AliVCluster* AliJetModelBaseTask::AddCluster(AliVCluster *oc) | |
540 | { | |
541 | // Add a cluster to the event. | |
542 | ||
543 | const Int_t nClusters = fOutClusters->GetEntriesFast(); | |
544 | AliVCluster *dc = static_cast<AliVCluster*>(fOutClusters->New(nClusters)); | |
545 | dc->SetType(AliVCluster::kEMCALClusterv1); | |
546 | dc->SetE(oc->E()); | |
547 | Float_t pos[3] = {0}; | |
548 | oc->GetPosition(pos); | |
549 | dc->SetPosition(pos); | |
550 | dc->SetNCells(oc->GetNCells()); | |
551 | dc->SetCellsAbsId(oc->GetCellsAbsId()); | |
552 | dc->SetCellsAmplitudeFraction(oc->GetCellsAmplitudeFraction()); | |
553 | dc->SetID(oc->GetID()); | |
554 | dc->SetDispersion(oc->GetDispersion()); | |
555 | dc->SetEmcCpvDistance(-1); | |
556 | dc->SetChi2(-1); | |
557 | dc->SetTOF(oc->GetTOF()); //time-of-flight | |
558 | dc->SetNExMax(oc->GetNExMax()); //number of local maxima | |
559 | dc->SetM02(oc->GetM02()); | |
560 | dc->SetM20(oc->GetM20()); | |
561 | dc->SetDistanceToBadChannel(oc->GetDistanceToBadChannel()); | |
562 | ||
563 | //MC | |
564 | UInt_t nlabels = oc->GetNLabels(); | |
565 | Int_t *labels = oc->GetLabels(); | |
a7477843 | 566 | TArrayI parents; |
567 | if (nlabels != 0 && labels) | |
568 | parents.Set(nlabels, labels); | |
569 | else { | |
570 | nlabels = 1; | |
571 | parents.Set(1); | |
572 | parents[0] = 0; | |
573 | } | |
507f74bc | 574 | |
a7477843 | 575 | if (fMarkMC+fMCLabelShift != 0) { |
576 | for (UInt_t i = 0; i < nlabels; i++) { | |
577 | parents[i] += fMarkMC+fMCLabelShift; | |
507f74bc | 578 | } |
579 | } | |
a7477843 | 580 | |
581 | AliESDCaloCluster *esdClus = dynamic_cast<AliESDCaloCluster*>(dc); | |
582 | if (esdClus) { | |
583 | esdClus->AddLabels(parents); | |
584 | } | |
585 | else { | |
586 | AliAODCaloCluster *aodClus = dynamic_cast<AliAODCaloCluster*>(dc); | |
587 | if (aodClus) | |
588 | aodClus->SetLabel(parents.GetArray(), nlabels); | |
fde82e42 | 589 | } |
507f74bc | 590 | |
591 | return dc; | |
592 | } | |
b16bb001 | 593 | |
ffe32451 | 594 | //________________________________________________________________________ |
787a3c4f | 595 | AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Double_t eta, Double_t phi, Int_t label) |
ffe32451 | 596 | { |
597 | // Add a cluster to the event. | |
598 | ||
599 | Int_t absId = 0; | |
600 | if (eta < -100 || phi < 0) { | |
601 | GetRandomCell(eta, phi, absId); | |
602 | } | |
603 | else { | |
604 | fGeom->EtaPhiFromIndex(absId, eta, phi); | |
605 | } | |
606 | ||
607 | if (absId == -1) { | |
608 | AliWarning(Form("Unable to embed cluster in eta = %f, phi = %f!" | |
609 | " Maybe the eta-phi range is not inside the EMCal acceptance (eta = [%f, %f], phi = [%f, %f])", | |
610 | eta, phi, fEtaMin, fEtaMax, fPhiMin, fPhiMax)); | |
611 | return 0; | |
612 | } | |
613 | ||
614 | if (e < 0) { | |
615 | Double_t pt = GetRandomPt(); | |
616 | TLorentzVector nPart; | |
617 | nPart.SetPtEtaPhiM(pt, eta, phi, 0); | |
618 | e = nPart.E(); | |
619 | } | |
620 | ||
787a3c4f | 621 | return AddCluster(e, absId, label); |
ffe32451 | 622 | } |
623 | ||
624 | //________________________________________________________________________ | |
787a3c4f | 625 | AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Int_t absId, Int_t label) |
ffe32451 | 626 | { |
627 | // Add a cluster to the event. | |
628 | ||
629 | const Int_t nClusters = fOutClusters->GetEntriesFast(); | |
630 | ||
631 | TClonesArray digits("AliEMCALDigit", 1); | |
632 | ||
633 | AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(digits.New(0)); | |
634 | digit->SetId(absId); | |
635 | digit->SetIndexInList(0); | |
636 | digit->SetType(AliEMCALDigit::kHG); | |
637 | digit->SetAmplitude(e); | |
638 | ||
639 | AliEMCALRecPoint recPoint(""); | |
640 | recPoint.AddDigit(*digit, e, kFALSE); | |
641 | recPoint.EvalGlobalPosition(0, &digits); | |
642 | ||
643 | TVector3 gpos; | |
644 | recPoint.GetGlobalPosition(gpos); | |
645 | Float_t g[3]; | |
646 | gpos.GetXYZ(g); | |
647 | ||
648 | AliVCluster *cluster = static_cast<AliVCluster*>(fOutClusters->New(nClusters)); | |
649 | cluster->SetType(AliVCluster::kEMCALClusterv1); | |
650 | cluster->SetE(recPoint.GetEnergy()); | |
651 | cluster->SetPosition(g); | |
652 | cluster->SetNCells(1); | |
653 | UShort_t shortAbsId = absId; | |
654 | cluster->SetCellsAbsId(&shortAbsId); | |
655 | Double32_t fract = 1; | |
656 | cluster->SetCellsAmplitudeFraction(&fract); | |
657 | cluster->SetID(nClusters); | |
658 | cluster->SetEmcCpvDistance(-1); | |
ffe32451 | 659 | |
f660c2d6 | 660 | //MC label |
a7477843 | 661 | if (label < 0) |
662 | label = 0; | |
663 | ||
664 | label += fMarkMC+fMCLabelShift; | |
787a3c4f | 665 | |
f660c2d6 | 666 | if (fEsdMode) { |
667 | AliESDCaloCluster *esdClus = static_cast<AliESDCaloCluster*>(cluster); | |
787a3c4f | 668 | TArrayI parents(1, &label); |
f660c2d6 | 669 | esdClus->AddLabels(parents); |
670 | } | |
671 | else { | |
672 | AliAODCaloCluster *aodClus = static_cast<AliAODCaloCluster*>(cluster); | |
787a3c4f | 673 | aodClus->SetLabel(&label, 1); |
f660c2d6 | 674 | } |
675 | ||
ffe32451 | 676 | return cluster; |
677 | } | |
678 | ||
679 | //________________________________________________________________________ | |
396ffc82 | 680 | AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t phi, Byte_t type, Double_t etaemc, Double_t phiemc, Double_t ptemc, Bool_t ise, Int_t label, Short_t charge, Double_t mass) |
ffe32451 | 681 | { |
682 | // Add a track to the event. | |
ffe32451 | 683 | |
43032ce2 | 684 | if (pt < 0 && eta < -100 && phi < 0) { |
685 | GetRandomParticle(pt,eta,phi); | |
686 | } | |
687 | else { | |
688 | if (pt < 0) | |
689 | pt = GetRandomPt(); | |
690 | if (eta < -100) | |
691 | eta = GetRandomEta(); | |
692 | if (phi < 0) | |
693 | phi = GetRandomPhi(pt); | |
694 | } | |
ffe32451 | 695 | |
a7477843 | 696 | if (label >= 0) |
697 | label += fMarkMC+fMCLabelShift; | |
698 | else if (label < 0) | |
699 | label -= fMarkMC+fMCLabelShift; | |
43032ce2 | 700 | |
50da5cf2 | 701 | if(fAddV2) AddV2(phi, pt); |
702 | ||
43032ce2 | 703 | const Int_t nTracks = fOutTracks->GetEntriesFast(); |
787a3c4f | 704 | |
ffe32451 | 705 | AliPicoTrack *track = new ((*fOutTracks)[nTracks]) AliPicoTrack(pt, |
706 | eta, | |
707 | phi, | |
6fd2f1a9 | 708 | charge, |
787a3c4f | 709 | label, |
b16bb001 | 710 | type, |
711 | etaemc, | |
56bd3193 | 712 | phiemc, |
713 | ptemc, | |
396ffc82 MV |
714 | ise, |
715 | mass); | |
b16bb001 | 716 | |
ffe32451 | 717 | return track; |
718 | } | |
719 | ||
787a3c4f | 720 | //________________________________________________________________________ |
721 | AliAODMCParticle* AliJetModelBaseTask::AddMCParticle(AliAODMCParticle *part, Int_t origIndex) | |
722 | { | |
723 | const Int_t nPart = fOutMCParticles->GetEntriesFast(); | |
724 | ||
725 | AliAODMCParticle *aodpart = new ((*fOutMCParticles)[nPart]) AliAODMCParticle(*part); | |
726 | ||
a7477843 | 727 | if (origIndex + fMCLabelShift >= fOutMCParticlesMap->GetSize()) |
728 | fOutMCParticlesMap->Set((origIndex + fMCLabelShift)*2); | |
729 | ||
5be3857d | 730 | fOutMCParticlesMap->AddAt(nPart, origIndex + fMCLabelShift); |
5ce8ae64 | 731 | AliDebug(2, Form("Setting bin %d to %d (fMCLabelShift=%d, origIndex=%d)", |
732 | origIndex + fMCLabelShift, fOutMCParticlesMap->At(origIndex + fMCLabelShift), fMCLabelShift, origIndex)); | |
787a3c4f | 733 | |
734 | return aodpart; | |
735 | } | |
736 | ||
50da5cf2 | 737 | //_____________________________________________________________________________ |
738 | void AliJetModelBaseTask::AddV2(Double_t &phi, Double_t &pt) const | |
739 | { | |
740 | // similar to AliFlowTrackSimple::AddV2, except for the flow fluctuations | |
741 | Double_t phi0(phi), v2(0.), f(0.), fp(0.), phiprev(0.); | |
742 | if(fDifferentialV2) v2 = fDifferentialV2->Eval(pt); | |
743 | if(TMath::AreEqualAbs(v2, 0, 1e-5)) return; | |
744 | // introduce flow fluctuations (gaussian) | |
745 | if(fFlowFluctuations) v2 += TMath::Sqrt(2*(v2*.25)*(v2*.25))*TMath::ErfInverse(2*(gRandom->Uniform(0, fFlowFluctuations))-1); | |
746 | for (Int_t i(0); i < 100; i++) { | |
747 | phiprev=phi; //store last value for comparison | |
748 | f = phi-phi0+v2*TMath::Sin(2.*(phi-fPsi)); | |
749 | fp = 1.0+2.0*v2*TMath::Cos(2.*(phi-fPsi)); //first derivative | |
750 | phi -= f/fp; | |
751 | if (TMath::AreEqualAbs(phiprev, phi, 1e-10)) break; | |
752 | } | |
753 | } | |
754 | ||
755 | //_____________________________________________________________________________ | |
4358e58a | 756 | void AliJetModelBaseTask::CopyCells() |
757 | { | |
758 | if (!fCaloCells) | |
759 | return; | |
760 | ||
fde82e42 | 761 | fAddedCells = 0; |
762 | fCaloCells->Sort(); | |
4358e58a | 763 | for (Short_t i = 0; i < fCaloCells->GetNumberOfCells(); i++) { |
764 | Int_t mclabel = 0; | |
765 | Double_t efrac = 0.; | |
766 | Double_t time = -1; | |
767 | Short_t cellNum = -1; | |
768 | Double_t amp = -1; | |
769 | ||
770 | fCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac); | |
771 | ||
772 | if (!fIsMC) | |
773 | mclabel = 0; | |
774 | ||
775 | fOutCaloCells->SetCell(i, cellNum, amp, time, mclabel, efrac); | |
fde82e42 | 776 | fAddedCells++; |
4358e58a | 777 | } |
778 | ||
4358e58a | 779 | AliDebug(2, Form("%d cells from the current event", fAddedCells)); |
780 | } | |
781 | ||
ffe32451 | 782 | //________________________________________________________________________ |
783 | void AliJetModelBaseTask::CopyClusters() | |
784 | { | |
785 | // Copy all the clusters in the new collection | |
787a3c4f | 786 | if (!fClusters) |
787 | return; | |
ffe32451 | 788 | |
ffe32451 | 789 | const Int_t nClusters = fClusters->GetEntriesFast(); |
56cc432c | 790 | Int_t nCopiedClusters = 0; |
ffe32451 | 791 | |
b16bb001 | 792 | if (fEsdMode) { |
ffe32451 | 793 | for (Int_t i = 0; i < nClusters; ++i) { |
794 | AliESDCaloCluster *esdcluster = static_cast<AliESDCaloCluster*>(fClusters->At(i)); | |
56cc432c | 795 | if (!esdcluster || !esdcluster->IsEMCAL()) |
ffe32451 | 796 | continue; |
4358e58a | 797 | AliESDCaloCluster *clus = new ((*fOutClusters)[nCopiedClusters]) AliESDCaloCluster(*esdcluster); |
798 | if (!fIsMC) { | |
799 | TArrayI *labels = clus->GetLabelsArray(); | |
800 | if (labels) | |
801 | labels->Reset(); | |
802 | } | |
56cc432c | 803 | nCopiedClusters++; |
ffe32451 | 804 | } |
805 | } | |
806 | else { | |
807 | for (Int_t i = 0; i < nClusters; ++i) { | |
808 | AliAODCaloCluster *aodcluster = static_cast<AliAODCaloCluster*>(fClusters->At(i)); | |
56cc432c | 809 | if (!aodcluster || !aodcluster->IsEMCAL()) |
ffe32451 | 810 | continue; |
4358e58a | 811 | AliAODCaloCluster *clus = new ((*fOutClusters)[nCopiedClusters]) AliAODCaloCluster(*aodcluster); |
812 | if (!fIsMC) | |
813 | clus->SetLabel(0,0); | |
56cc432c | 814 | nCopiedClusters++; |
ffe32451 | 815 | } |
816 | } | |
817 | } | |
818 | ||
819 | //________________________________________________________________________ | |
820 | void AliJetModelBaseTask::CopyTracks() | |
821 | { | |
787a3c4f | 822 | // Copy all the tracks in the new collection |
823 | ||
824 | if (!fTracks) | |
825 | return; | |
826 | ||
ffe32451 | 827 | const Int_t nTracks = fTracks->GetEntriesFast(); |
56cc432c | 828 | Int_t nCopiedTracks = 0; |
ffe32451 | 829 | for (Int_t i = 0; i < nTracks; ++i) { |
4358e58a | 830 | AliPicoTrack *picotrack = static_cast<AliPicoTrack*>(fTracks->At(i)); |
831 | if (!picotrack) | |
ffe32451 | 832 | continue; |
4358e58a | 833 | AliPicoTrack *track = new ((*fOutTracks)[nCopiedTracks]) AliPicoTrack(*picotrack); |
834 | if (!fIsMC && track->GetLabel() != 0) | |
835 | track->SetLabel(0); | |
56cc432c | 836 | nCopiedTracks++; |
ffe32451 | 837 | } |
838 | } | |
839 | ||
787a3c4f | 840 | //________________________________________________________________________ |
841 | void AliJetModelBaseTask::CopyMCParticles() | |
842 | { | |
843 | // Copy all the MC particles in the new collection | |
844 | ||
845 | if (!fMCParticles) | |
846 | return; | |
847 | ||
848 | const Int_t nPart = fMCParticles->GetEntriesFast(); | |
849 | Int_t nCopiedPart = 0; | |
850 | for (Int_t i = 0; i < nPart; ++i) { | |
851 | AliAODMCParticle *part = static_cast<AliAODMCParticle*>(fMCParticles->At(i)); | |
852 | if (!part) | |
853 | continue; | |
854 | new ((*fOutMCParticles)[nCopiedPart]) AliAODMCParticle(*part); | |
855 | ||
856 | nCopiedPart++; | |
857 | } | |
858 | ||
a7477843 | 859 | if (!fMCParticlesMap || !fOutMCParticlesMap) |
787a3c4f | 860 | return; |
861 | ||
a7477843 | 862 | if (fOutMCParticlesMap->GetSize() < fMCParticlesMap->GetSize()) |
863 | fOutMCParticlesMap->Set(fMCParticlesMap->GetSize() * 2); | |
864 | ||
5be3857d | 865 | for (Int_t i = 0; i < fMCParticlesMap->GetSize(); i++) { |
866 | fOutMCParticlesMap->AddAt(fMCParticlesMap->At(i), i); | |
867 | if (fMCParticlesMap->At(i) >= 0) | |
5ce8ae64 | 868 | fMCLabelShift = i; |
787a3c4f | 869 | } |
870 | ||
5ce8ae64 | 871 | AliDebug(2,Form("MC particles copied. fMCLabelShift=%d",fMCLabelShift)); |
787a3c4f | 872 | } |
873 | ||
762e8424 | 874 | //________________________________________________________________________ |
875 | void AliJetModelBaseTask::GetRandomCell(Double_t &eta, Double_t &phi, Int_t &absId) | |
876 | { | |
acd859a2 | 877 | // Get random cell. |
878 | ||
762e8424 | 879 | Int_t repeats = 0; |
05debff6 | 880 | Double_t rndEta = eta; |
881 | Double_t rndPhi = phi; | |
762e8424 | 882 | do { |
05debff6 | 883 | if (eta < -100) |
2103dc6a | 884 | rndEta = GetRandomEta(kTRUE); |
05debff6 | 885 | if (phi < 0) |
2103dc6a | 886 | rndPhi = GetRandomPhi(kTRUE); |
05debff6 | 887 | fGeom->GetAbsCellIdFromEtaPhi(rndEta, rndPhi, absId); |
762e8424 | 888 | repeats++; |
889 | } while (absId == -1 && repeats < 100); | |
890 | ||
891 | if (!(absId > -1)) { | |
892 | AliWarning(Form("Could not extract random cluster! Random eta-phi extracted more than 100 times!\n" | |
893 | "eta [%f, %f], phi [%f, %f]\n", fEtaMin, fEtaMax, fPhiMin, fPhiMax)); | |
894 | } | |
05debff6 | 895 | else { |
896 | eta = rndEta; | |
897 | phi = rndPhi; | |
898 | } | |
762e8424 | 899 | } |
900 | ||
901 | //________________________________________________________________________ | |
2103dc6a | 902 | Double_t AliJetModelBaseTask::GetRandomEta(Bool_t emcal) |
762e8424 | 903 | { |
acd859a2 | 904 | // Get random eta. |
905 | ||
2103dc6a | 906 | Double_t etamax = fEtaMax; |
907 | Double_t etamin = fEtaMin; | |
908 | ||
909 | if (emcal) { | |
910 | const Double_t EmcalMinEta = fGeom->GetArm1EtaMin(); | |
911 | const Double_t EmcalMaxEta = fGeom->GetArm1EtaMax(); | |
912 | ||
913 | if (etamax > EmcalMaxEta) etamax = EmcalMaxEta; | |
914 | if (etamax < EmcalMinEta) etamax = EmcalMinEta; | |
915 | if (etamin > EmcalMaxEta) etamin = EmcalMaxEta; | |
916 | if (etamin < EmcalMinEta) etamin = EmcalMinEta; | |
917 | } | |
918 | ||
919 | return gRandom->Rndm() * (etamax - etamin) + etamin; | |
762e8424 | 920 | } |
921 | ||
922 | //________________________________________________________________________ | |
2103dc6a | 923 | Double_t AliJetModelBaseTask::GetRandomPhi(Bool_t emcal) |
762e8424 | 924 | { |
acd859a2 | 925 | // Get random phi. |
2103dc6a | 926 | |
927 | Double_t phimax = fPhiMax; | |
928 | Double_t phimin = fPhiMin; | |
929 | ||
930 | if (emcal) { | |
931 | const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad(); | |
932 | const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad(); | |
933 | ||
934 | if (phimax > EmcalMaxPhi) phimax = EmcalMaxPhi; | |
935 | if (phimax < EmcalMinPhi) phimax = EmcalMinPhi; | |
936 | if (phimin > EmcalMaxPhi) phimin = EmcalMaxPhi; | |
937 | if (phimin < EmcalMinPhi) phimin = EmcalMinPhi; | |
938 | } | |
acd859a2 | 939 | |
43032ce2 | 940 | Double_t result = gRandom->Rndm() * (phimax - phimin) + phimin; |
941 | ||
942 | return result; | |
762e8424 | 943 | } |
944 | ||
945 | //________________________________________________________________________ | |
946 | Double_t AliJetModelBaseTask::GetRandomPt() | |
947 | { | |
acd859a2 | 948 | // Get random pt. |
949 | ||
e44e8726 | 950 | if (fPtSpectrum) |
951 | return fPtSpectrum->GetRandom(); | |
952 | else | |
953 | return gRandom->Rndm() * (fPtMax - fPtMin) + fPtMin; | |
762e8424 | 954 | } |
955 | ||
43032ce2 | 956 | //________________________________________________________________________ |
957 | void AliJetModelBaseTask::GetRandomParticle(Double_t &pt, Double_t &eta, Double_t &phi, Bool_t emcal) | |
958 | { | |
959 | // Get a random particle. | |
960 | ||
961 | eta = GetRandomEta(emcal); | |
962 | ||
963 | if (fPtPhiEvPlDistribution) { | |
964 | Double_t phimax = fPhiMax; | |
965 | Double_t phimin = fPhiMin; | |
966 | ||
967 | if (emcal) { | |
968 | const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad(); | |
969 | const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad(); | |
970 | ||
971 | if (phimax > EmcalMaxPhi) phimax = EmcalMaxPhi; | |
972 | if (phimax < EmcalMinPhi) phimax = EmcalMinPhi; | |
973 | if (phimin > EmcalMaxPhi) phimin = EmcalMaxPhi; | |
974 | if (phimin < EmcalMinPhi) phimin = EmcalMinPhi; | |
975 | } | |
976 | ||
977 | if (fPtPhiEvPlDistribution->GetXmin() > phimax || fPtPhiEvPlDistribution->GetXmax() < phimin) { | |
978 | AliWarning(Form("The hisogram %s does not overlap with the EMCal acceptance limits. It will be ignored.",fPtPhiEvPlDistribution->GetName())); | |
979 | pt = GetRandomPt(); | |
980 | phi = GetRandomPhi(emcal); | |
981 | } | |
982 | else { | |
983 | do { | |
984 | fPtPhiEvPlDistribution->GetRandom2(pt,phi); | |
985 | phi += fPsi; | |
986 | if (phi > TMath::Pi() * 2) phi -= TMath::Pi() * 2; | |
987 | } while (phi > phimax || phi < phimin); | |
988 | } | |
989 | } | |
990 | else { | |
991 | pt = GetRandomPt(); | |
992 | phi = GetRandomPhi(emcal); | |
993 | } | |
994 | } | |
995 | ||
762e8424 | 996 | //________________________________________________________________________ |
997 | void AliJetModelBaseTask::Run() | |
998 | { | |
acd859a2 | 999 | // Run. |
762e8424 | 1000 | } |