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