]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/UserTasks/EmcalTasks/AliAnalysisTaskEMCALPi0PbPb.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Protection for tracks with negative ID
[u/mrichter/AliRoot.git] / PWG4 / UserTasks / EmcalTasks / AliAnalysisTaskEMCALPi0PbPb.cxx
CommitLineData
ea3fd2d5 1// $Id$
2
ea3fd2d5 3#include "AliAnalysisTaskEMCALPi0PbPb.h"
fa443410 4#include <TAxis.h>
717fe7de 5#include <TChain.h>
6#include <TClonesArray.h>
f5d4ab70 7#include <TFile.h>
296ea9b4 8#include <TGeoGlobalMagField.h>
717fe7de 9#include <TH1F.h>
10#include <TH2F.h>
11#include <TList.h>
12#include <TLorentzVector.h>
f5d4ab70 13#include <TNtuple.h>
296ea9b4 14#include <TProfile.h>
0ceca1f8 15#include <TRegexp.h>
b3ee6797 16#include <TString.h>
296ea9b4 17#include <TVector2.h>
717fe7de 18#include "AliAODEvent.h"
56fd6cb2 19#include "AliAODMCParticle.h"
717fe7de 20#include "AliAODVertex.h"
21#include "AliAnalysisManager.h"
22#include "AliAnalysisTaskEMCALClusterizeFast.h"
296ea9b4 23#include "AliCDBManager.h"
ea3fd2d5 24#include "AliCentrality.h"
c2fe5f0e 25#include "AliEMCALDigit.h"
2ef5608f 26#include "AliEMCALGeometry.h"
0fbe8d4f 27#include "AliEMCALRecPoint.h"
296ea9b4 28#include "AliEMCALRecoUtils.h"
3a952328 29#include "AliESDCaloTrigger.h"
717fe7de 30#include "AliESDEvent.h"
5fe1ca23 31#include "AliESDUtils.h"
717fe7de 32#include "AliESDVertex.h"
296ea9b4 33#include "AliESDtrack.h"
0ec74551 34#include "AliESDtrackCuts.h"
b6c599fe 35#include "AliEventplane.h"
296ea9b4 36#include "AliGeomManager.h"
b3ee6797 37#include "AliInputEventHandler.h"
43cfaa06 38#include "AliLog.h"
38727e64 39#include "AliMCEvent.h"
40#include "AliMCParticle.h"
296ea9b4 41#include "AliMagF.h"
5fe1ca23 42#include "AliMultiplicity.h"
38727e64 43#include "AliStack.h"
0ec74551 44#include "AliTrackerBase.h"
c0d2e1f2 45#include "AliTriggerAnalysis.h"
ea3fd2d5 46
47ClassImp(AliAnalysisTaskEMCALPi0PbPb)
717fe7de 48
4ea96211 49//________________________________________________________________________
50AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb()
51 : AliAnalysisTaskSE(),
52 fCentVar("V0M"),
53 fCentFrom(0),
54 fCentTo(100),
55 fVtxZMin(-10),
56 fVtxZMax(+10),
57 fUseQualFlag(1),
58 fClusName(),
59 fDoNtuple(0),
60 fDoAfterburner(0),
61 fAsymMax(1),
62 fNminCells(1),
63 fMinE(0.100),
64 fMinErat(0),
65 fMinEcc(0),
66 fGeoName("EMCAL_FIRSTYEARV1"),
67 fMinNClusPerTr(50),
68 fIsoDist(0.2),
69 fTrClassNames(""),
70 fTrCuts(0),
71 fPrimTrCuts(0),
e0e1022c 72 fPrimTracksName(""),
4ea96211 73 fDoTrMatGeom(0),
74 fTrainMode(0),
75 fMarkCells(),
76 fMinL0Time(-1),
77 fMaxL0Time(1024),
78 fMcMode(0),
79 fEmbedMode(0),
80 fGeom(0),
81 fReco(0),
469b2bff 82 fTrigName(),
4ea96211 83 fDoPSel(kTRUE),
84 fIsGeoMatsSet(0),
85 fNEvs(0),
86 fOutput(0),
87 fTrClassNamesArr(0),
88 fEsdEv(0),
89 fAodEv(0),
90 fRecPoints(0),
91 fDigits(0),
92 fEsdClusters(0),
93 fEsdCells(0),
94 fAodClusters(0),
95 fAodCells(0),
96 fPtRanges(0),
97 fSelTracks(0),
98 fSelPrimTracks(0),
4ea96211 99 fNtuple(0),
100 fHeader(0),
101 fPrimVert(0),
102 fSpdVert(0),
103 fTpcVert(0),
104 fClusters(0),
105 fTriggers(0),
106 fMcParts(0),
107 fHCuts(0x0),
108 fHVertexZ(0x0),
109 fHVertexZ2(0x0),
110 fHCent(0x0),
111 fHCentQual(0x0),
112 fHTclsBeforeCuts(0x0),
113 fHTclsAfterCuts(0x0),
114 fHColuRow(0x0),
115 fHColuRowE(0x0),
116 fHCellMult(0x0),
117 fHCellE(0x0),
118 fHCellH(0x0),
119 fHCellM(0x0),
120 fHCellM2(0x0),
121 fHCellFreqNoCut(0x0),
122 fHCellFreqCut100M(0x0),
123 fHCellFreqCut300M(0x0),
124 fHCellFreqE(0x0),
125 fHCellCheckE(0x0),
126 fHClustEccentricity(0),
127 fHClustEtaPhi(0x0),
128 fHClustEnergyPt(0x0),
129 fHClustEnergySigma(0x0),
130 fHClustSigmaSigma(0x0),
131 fHClustNCellEnergyRatio(0x0),
132 fHClustEnergyNCell(0x0),
133 fHPrimTrackPt(0x0),
134 fHPrimTrackEta(0x0),
135 fHPrimTrackPhi(0x0),
136 fHMatchDr(0x0),
137 fHMatchDz(0x0),
138 fHMatchEp(0x0),
139 fHPionEtaPhi(0x0),
140 fHPionMggPt(0x0),
141 fHPionMggAsym(0x0),
142 fHPionMggDgg(0x0)
143{
144 // Constructor.
145}
146
ea3fd2d5 147//________________________________________________________________________
148AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name)
149 : AliAnalysisTaskSE(name),
717fe7de 150 fCentVar("V0M"),
151 fCentFrom(0),
152 fCentTo(100),
76332037 153 fVtxZMin(-10),
154 fVtxZMax(+10),
43cfaa06 155 fUseQualFlag(1),
717fe7de 156 fClusName(),
f5d4ab70 157 fDoNtuple(0),
a49742b5 158 fDoAfterburner(0),
f224d35b 159 fAsymMax(1),
a49742b5 160 fNminCells(1),
3c661da5 161 fMinE(0.100),
f224d35b 162 fMinErat(0),
163 fMinEcc(0),
6bf90832 164 fGeoName("EMCAL_FIRSTYEARV1"),
b6c599fe 165 fMinNClusPerTr(50),
296ea9b4 166 fIsoDist(0.2),
b3ee6797 167 fTrClassNames(""),
0ec74551 168 fTrCuts(0),
b6c599fe 169 fPrimTrCuts(0),
e0e1022c 170 fPrimTracksName(""),
b6c599fe 171 fDoTrMatGeom(0),
3a952328 172 fTrainMode(0),
173 fMarkCells(),
174 fMinL0Time(-1),
175 fMaxL0Time(1024),
38727e64 176 fMcMode(0),
cfd7d5b2 177 fEmbedMode(0),
d595acbb 178 fGeom(0),
296ea9b4 179 fReco(0),
469b2bff 180 fTrigName(),
807016ea 181 fDoPSel(kTRUE),
38727e64 182 fIsGeoMatsSet(0),
183 fNEvs(0),
717fe7de 184 fOutput(0),
b3ee6797 185 fTrClassNamesArr(0),
717fe7de 186 fEsdEv(0),
187 fAodEv(0),
43cfaa06 188 fRecPoints(0),
c2fe5f0e 189 fDigits(0),
717fe7de 190 fEsdClusters(0),
191 fEsdCells(0),
192 fAodClusters(0),
286b47a5 193 fAodCells(0),
fa443410 194 fPtRanges(0),
296ea9b4 195 fSelTracks(0),
b6c599fe 196 fSelPrimTracks(0),
788ca675 197 fNtuple(0),
198 fHeader(0),
199 fPrimVert(0),
200 fSpdVert(0),
201 fTpcVert(0),
202 fClusters(0),
3a952328 203 fTriggers(0),
807016ea 204 fMcParts(0),
fa443410 205 fHCuts(0x0),
206 fHVertexZ(0x0),
76332037 207 fHVertexZ2(0x0),
fa443410 208 fHCent(0x0),
209 fHCentQual(0x0),
b3ee6797 210 fHTclsBeforeCuts(0x0),
211 fHTclsAfterCuts(0x0),
d595acbb 212 fHColuRow(0x0),
213 fHColuRowE(0x0),
214 fHCellMult(0x0),
215 fHCellE(0x0),
216 fHCellH(0x0),
6eb6260e 217 fHCellM(0x0),
a49742b5 218 fHCellM2(0x0),
296ea9b4 219 fHCellFreqNoCut(0x0),
2e4d8148 220 fHCellFreqCut100M(0x0),
221 fHCellFreqCut300M(0x0),
222 fHCellFreqE(0x0),
296ea9b4 223 fHCellCheckE(0x0),
6eb6260e 224 fHClustEccentricity(0),
fa443410 225 fHClustEtaPhi(0x0),
226 fHClustEnergyPt(0x0),
227 fHClustEnergySigma(0x0),
d595acbb 228 fHClustSigmaSigma(0x0),
6eb6260e 229 fHClustNCellEnergyRatio(0x0),
f5e0f1e2 230 fHClustEnergyNCell(0x0),
231 fHPrimTrackPt(0x0),
232 fHPrimTrackEta(0x0),
233 fHPrimTrackPhi(0x0),
b6c599fe 234 fHMatchDr(0x0),
235 fHMatchDz(0x0),
236 fHMatchEp(0x0),
fa443410 237 fHPionEtaPhi(0x0),
238 fHPionMggPt(0x0),
6eb6260e 239 fHPionMggAsym(0x0),
240 fHPionMggDgg(0x0)
ea3fd2d5 241{
717fe7de 242 // Constructor.
ea3fd2d5 243
ea3fd2d5 244 DefineOutput(1, TList::Class());
469b2bff 245 fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,EMCALCells.,Tracks,EMCALTrigger.,SPDPileupVertices,TrkPileupVertices "
296ea9b4 246 "AOD:header,vertices,emcalCells,tracks";
ea3fd2d5 247}
ea3fd2d5 248
717fe7de 249//________________________________________________________________________
250AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb()
251{
252 // Destructor.
ea3fd2d5 253
b3ee6797 254 if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
255 delete fOutput; fOutput = 0;
256 }
fa443410 257 delete fPtRanges; fPtRanges = 0;
a2de5ca1 258 fGeom = 0; // do not delete geometry when using instance
296ea9b4 259 delete fReco; fReco = 0;
b3ee6797 260 delete fTrClassNamesArr;
296ea9b4 261 delete fSelTracks;
b6c599fe 262 delete fSelPrimTracks;
d595acbb 263 delete [] fHColuRow;
264 delete [] fHColuRowE;
265 delete [] fHCellMult;
296ea9b4 266 delete [] fHCellFreqNoCut;
2e4d8148 267 delete [] fHCellFreqCut100M;
268 delete [] fHCellFreqCut300M;
ea3fd2d5 269}
717fe7de 270
ea3fd2d5 271//________________________________________________________________________
272void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects()
273{
717fe7de 274 // Create user objects here.
ea3fd2d5 275
f3582e89 276 cout << "AliAnalysisTaskEMCALPi0PbPb: Input settings" << endl;
277 cout << " fCentVar: " << fCentVar << endl;
278 cout << " fCentFrom: " << fCentFrom << endl;
279 cout << " fCentTo: " << fCentTo << endl;
280 cout << " fVtxZMin: " << fVtxZMin << endl;
281 cout << " fVtxZMax: " << fVtxZMax << endl;
282 cout << " fUseQualFlag: " << fUseQualFlag << endl;
283 cout << " fClusName: \"" << fClusName << "\"" << endl;
284 cout << " fDoNtuple: " << fDoNtuple << endl;
285 cout << " fDoAfterburner: " << fDoAfterburner << endl;
286 cout << " fAsymMax: " << fAsymMax << endl;
287 cout << " fNminCells: " << fNminCells << endl;
288 cout << " fMinE: " << fMinE << endl;
289 cout << " fMinErat: " << fMinErat << endl;
290 cout << " fMinEcc: " << fMinEcc << endl;
291 cout << " fGeoName: \"" << fGeoName << "\"" << endl;
292 cout << " fMinNClusPerTr: " << fMinNClusPerTr << endl;
38727e64 293 cout << " fIsoDist: " << fIsoDist << endl;
f3582e89 294 cout << " fTrClassNames: \"" << fTrClassNames << "\"" << endl;
f3582e89 295 cout << " fTrCuts: " << fTrCuts << endl;
296 cout << " fPrimTrCuts: " << fPrimTrCuts << endl;
297 cout << " fDoTrMatGeom: " << fDoTrMatGeom << endl;
298 cout << " fTrainMode: " << fTrainMode << endl;
299 cout << " fMarkCells: " << fMarkCells << endl;
300 cout << " fMinL0Time: " << fMinL0Time << endl;
301 cout << " fMaxL0Time: " << fMaxL0Time << endl;
38727e64 302 cout << " fMcMode: " << fMcMode << endl;
cfd7d5b2 303 cout << " fEmbedMode: " << fEmbedMode << endl;
38727e64 304 cout << " fGeom: " << fGeom << endl;
305 cout << " fReco: " << fReco << endl;
2ee7bda4 306 cout << " fTrigName: " << fTrigName << endl;
807016ea 307 cout << " fDoPSel: " << fDoPSel << endl;
38727e64 308
309 if (!fGeom)
80d69b70 310 fGeom = AliEMCALGeometry::GetInstance(fGeoName);
38727e64 311 else {
312 if (fGeom->GetMatrixForSuperModule(0))
313 fIsGeoMatsSet = kTRUE;
314 }
315 if (!fReco)
316 fReco = new AliEMCALRecoUtils();
b3ee6797 317 fTrClassNamesArr = fTrClassNames.Tokenize(" ");
717fe7de 318 fOutput = new TList();
4ea96211 319 fOutput->SetOwner(1);
296ea9b4 320 fSelTracks = new TObjArray;
b6c599fe 321 fSelPrimTracks = new TObjArray;
ea3fd2d5 322
f5d4ab70 323 if (fDoNtuple) {
a9c1389d 324 TFile *f = OpenFile(1);
325 TDirectory::TContext context(f);
6bf90832 326 if (f) {
327 f->SetCompressionLevel(2);
788ca675 328 fNtuple = new TTree(Form("tree%.0fto%.0f",fCentFrom,fCentTo), "StandaloneTree");
6bf90832 329 fNtuple->SetDirectory(f);
3a952328 330 if (fTrainMode) {
331 fNtuple->SetAutoFlush(-2*1024*1024);
332 fNtuple->SetAutoSave(0);
333 } else {
334 fNtuple->SetAutoFlush(-32*1024*1024);
335 fNtuple->SetAutoSave(0);
336 }
337
788ca675 338 fHeader = new AliStaHeader;
339 fNtuple->Branch("header", &fHeader, 16*1024, 99);
340 fPrimVert = new AliStaVertex;
341 fNtuple->Branch("primv", &fPrimVert, 16*1024, 99);
342 fSpdVert = new AliStaVertex;
343 fNtuple->Branch("spdv", &fSpdVert, 16*1024, 99);
344 fTpcVert = new AliStaVertex;
345 fNtuple->Branch("tpcv", &fTpcVert, 16*1024, 99);
346 if (TClass::GetClass("AliStaCluster"))
347 TClass::GetClass("AliStaCluster")->IgnoreTObjectStreamer();
3a952328 348 fClusters = new TClonesArray("AliStaCluster");
788ca675 349 fNtuple->Branch("clusters", &fClusters, 8*16*1024, 99);
3a952328 350 if (TClass::GetClass("AliStaTrigger"))
351 TClass::GetClass("AliStaTrigger")->IgnoreTObjectStreamer();
352 fTriggers = new TClonesArray("AliStaTrigger");
353 fNtuple->Branch("l0prim", &fTriggers, 16*1024, 99);
cfd7d5b2 354 if (fMcMode||fEmbedMode) {
807016ea 355 if (TClass::GetClass("AliStaPart"))
356 TClass::GetClass("AliStaPart")->IgnoreTObjectStreamer();
357 fMcParts = new TClonesArray("AliStaPart");
358 fNtuple->Branch("mcparts", &fMcParts, 8*16*1024, 99);
359 }
6bf90832 360 }
f5d4ab70 361 }
362
002dcebe 363 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
b6c599fe 364 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
365 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
002dcebe 366
d595acbb 367 // histograms
587ce6c6 368 Bool_t th1 = TH1::GetDefaultSumw2();
a49742b5 369 TH1::SetDefaultSumw2(kTRUE);
587ce6c6 370 Bool_t th2 = TH2::GetDefaultSumw2();
a49742b5 371 TH2::SetDefaultSumw2(kTRUE);
b3ee6797 372 fHCuts = new TH1F("hEventCuts","",5,0.5,5.5);
373 fHCuts->GetXaxis()->SetBinLabel(1,"All");
374 fHCuts->GetXaxis()->SetBinLabel(2,"PS");
375 fHCuts->GetXaxis()->SetBinLabel(3,Form("%s: %.0f-%.0f",fCentVar.Data(),fCentFrom,fCentTo));
376 fHCuts->GetXaxis()->SetBinLabel(4,"QFlag");
377 fHCuts->GetXaxis()->SetBinLabel(5,Form("zvtx: %.0f-%.0f",fVtxZMin,fVtxZMax));
fa443410 378 fOutput->Add(fHCuts);
379 fHVertexZ = new TH1F("hVertexZBeforeCut","",100,-25,25);
380 fHVertexZ->SetXTitle("z [cm]");
381 fOutput->Add(fHVertexZ);
76332037 382 fHVertexZ2 = new TH1F("hVertexZAfterCut","",100,-25,25);
383 fHVertexZ2->SetXTitle("z [cm]");
384 fOutput->Add(fHVertexZ2);
f2b8fca6 385 fHCent = new TH1F("hCentBeforeCut","",102,-1,101);
fa443410 386 fHCent->SetXTitle(fCentVar.Data());
76332037 387 fOutput->Add(fHCent);
f2b8fca6 388 fHCentQual = new TH1F("hCentAfterCut","",102,-1,101);
fa443410 389 fHCentQual->SetXTitle(fCentVar.Data());
390 fOutput->Add(fHCentQual);
2e4d8148 391 fHTclsBeforeCuts = new TH1F("hTclsBeforeCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
392 fHTclsAfterCuts = new TH1F("hTclsAfterCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
b3ee6797 393 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
394 const char *name = fTrClassNamesArr->At(i)->GetName();
395 fHTclsBeforeCuts->GetXaxis()->SetBinLabel(1+i,name);
396 fHTclsAfterCuts->GetXaxis()->SetBinLabel(1+i,name);
397 }
398 fOutput->Add(fHTclsBeforeCuts);
399 fOutput->Add(fHTclsAfterCuts);
90d5b88b 400
d595acbb 401 // histograms for cells
296ea9b4 402 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
403 fHColuRow = new TH2*[nsm];
404 fHColuRowE = new TH2*[nsm];
405 fHCellMult = new TH1*[nsm];
406 for (Int_t i = 0; i<nsm; ++i) {
2e4d8148 407 fHColuRow[i] = new TH2F(Form("hColRow_Mod%d",i),"",48,0,48,24,0.,24);
d595acbb 408 fHColuRow[i]->SetTitle(Form("Module %d: Occupancy", i));
409 fHColuRow[i]->SetXTitle("col (i#eta)");
410 fHColuRow[i]->SetYTitle("row (i#phi)");
2e4d8148 411 fHColuRowE[i] = new TH2F(Form("hColRowE_Mod%d", i),"",48,0,48,24,0,24);
90d5b88b 412 fHColuRowE[i]->SetTitle(Form("Module %d: Cell energy",i));
d595acbb 413 fHColuRowE[i]->SetXTitle("col (i#eta)");
414 fHColuRowE[i]->SetYTitle("row (i#phi)");
415 fHCellMult[i] = new TH1F(Form("hCellMult_Mod%d",i),"",1000,0,1000);
416 fHCellMult[i]->SetTitle(Form("Module %d: Cell multiplicity",i));
90d5b88b 417 fHCellMult[i]->SetXTitle("# of cells");
d595acbb 418 fOutput->Add(fHColuRow[i]);
419 fOutput->Add(fHColuRowE[i]);
420 fOutput->Add(fHCellMult[i]);
421 }
2e4d8148 422 fHCellE = new TH1F("hCellE","",250,0.,25.);
d595acbb 423 fHCellE->SetXTitle("E_{cell} [GeV]");
424 fOutput->Add(fHCellE);
fb5a0b69 425 fHCellH = new TH1F ("hCellHighestE","",250,0.,25.);
6eb6260e 426 fHCellH->SetXTitle("E^{max}_{cell} [GeV]");
d595acbb 427 fOutput->Add(fHCellH);
fb5a0b69 428 fHCellM = new TH1F ("hCellMeanEperHitCell","",250,0.,2.5);
6eb6260e 429 fHCellM->SetXTitle("#LT E_{cell}#GT [GeV]");
430 fOutput->Add(fHCellM);
fb5a0b69 431 fHCellM2 = new TH1F ("hCellMeanEperAllCells","",250,0.,1);
f4ec884e 432 fHCellM2->SetXTitle("1/N_{cells} #Sigma E_{cell} [GeV]");
a49742b5 433 fOutput->Add(fHCellM2);
90d5b88b 434
2e4d8148 435 fHCellFreqNoCut = new TH1*[nsm];
436 fHCellFreqCut100M = new TH1*[nsm];
437 fHCellFreqCut300M = new TH1*[nsm];
438 fHCellFreqE = new TH1*[nsm];
296ea9b4 439 for (Int_t i = 0; i<nsm; ++i){
440 Double_t lbin = i*24*48-0.5;
441 Double_t hbin = lbin+24*48;
2e4d8148 442 fHCellFreqNoCut[i] = new TH1F(Form("hCellFreqNoCut_SM%d",i),
443 Form("Frequency SM%d (no cut);id;#",i), 1152, lbin, hbin);
444 fHCellFreqCut100M[i] = new TH1F(Form("hCellFreqCut100M_SM%d",i),
445 Form("Frequency SM%d (>0.1GeV);id;#",i), 1152, lbin, hbin);
446 fHCellFreqCut300M[i] = new TH1F(Form("hCellFreqCut300M_SM%d",i),
447 Form("Frequency SM%d (>0.3GeV);id;#",i), 1152, lbin, hbin);
448 fHCellFreqE[i] = new TH1F(Form("hCellFreqE_SM%d",i),
449 Form("Frequency SM%d (E weighted);id;#",i), 1152, lbin, hbin);
296ea9b4 450 fOutput->Add(fHCellFreqNoCut[i]);
2e4d8148 451 fOutput->Add(fHCellFreqCut100M[i]);
452 fOutput->Add(fHCellFreqCut300M[i]);
453 fOutput->Add(fHCellFreqE[i]);
296ea9b4 454 }
3a952328 455 if (!fMarkCells.IsNull()) {
296ea9b4 456 fHCellCheckE = new TH1*[24*48*nsm];
408dc04c 457 memset(fHCellCheckE,0,24*48*nsm*sizeof(TH1*));
3a952328 458 TObjArray *cells = fMarkCells.Tokenize(" ");
459 Int_t n = cells->GetEntries();
460 Int_t *tcs = new Int_t[n];
461 for (Int_t i=0;i<n;++i) {
462 TString name(cells->At(i)->GetName());
463 tcs[i]=name.Atoi();
464 }
465 for (Int_t i = 0; i<n; ++i) {
296ea9b4 466 Int_t c=tcs[i];
467 if (c<24*48*nsm) {
3a952328 468 fHCellCheckE[i] = new TH1F(Form("hCellE_id%d",c), Form("Cell %d;E [GeV/c];#",c), 1000, 0, 10);
296ea9b4 469 fOutput->Add(fHCellCheckE[i]);
470 }
471 }
3a952328 472 delete cells;
473 delete [] tcs;
296ea9b4 474 }
475
d595acbb 476 // histograms for clusters
3a952328 477 if (!fTrainMode) {
478 fHClustEccentricity = new TH1F("hClustEccentricity","",100,-0.1,1.1);
479 fHClustEccentricity->SetXTitle("#epsilon_{C}");
480 fOutput->Add(fHClustEccentricity);
481 fHClustEtaPhi = new TH2F("hClustEtaPhi","",500,-0.8,0.8,100*nsm,phimin,phimax);
482 fHClustEtaPhi->SetXTitle("#eta");
483 fHClustEtaPhi->SetYTitle("#varphi");
484 fOutput->Add(fHClustEtaPhi);
485 fHClustEnergyPt = new TH2F("hClustEnergyPt","",250,0,50,250,0,50);
486 fHClustEnergyPt->SetXTitle("E [GeV]");
487 fHClustEnergyPt->SetYTitle("p_{T} [GeV/c]");
488 fOutput->Add(fHClustEnergyPt);
489 fHClustEnergySigma = new TH2F("hClustEnergySigma","",100,0,100,500,0,50);
490 fHClustEnergySigma->SetXTitle("E_{C} * #sigma_{max} [GeV*cm]");
491 fHClustEnergySigma->SetYTitle("E_{C} [GeV]");
492 fOutput->Add(fHClustEnergySigma);
493 fHClustSigmaSigma = new TH2F("hClustSigmaSigma","",500,0,50,500,0,50);
494 fHClustSigmaSigma->SetXTitle("#lambda_{0} [cm]");
495 fHClustSigmaSigma->SetYTitle("#sigma_{max} [cm]");
496 fOutput->Add(fHClustSigmaSigma);
497 fHClustNCellEnergyRatio = new TH2F("hClustNCellEnergyRatio","",27,-0.5,26.5,101,-0.05,1.05);
498 fHClustNCellEnergyRatio->SetXTitle("N_{cells}");
499 fHClustNCellEnergyRatio->SetYTitle("E^{max}_{cell}/E_{clus}");
500 fOutput->Add(fHClustNCellEnergyRatio);
f5e0f1e2 501 fHClustEnergyNCell = new TH2F("hClustEnergyNCell","",200,0,100,50,0,50);
502 fHClustEnergyNCell->SetXTitle("E_{clus}");
503 fHClustEnergyNCell->SetYTitle("N_{cells}");
504 fOutput->Add(fHClustEnergyNCell);
3a952328 505 }
90d5b88b 506
f5e0f1e2 507 // histograms for primary tracks
508 fHPrimTrackPt = new TH1F("hPrimTrackPt",";p_{T} [GeV/c]",500,0,50);
509 fOutput->Add(fHPrimTrackPt);
510 fHPrimTrackEta = new TH1F("hPrimTrackEta",";#eta",40,-2,2);
511 fOutput->Add(fHPrimTrackEta);
512 fHPrimTrackPhi = new TH1F("hPrimTrackPhi",";#varPhi [rad]",63,0,6.3);
513 fOutput->Add(fHPrimTrackPhi);
a186e444 514
b6c599fe 515 // histograms for track matching
a186e444 516 if (fDoTrMatGeom) {
517 fHMatchDr = new TH1F("hMatchDrDist",";dR [cm]",500,0,200);
518 fOutput->Add(fHMatchDr);
519 fHMatchDz = new TH1F("hMatchDzDist",";dZ [cm]",500,-100,100);
520 fOutput->Add(fHMatchDz);
521 fHMatchEp = new TH1F("hMatchEpDist",";E/p",100,0,10);
522 fOutput->Add(fHMatchEp);
523 }
3a952328 524
d595acbb 525 // histograms for pion candidates
3a952328 526 if (!fTrainMode) {
527 fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100*nsm,phimin,phimax);
528 fHPionEtaPhi->SetXTitle("#eta_{#gamma#gamma}");
529 fHPionEtaPhi->SetYTitle("#varphi_{#gamma#gamma}");
530 fOutput->Add(fHPionEtaPhi);
531 fHPionMggPt = new TH2F("hPionMggPt","",1000,0,2,100,0,20.0);
532 fHPionMggPt->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
533 fHPionMggPt->SetYTitle("p_{T}^{#gamma#gamma} [GeV/c]");
534 fOutput->Add(fHPionMggPt);
535 fHPionMggAsym = new TH2F("hPionMggAsym","",1000,0,2,100,0,1);
536 fHPionMggAsym->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
537 fHPionMggAsym->SetYTitle("Z_{#gamma#gamma} [GeV]");
538 fOutput->Add(fHPionMggAsym);
539 fHPionMggDgg = new TH2F("hPionMggDgg","",1000,0,2,100,0,10);
540 fHPionMggDgg->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
541 fHPionMggDgg->SetYTitle("opening angle [grad]");
542 fOutput->Add(fHPionMggDgg);
543 const Int_t nbins = 20;
544 Double_t xbins[nbins] = {0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,6,7,8,9,10,12.5,15,20,25,50};
545 fPtRanges = new TAxis(nbins-1,xbins);
546 for (Int_t i = 0; i<=nbins; ++i) {
547 fHPionInvMasses[i] = new TH1F(Form("hPionInvMass%d",i),"",1000,0,2);
548 fHPionInvMasses[i]->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
549 if (i==0)
550 fHPionInvMasses[i]->SetTitle(Form("0 < p_{T}^{#gamma#gamma} <%.1f",xbins[0]));
551 else if (i==nbins)
552 fHPionInvMasses[i]->SetTitle(Form("p_{T}^{#gamma#gamma} > 50"));
553 else
554 fHPionInvMasses[i]->SetTitle(Form("%.1f < p_{T}^{#gamma#gamma} <%.1f",xbins[i-1],xbins[i]));
555 fOutput->Add(fHPionInvMasses[i]);
556 }
fa443410 557 }
296ea9b4 558
587ce6c6 559 TH1::SetDefaultSumw2(th1);
560 TH2::SetDefaultSumw2(th2);
717fe7de 561 PostData(1, fOutput);
ea3fd2d5 562}
563
564//________________________________________________________________________
565void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *)
566{
717fe7de 567 // Called for each event.
568
569 if (!InputEvent())
570 return;
571
43cfaa06 572 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
573 fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent());
44310bce 574 UInt_t offtrigger = 0;
43cfaa06 575 if (fEsdEv) {
576 am->LoadBranch("AliESDRun.");
577 am->LoadBranch("AliESDHeader.");
9809ed8c 578 UInt_t mask1 = fEsdEv->GetESDRun()->GetDetectorsInDAQ();
579 UInt_t mask2 = fEsdEv->GetESDRun()->GetDetectorsInReco();
580 Bool_t desc1 = (mask1 >> 18) & 0x1;
581 Bool_t desc2 = (mask2 >> 18) & 0x1;
2cee9afd 582 if (desc1==0 || desc2==0) { //AliDAQ::OfflineModuleName(18)=="EMCAL"
9809ed8c 583 AliError(Form("EMCAL not in DAQ/RECO: %u (%u)/%u (%u)",
584 mask1, fEsdEv->GetESDRun()->GetDetectorsInReco(),
585 mask2, fEsdEv->GetESDRun()->GetDetectorsInDAQ()));
44310bce 586 return;
587 }
588 offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
43cfaa06 589 } else {
590 fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
408dc04c 591 if (!fAodEv) {
592 AliFatal("Neither ESD nor AOD event found");
593 return;
594 }
43cfaa06 595 am->LoadBranch("header");
44310bce 596 offtrigger = fAodEv->GetHeader()->GetOfflineTrigger();
597 }
38727e64 598 if (!fMcMode && (offtrigger & AliVEvent::kFastOnly)) {
44310bce 599 AliWarning(Form("EMCAL not in fast only partition"));
600 return;
43cfaa06 601 }
b6c599fe 602 if (fDoTrMatGeom && !AliGeomManager::GetGeometry()) { // get geometry
27c2e3d9 603 AliWarning("Accessing geometry from OCDB, this is not very efficient!");
604 AliCDBManager *cdb = AliCDBManager::Instance();
605 if (!cdb->IsDefaultStorageSet())
606 cdb->SetDefaultStorage("raw://");
607 Int_t runno = InputEvent()->GetRunNumber();
608 if (runno != cdb->GetRun())
609 cdb->SetRun(runno);
610 AliGeomManager::LoadGeometry();
611 }
612
613 if (!AliGeomManager::GetGeometry()&&!fIsGeoMatsSet) { // set misalignment matrices (stored in first event)
614 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
9809ed8c 615 for (Int_t i=0; i<nsm; ++i) {
616 const TGeoHMatrix *geom = 0;
617 if (fEsdEv)
618 geom = fEsdEv->GetESDRun()->GetEMCALMatrix(i);
619 else
620 geom = fAodEv->GetHeader()->GetEMCALMatrix(i);
621 if (!geom)
622 continue;
623 geom->Print();
624 fGeom->SetMisalMatrix(geom,i);
27c2e3d9 625 }
626 fIsGeoMatsSet = kTRUE;
627 }
628
629 if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
630 if (fEsdEv) {
631 const AliESDRun *erun = fEsdEv->GetESDRun();
632 AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(),
633 erun->GetCurrentDip(),
634 AliMagF::kConvLHC,
635 kFALSE,
636 erun->GetBeamEnergy(),
637 erun->GetBeamType());
638 TGeoGlobalMagField::Instance()->SetField(field);
639 } else {
640 Double_t pol = -1; //polarity
641 Double_t be = -1; //beam energy
642 AliMagF::BeamType_t btype = AliMagF::kBeamTypepp;
643 Int_t runno = fAodEv->GetRunNumber();
644 if (runno>=136851 && runno<138275) {
645 pol = -1;
646 be = 2760;
647 btype = AliMagF::kBeamTypeAA;
648 } else if (runno>=138275 && runno<=139517) {
649 pol = +1;
650 be = 2760;
651 btype = AliMagF::kBeamTypeAA;
652 } else {
653 AliError(Form("Do not know the bfield parameters for run %d! Using defaults!!!", runno));
654 }
655 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", pol, pol, AliMagF::k5kG, btype, be));
656 }
657 }
658
b3ee6797 659 Int_t cut = 1;
660 fHCuts->Fill(cut++);
661
662 TString trgclasses;
663 AliESDHeader *h = dynamic_cast<AliESDHeader*>(InputEvent()->GetHeader());
664 if (h) {
665 trgclasses = fEsdEv->GetFiredTriggerClasses();
666 } else {
667 AliAODHeader *h2 = dynamic_cast<AliAODHeader*>(InputEvent()->GetHeader());
668 if (h2)
669 trgclasses = h2->GetFiredTriggerClasses();
670 }
671 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
672 const char *name = fTrClassNamesArr->At(i)->GetName();
0ceca1f8 673 TRegexp regexp(name);
674 if (trgclasses.Contains(regexp))
b3ee6797 675 fHTclsBeforeCuts->Fill(1+i);
676 }
677
807016ea 678 if (fDoPSel && offtrigger==0)
b3ee6797 679 return;
680
fa443410 681 fHCuts->Fill(cut++);
286b47a5 682
717fe7de 683 const AliCentrality *centP = InputEvent()->GetCentrality();
684 Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar);
fa443410 685 fHCent->Fill(cent);
717fe7de 686 if (cent<fCentFrom||cent>fCentTo)
286b47a5 687 return;
43cfaa06 688
fa443410 689 fHCuts->Fill(cut++);
286b47a5 690
43cfaa06 691 if (fUseQualFlag) {
692 if (centP->GetQuality()>0)
693 return;
694 }
286b47a5 695
fa443410 696 fHCentQual->Fill(cent);
697 fHCuts->Fill(cut++);
717fe7de 698
43cfaa06 699 if (fEsdEv) {
fa443410 700 am->LoadBranch("PrimaryVertex.");
701 am->LoadBranch("SPDVertex.");
702 am->LoadBranch("TPCVertex.");
43cfaa06 703 } else {
704 fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
705 am->LoadBranch("vertices");
3c661da5 706 if (!fAodEv) return;
43cfaa06 707 }
708
709 const AliVVertex *vertex = InputEvent()->GetPrimaryVertex();
717fe7de 710 if (!vertex)
711 return;
712
fa443410 713 fHVertexZ->Fill(vertex->GetZ());
286b47a5 714
717fe7de 715 if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax)
716 return;
286b47a5 717
fa443410 718 fHCuts->Fill(cut++);
76332037 719 fHVertexZ2->Fill(vertex->GetZ());
717fe7de 720
b3ee6797 721 // count number of accepted events
722 ++fNEvs;
723
724 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
725 const char *name = fTrClassNamesArr->At(i)->GetName();
0ceca1f8 726 TRegexp regexp(name);
727 if (trgclasses.Contains(regexp))
b3ee6797 728 fHTclsAfterCuts->Fill(1+i);
729 }
730
43cfaa06 731 fRecPoints = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
c2fe5f0e 732 fDigits = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
b6c599fe 733 fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set or if clusters are attached
43cfaa06 734 fEsdCells = 0; // will be set if ESD input used
b6c599fe 735 fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set or if clusters are attached
43cfaa06 736 // or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode
737 fAodCells = 0; // will be set if AOD input used
717fe7de 738
43cfaa06 739 // deal with special output from AliAnalysisTaskEMCALClusterizeFast first
2ee7bda4 740 Bool_t overwrite = 0;
43cfaa06 741 Bool_t clusattached = 0;
742 Bool_t recalibrated = 0;
743 if (1 && !fClusName.IsNull()) {
744 AliAnalysisTaskEMCALClusterizeFast *cltask = 0;
745 TObjArray *ts = am->GetTasks();
746 cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName));
747 if (cltask && cltask->GetClusters()) {
2ee7bda4 748 fRecPoints = cltask->GetClusters();
749 fDigits = cltask->GetDigits();
43cfaa06 750 clusattached = cltask->GetAttachClusters();
2ee7bda4 751 overwrite = cltask->GetOverwrite();
43cfaa06 752 if (cltask->GetCalibData()!=0)
753 recalibrated = kTRUE;
754 }
755 }
38727e64 756 if (1 && !fClusName.IsNull()) {
757 TList *l = 0;
758 if (AODEvent())
759 l = AODEvent()->GetList();
760 else if (fAodEv)
761 l = fAodEv->GetList();
717fe7de 762 if (l) {
38727e64 763 fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
717fe7de 764 }
ea3fd2d5 765 }
717fe7de 766
767 if (fEsdEv) { // ESD input mode
43cfaa06 768 if (1 && (!fRecPoints||clusattached)) {
2ee7bda4 769 if (!clusattached && !overwrite)
43cfaa06 770 am->LoadBranch("CaloClusters");
717fe7de 771 TList *l = fEsdEv->GetList();
2ee7bda4 772 if (clusattached) {
773 fEsdClusters = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
1829ebfe 774 } else {
38727e64 775 fEsdClusters = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
ea3fd2d5 776 }
777 }
43cfaa06 778 if (1) {
779 if (!recalibrated)
780 am->LoadBranch("EMCALCells.");
717fe7de 781 fEsdCells = fEsdEv->GetEMCALCells();
782 }
783 } else if (fAodEv) { // AOD input mode
43cfaa06 784 if (1 && (!fAodClusters || clusattached)) {
785 if (!clusattached)
786 am->LoadBranch("caloClusters");
717fe7de 787 TList *l = fAodEv->GetList();
717fe7de 788 if (l) {
38727e64 789 fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
ea3fd2d5 790 }
791 }
43cfaa06 792 if (1) {
793 if (!recalibrated)
794 am->LoadBranch("emcalCells");
717fe7de 795 fAodCells = fAodEv->GetEMCALCells();
796 }
797 } else {
798 AliFatal("Impossible to not have either pointer to ESD or AOD event");
799 }
ea3fd2d5 800
43cfaa06 801 if (1) {
802 AliDebug(2,Form("fRecPoints set: %p", fRecPoints));
c2fe5f0e 803 AliDebug(2,Form("fDigits set: %p", fDigits));
43cfaa06 804 AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters));
805 AliDebug(2,Form("fEsdCells set: %p", fEsdCells));
806 AliDebug(2,Form("fAodClusters set: %p", fAodClusters));
807 AliDebug(2,Form("fAodCells set: %p", fAodCells));
808 }
809
a49742b5 810 if (fDoAfterburner)
811 ClusterAfterburner();
6eb6260e 812
38727e64 813 CalcMcInfo();
3a952328 814 CalcCaloTriggers();
b6c599fe 815 CalcPrimTracks();
3a952328 816 CalcTracks();
296ea9b4 817 CalcClusterProps();
818
286b47a5 819 FillCellHists();
3a952328 820 if (!fTrainMode) {
821 FillClusHists();
822 FillPionHists();
823 FillOtherHists();
824 }
38727e64 825 FillMcHists();
788ca675 826 FillNtuple();
ea3fd2d5 827
3a952328 828 if (fTrainMode) {
829 fSelTracks->Clear();
830 fSelPrimTracks->Clear();
cfd7d5b2 831 if (fMcParts)
807016ea 832 fMcParts->Clear();
4ea96211 833 if (fTriggers)
834 fTriggers->Clear();
835 if (fClusters)
836 fClusters->Clear();
3a952328 837 }
de4cee41 838
717fe7de 839 PostData(1, fOutput);
ea3fd2d5 840}
841
842//________________________________________________________________________
843void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *)
844{
717fe7de 845 // Terminate called at the end of analysis.
f5d4ab70 846
847 if (fNtuple) {
4ea96211 848 TFile *f = OpenFile(1);
a9c1389d 849 TDirectory::TContext context(f);
f5d4ab70 850 if (f)
851 fNtuple->Write();
852 }
d9f26424 853
bc8a45bd 854 AliInfo(Form("%s: Accepted %lld events ", GetName(), fNEvs));
ea3fd2d5 855}
ea3fd2d5 856
296ea9b4 857//________________________________________________________________________
3a952328 858void AliAnalysisTaskEMCALPi0PbPb::CalcCaloTriggers()
296ea9b4 859{
3a952328 860 // Calculate triggers
296ea9b4 861
38727e64 862 if (fAodEv)
863 return; // information not available in AOD
864
4ea96211 865 if (!fTriggers)
866 return;
867
3a952328 868 fTriggers->Clear();
296ea9b4 869
2ee7bda4 870 if (fTrigName.Length()<=0)
3a952328 871 return;
872
2ee7bda4 873 TClonesArray *arr = dynamic_cast<TClonesArray*>(fEsdEv->FindListObject(fTrigName));
874 if (!arr) {
875 AliError(Form("Could not get array with name %s", fTrigName.Data()));
3a952328 876 return;
3a952328 877 }
878
2ee7bda4 879 Int_t nNumberOfCaloClusters = arr->GetEntries();
880 for(Int_t j = 0, ntrigs = 0; j < nNumberOfCaloClusters; ++j) {
881 AliVCluster *cl = dynamic_cast<AliVCluster*>(arr->At(j));
882 if (!cl)
3a952328 883 continue;
2ee7bda4 884 if (!cl->IsEMCAL())
3a952328 885 continue;
2ee7bda4 886 if (cl->E()<1)
c47674e0 887 continue;
3a952328 888 AliStaTrigger *trignew = static_cast<AliStaTrigger*>(fTriggers->New(ntrigs++));
2ee7bda4 889 Float_t pos[3] = {0,0,0};
890 cl->GetPosition(pos);
891 TVector3 vpos(pos);
892 trignew->fE = cl->E();
893 trignew->fEta = vpos.Eta();
894 trignew->fPhi = vpos.Phi();
895 Short_t id = -1;
896 GetMaxCellEnergy(cl, id);
897 trignew->fIdMax = id;
3a952328 898 }
899}
900
901//________________________________________________________________________
902void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
903{
904 // Calculate cluster properties
905
4ea96211 906 if (!fClusters)
907 return;
908
3a952328 909 fClusters->Clear();
910
911 AliVCaloCells *cells = fEsdCells;
912 if (!cells)
913 cells = fAodCells;
914 if (!cells)
915 return;
916
917 TObjArray *clusters = fEsdClusters;
918 if (!clusters)
919 clusters = fAodClusters;
920 if (!clusters)
921 return;
922
a186e444 923 const Int_t ncells = cells->GetNumberOfCells();
924 const Int_t nclus = clusters->GetEntries();
925 const Int_t ntrks = fSelTracks->GetEntries();
3a952328 926
927 std::map<Short_t,Short_t> map;
a186e444 928 for (Short_t pos=0; pos<ncells; ++pos) {
3a952328 929 Short_t id = cells->GetCellNumber(pos);
a186e444 930 if (id>24*48*fGeom->GetNumberOfSuperModules()) {
931 AliError(Form("Id of cell found to be too large: %d", id));
932 continue;
933 }
934 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigits->At(pos));
935 if (digit->GetId() != id) {
936 AliError(Form("Ids should be equal: %d %d", id, digit->GetId()));
937 return;
938 }
939 map[id] = pos;
3a952328 940 }
941
8c56d760 942 TObjArray filtMcParts;
cfd7d5b2 943 if (fMcParts) {
8c56d760 944 Int_t nmc = fMcParts->GetEntries();
945 for (Int_t i=0; i<nmc; ++i) {
946 AliStaPart *pa = static_cast<AliStaPart*>(fMcParts->At(i));
947 if (pa->OnEmcal())
948 filtMcParts.Add(pa);
949 }
950 }
951
fe38566a 952 Double_t vertex[3] = {0,0,0};
953 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
954
3a952328 955 for(Int_t i=0, ncl=0; i<nclus; ++i) {
956 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
b5ba4932 957
3a952328 958 if (!clus)
959 continue;
960 if (!clus->IsEMCAL())
961 continue;
962 if (clus->E()<fMinE)
963 continue;
964
fe38566a 965 Float_t clsPos[3] = {0,0,0};
3a952328 966 clus->GetPosition(clsPos);
967 TVector3 clsVec(clsPos);
3a952328 968 TLorentzVector clusterVec;
969 clus->GetMomentum(clusterVec,vertex);
970 Double_t clsEta = clusterVec.Eta();
971
972 AliStaCluster *cl = static_cast<AliStaCluster*>(fClusters->New(ncl++));
973 cl->fE = clus->E();
974 cl->fR = clsVec.Perp();
975 cl->fEta = clsVec.Eta();
976 cl->fPhi = clsVec.Phi();
977 cl->fN = clus->GetNCells();
978 cl->fN1 = GetNCells(clus,0.100);
979 cl->fN3 = GetNCells(clus,0.300);
980 Short_t id = -1;
981 Double_t emax = GetMaxCellEnergy(clus, id);
982 cl->fIdMax = id;
2ee7bda4 983 cl->fSM = fGeom->GetSuperModuleNumber(id);
3a952328 984 cl->fEmax = emax;
1f41ed3d 985 Short_t id2 = -1;
986 cl->fE2max = GetSecondMaxCellEnergy(clus,id2);
2ee7bda4 987 cl->fTmax = cells->GetCellTime(id);
3a952328 988 if (clus->GetDistanceToBadChannel()<10000)
989 cl->fDbc = clus->GetDistanceToBadChannel();
990 if (!TMath::IsNaN(clus->GetDispersion()))
991 cl->fDisp = clus->GetDispersion();
992 if (!TMath::IsNaN(clus->GetM20()))
0fbe8d4f 993 cl->fM20 = clus->GetM20();
3a952328 994 if (!TMath::IsNaN(clus->GetM02()))
0fbe8d4f 995 cl->fM02 = clus->GetM02();
2ee7bda4 996 Double_t maxAxis = -1, minAxis = -1;
3a952328 997 GetSigma(clus,maxAxis,minAxis);
2ee7bda4 998 clus->SetTOF(maxAxis); // store sigma in TOF for later plotting
3a952328 999 cl->fSig = maxAxis;
2ee7bda4 1000 Double_t sEtaEta = -1;
1001 Double_t sPhiPhi = -1;
5fc7508c 1002 GetSigmaEtaEta(clus, sEtaEta, sPhiPhi);
1003 cl->fSigEtaEta = sEtaEta;
1004 cl->fSigPhiPhi = sPhiPhi;
2ee7bda4 1005 Double_t clusterEcc = -1;
3a952328 1006 if (maxAxis > 0)
1007 clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
2ee7bda4 1008 clus->SetChi2(clusterEcc); // store ecc in chi2 for later plotting
1009 cl->fEcc = clusterEcc;
1010 cl->fTrIso = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
1011 cl->fTrIso1 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 1);
1012 cl->fTrIso2 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 2);
5fc7508c 1013 cl->fTrIsoD1 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist-0.1);
1014 cl->fTrIso1D1 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist-0.1, 1);
1015 cl->fTrIso2D1 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist-0.1, 2);
1016 cl->fTrIsoD3 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.1);
1017 cl->fTrIso1D3 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.1, 1);
1018 cl->fTrIso2D3 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.1, 2);
1f41ed3d 1019 cl->fTrIsoD4 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.2);
1020 cl->fTrIso1D4 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.2, 1);
1021 cl->fTrIso2D4 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.2, 2);
5fc7508c 1022 cl->fTrIsoStrip = GetTrackIsoStrip(clusterVec.Eta(), clusterVec.Phi());
2ee7bda4 1023 cl->fCeCore = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.05);
1024 cl->fCeIso = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
1025 cl->fCeIso1 = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.10);
1026 cl->fCeIso3 = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.30);
1f41ed3d 1027 cl->fCeIso4 = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.40);
a186e444 1028 cl->fCeIso3x3 = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 3, 3);
2ee7bda4 1029 cl->fCeIso4x4 = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 4, 4);
1030 cl->fCeIso5x5 = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 5, 5);
1031 cl->fCeIso3x22 = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 3, 22);
1032 cl->fIsShared = IsShared(clus);
1033 cl->fTrigId = -1;
1034 cl->fTrigE = 0;
e0e1022c 1035
2ee7bda4 1036 if (fTriggers) {
1037 Int_t ntrig = fTriggers->GetEntries();
1038 for (Int_t j = 0; j<ntrig; ++j) {
1039 AliStaTrigger *sta = static_cast<AliStaTrigger*>(fTriggers->At(j));
1040 if (!sta)
38727e64 1041 continue;
2ee7bda4 1042 Short_t idmax = sta->fIdMax;
1043 Bool_t inc = IsIdPartOfCluster(clus, idmax);
1044 if (inc) {
1045 cl->fTrigId = j;
1046 cl->fTrigE = sta->fE;
1047 break;
1048 }
38727e64 1049 }
3a952328 1050 }
1051
1052 // track matching
a186e444 1053 if (fDoTrMatGeom) {
1054 Double_t mind2 = 1e10;
1055 for(Int_t j = 0; j<ntrks; ++j) {
1056 AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
1057 if (!track)
1058 continue;
1059 if (TMath::Abs(clsEta-track->Eta())>0.5)
1060 continue;
1061 TVector3 vec(clsPos);
1062 Float_t tmpEta=-1, tmpPhi=-1, tmpR2=1e10;
1063 Double_t dedx = -1;
1064 if (fEsdEv) {
1065 AliESDtrack *esdTrack = static_cast<AliESDtrack*>(track);
1066 dedx = esdTrack->GetTPCsignal();
1067 }
3a952328 1068 if (TMath::Abs(clsEta-track->Eta())>fIsoDist)
1069 continue;
e0e1022c 1070 AliExternalTrackParam tParam;
1071 tParam.CopyFromVTrack(track);
1072 if (!fReco->ExtrapolateTrackToCluster(&tParam, clus, tmpEta, tmpPhi))
3a952328 1073 continue;
e0e1022c 1074 tmpR2 = tmpEta*tmpEta + tmpPhi*tmpPhi;
a186e444 1075 Double_t d2 = tmpR2;
1076 if (mind2>d2) {
1077 mind2=d2;
1078 cl->fTrDz = tmpEta;
1079 cl->fTrDr = tmpPhi;
1080 cl->fTrEp = clus->E()/track->P();
1081 cl->fTrDedx = dedx;
1082 cl->fIsTrackM = 1;
1083 }
b6c599fe 1084 }
a186e444 1085 if (cl->fIsTrackM) {
1086 fHMatchDr->Fill(cl->fTrDr);
1087 fHMatchDz->Fill(cl->fTrDz);
1088 fHMatchEp->Fill(cl->fTrEp);
3a952328 1089 }
1090 }
8c56d760 1091
1092 //mc matching
cfd7d5b2 1093 if (fMcParts) {
8c56d760 1094 Int_t nmc = filtMcParts.GetEntries();
1095 Double_t diffR2 = 1e9;
1096 AliStaPart *msta = 0;
1097 for (Int_t j=0; j<nmc; ++j) {
1098 AliStaPart *pa = static_cast<AliStaPart*>(filtMcParts.At(j));
1099 Double_t t1=clsVec.Eta()-pa->fVEta;
1100 Double_t t2=TVector2::Phi_mpi_pi(clsVec.Phi()-pa->fVPhi);
1101 Double_t tmp = t1*t1+t2*t2;
1102 if (tmp<diffR2) {
1103 diffR2 = tmp;
1104 msta = pa;
1105 }
1106 }
c2fe5f0e 1107 if (diffR2<10 && msta!=0) {
1108 cl->fMcLabel = msta->fLab;
1109 }
1110 }
1111
1112 cl->fEmbE = 0;
1113 if (fDigits && fEmbedMode) {
1114 for(Int_t j=0; j<cl->fN; ++j) {
1115 Short_t cid = TMath::Abs(clus->GetCellAbsId(j));
1116 Short_t pos = -1;
1117 std::map<Short_t,Short_t>::iterator it = map.find(cid);
1118 if (it!=map.end())
1119 pos = it->second;
1120 if (pos<0)
1121 continue;
1122 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigits->At(pos));
1123 if (!digit)
1124 continue;
c2fe5f0e 1125 if (digit->GetType()<-1) {
1126 cl->fEmbE += digit->GetChi2();
1127 }
1128 }
8c56d760 1129 }
b6c599fe 1130 }
1131}
1132
1133//________________________________________________________________________
1134void AliAnalysisTaskEMCALPi0PbPb::CalcPrimTracks()
1135{
a2de5ca1 1136 // Calculate track properties for primary tracks.
b6c599fe 1137
1138 fSelPrimTracks->Clear();
1139
1140 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1141 TClonesArray *tracks = 0;
e0e1022c 1142 Bool_t doConstrain = kFALSE;
1143 TString trkname(fPrimTracksName);
b6c599fe 1144 if (fEsdEv) {
e0e1022c 1145 if (trkname.IsNull()) {
1146 trkname = "Tracks";
1147 am->LoadBranch("Tracks");
1148 fSelPrimTracks->SetOwner(kTRUE);
1149 doConstrain = kTRUE;
1150 }
b6c599fe 1151 TList *l = fEsdEv->GetList();
e0e1022c 1152 tracks = dynamic_cast<TClonesArray*>(l->FindObject(trkname));
b6c599fe 1153 } else {
e0e1022c 1154 trkname = "tracks";
b6c599fe 1155 am->LoadBranch("tracks");
1156 TList *l = fAodEv->GetList();
1157 tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1158 }
1159
e0e1022c 1160 if (!tracks) {
1161 AliError(Form("Pointer to tracks %s == 0", trkname.Data() ));
296ea9b4 1162 return;
e0e1022c 1163 }
296ea9b4 1164
0ec74551 1165 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
b6c599fe 1166 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-fIsoDist*1.25;
1167 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+fIsoDist*1.25;
0ec74551 1168
296ea9b4 1169 if (fEsdEv) {
0ec74551 1170 am->LoadBranch("PrimaryVertex.");
1171 const AliESDVertex *vtx = fEsdEv->GetPrimaryVertexSPD();
1172 am->LoadBranch("SPDVertex.");
1173 const AliESDVertex *vtxSPD = fEsdEv->GetPrimaryVertexSPD();
1174 am->LoadBranch("Tracks");
1175 const Int_t Ntracks = tracks->GetEntries();
1176 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1177 AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1178 if (!track) {
1179 AliWarning(Form("Could not receive track %d\n", iTracks));
1180 continue;
1181 }
e0e1022c 1182 if (fPrimTrCuts && !fPrimTrCuts->IsSelected(track))
296ea9b4 1183 continue;
1184 Double_t eta = track->Eta();
1185 if (eta<-1||eta>1)
1186 continue;
0ec74551 1187 if (track->Phi()<phimin||track->Phi()>phimax)
1188 continue;
e0e1022c 1189 if (!doConstrain) {
1190 fSelPrimTracks->Add(track);
1191 continue;
1192 }
0ec74551 1193 AliESDtrack copyt(*track);
1194 Double_t bfield[3];
1195 copyt.GetBxByBz(bfield);
1196 AliExternalTrackParam tParam;
1197 Bool_t relate = copyt.RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&tParam);
1198 if (!relate)
1199 continue;
1200 copyt.Set(tParam.GetX(),tParam.GetAlpha(),tParam.GetParameter(),tParam.GetCovariance());
1201 Double_t p[3] = { 0. };
1202 copyt.GetPxPyPz(p);
1203 Double_t pos[3] = { 0. };
1204 copyt.GetXYZ(pos);
1205 Double_t covTr[21] = { 0. };
1206 copyt.GetCovarianceXYZPxPyPz(covTr);
1207 Double_t pid[10] = { 0. };
1208 copyt.GetESDpid(pid);
1209 AliAODTrack *aTrack = new AliAODTrack(copyt.GetID(),
1210 copyt.GetLabel(),
1211 p,
1212 kTRUE,
1213 pos,
1214 kFALSE,
1215 covTr,
1216 (Short_t)copyt.GetSign(),
1217 copyt.GetITSClusterMap(),
1218 pid,
1219 0,/*fPrimaryVertex,*/
1220 kTRUE, // check if this is right
1221 vtx->UsesTrack(copyt.GetID()));
1222 aTrack->SetTPCClusterMap(copyt.GetTPCClusterMap());
1223 aTrack->SetTPCSharedMap (copyt.GetTPCSharedMap());
1224 Float_t ndf = copyt.GetTPCNcls() + 1 - 5 ;
1225 if(ndf>0)
1226 aTrack->SetChi2perNDF(copyt.GetTPCchi2()/ndf);
1227 else
1228 aTrack->SetChi2perNDF(-1);
1229 aTrack->SetFlags(copyt.GetStatus());
1230 aTrack->SetTPCPointsF(copyt.GetTPCNclsF());
b6c599fe 1231 fSelPrimTracks->Add(aTrack);
296ea9b4 1232 }
1233 } else {
1234 Int_t ntracks = tracks->GetEntries();
1235 for (Int_t i=0; i<ntracks; ++i) {
1236 AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
1237 if (!track)
1238 continue;
1239 Double_t eta = track->Eta();
1240 if (eta<-1||eta>1)
1241 continue;
0ec74551 1242 if (track->Phi()<phimin||track->Phi()>phimax)
1243 continue;
b6c599fe 1244 if(track->GetTPCNcls()<fMinNClusPerTr)
296ea9b4 1245 continue;
b6c599fe 1246 //todo: Learn how to set/filter AODs for prim/sec tracks
1247 fSelPrimTracks->Add(track);
296ea9b4 1248 }
1249 }
1250}
1251
1252//________________________________________________________________________
3a952328 1253void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
296ea9b4 1254{
3a952328 1255 // Calculate track properties (including secondaries).
b6c599fe 1256
3a952328 1257 fSelTracks->Clear();
296ea9b4 1258
3a952328 1259 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1260 TClonesArray *tracks = 0;
1261 if (fEsdEv) {
1262 am->LoadBranch("Tracks");
1263 TList *l = fEsdEv->GetList();
1264 tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
1265 } else {
1266 am->LoadBranch("tracks");
1267 TList *l = fAodEv->GetList();
1268 tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1269 }
296ea9b4 1270
3a952328 1271 if (!tracks)
1272 return;
296ea9b4 1273
3a952328 1274 if (fEsdEv) {
1275 const Int_t Ntracks = tracks->GetEntries();
1276 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1277 AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1278 if (!track) {
1279 AliWarning(Form("Could not receive track %d\n", iTracks));
1280 continue;
1281 }
1282 if (fTrCuts && !fTrCuts->IsSelected(track))
1283 continue;
1284 Double_t eta = track->Eta();
1285 if (eta<-1||eta>1)
1286 continue;
1287 fSelTracks->Add(track);
1288 }
1289 } else {
1290 Int_t ntracks = tracks->GetEntries();
1291 for (Int_t i=0; i<ntracks; ++i) {
1292 AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
296ea9b4 1293 if (!track)
1294 continue;
3a952328 1295 Double_t eta = track->Eta();
1296 if (eta<-1||eta>1)
c4236a58 1297 continue;
3a952328 1298 if(track->GetTPCNcls()<fMinNClusPerTr)
b6c599fe 1299 continue;
2e4d8148 1300
3a952328 1301 if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL
e0e1022c 1302 AliExternalTrackParam tParam;
1303 tParam.CopyFromVTrack(track);
a186e444 1304 if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 430, 0.139, 1, kTRUE)) {
3a952328 1305 Double_t trkPos[3];
1306 tParam.GetXYZ(trkPos);
1307 track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]);
1308 }
b6c599fe 1309 }
3a952328 1310 fSelTracks->Add(track);
c4236a58 1311 }
296ea9b4 1312 }
1313}
1314
323834f0 1315//________________________________________________________________________
3a952328 1316void AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner()
323834f0 1317{
296ea9b4 1318 // Run custer reconstruction afterburner.
323834f0 1319
1320 AliVCaloCells *cells = fEsdCells;
1321 if (!cells)
1322 cells = fAodCells;
1323
1324 if (!cells)
1325 return;
1326
1327 Int_t ncells = cells->GetNumberOfCells();
1328 if (ncells<=0)
1329 return;
1330
1331 Double_t cellMeanE = 0, cellSigE = 0;
1332 for (Int_t i = 0; i<ncells; ++i) {
1333 Double_t cellE = cells->GetAmplitude(i);
1334 cellMeanE += cellE;
1335 cellSigE += cellE*cellE;
1336 }
1337 cellMeanE /= ncells;
1338 cellSigE /= ncells;
1339 cellSigE -= cellMeanE*cellMeanE;
1340 if (cellSigE<0)
1341 cellSigE = 0;
1342 cellSigE = TMath::Sqrt(cellSigE / ncells);
1343
1344 Double_t subE = cellMeanE - 7*cellSigE;
1345 if (subE<0)
1346 return;
1347
1348 for (Int_t i = 0; i<ncells; ++i) {
1349 Short_t id=-1;
1350 Double_t amp=0,time=0;
1351 if (!cells->GetCell(i, id, amp, time))
1352 continue;
1353 amp -= cellMeanE;
1354 if (amp<0.001)
1355 amp = 0;
1356 cells->SetCell(i, id, amp, time);
1357 }
1358
1359 TObjArray *clusters = fEsdClusters;
1360 if (!clusters)
1361 clusters = fAodClusters;
323834f0 1362 if (!clusters)
1363 return;
1364
1365 Int_t nclus = clusters->GetEntries();
1366 for (Int_t i = 0; i<nclus; ++i) {
1367 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1368 if (!clus->IsEMCAL())
1369 continue;
1370 Int_t nc = clus->GetNCells();
1371 Double_t clusE = 0;
1372 UShort_t ids[100] = {0};
1373 Double_t fra[100] = {0};
1374 for (Int_t j = 0; j<nc; ++j) {
1375 Short_t id = TMath::Abs(clus->GetCellAbsId(j));
1376 Double_t cen = cells->GetCellAmplitude(id);
1377 clusE += cen;
1378 if (cen>0) {
1379 ids[nc] = id;
1380 ++nc;
1381 }
1382 }
1383 if (clusE<=0) {
1384 clusters->RemoveAt(i);
1385 continue;
1386 }
1387
1388 for (Int_t j = 0; j<nc; ++j) {
1389 Short_t id = ids[j];
1390 Double_t cen = cells->GetCellAmplitude(id);
1391 fra[j] = cen/clusE;
1392 }
1393 clus->SetE(clusE);
1394 AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus);
1395 if (aodclus) {
1396 aodclus->Clear("");
1397 aodclus->SetNCells(nc);
1398 aodclus->SetCellsAmplitudeFraction(fra);
1399 aodclus->SetCellsAbsId(ids);
1400 }
1401 }
1402 clusters->Compress();
1403}
1404
286b47a5 1405//________________________________________________________________________
1406void AliAnalysisTaskEMCALPi0PbPb::FillCellHists()
1407{
1408 // Fill histograms related to cell properties.
d595acbb 1409
90d5b88b 1410 AliVCaloCells *cells = fEsdCells;
1411 if (!cells)
1412 cells = fAodCells;
1413
296ea9b4 1414 if (!cells)
1415 return;
90d5b88b 1416
296ea9b4 1417 Int_t cellModCount[12] = {0};
1418 Double_t cellMaxE = 0;
1419 Double_t cellMeanE = 0;
1420 Int_t ncells = cells->GetNumberOfCells();
1421 if (ncells==0)
1422 return;
1423
1424 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
1425
1426 for (Int_t i = 0; i<ncells; ++i) {
1427 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
1428 Double_t cellE = cells->GetAmplitude(i);
1429 fHCellE->Fill(cellE);
1430 if (cellE>cellMaxE)
1431 cellMaxE = cellE;
1432 cellMeanE += cellE;
1433
1434 Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1;
1435 Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta);
1436 if (!ret) {
1437 AliError(Form("Could not get cell index for %d", absID));
1438 continue;
1439 }
1440 ++cellModCount[iSM];
1441 Int_t iPhi=-1, iEta=-1;
1442 fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta);
1443 fHColuRow[iSM]->Fill(iEta,iPhi,1);
1444 fHColuRowE[iSM]->Fill(iEta,iPhi,cellE);
1445 fHCellFreqNoCut[iSM]->Fill(absID);
2e4d8148 1446 if (cellE > 0.1) fHCellFreqCut100M[iSM]->Fill(absID);
1447 if (cellE > 0.3) fHCellFreqCut300M[iSM]->Fill(absID);
296ea9b4 1448 if (fHCellCheckE && fHCellCheckE[absID])
1449 fHCellCheckE[absID]->Fill(cellE);
2e4d8148 1450 fHCellFreqE[iSM]->Fill(absID, cellE);
296ea9b4 1451 }
1452 fHCellH->Fill(cellMaxE);
1453 cellMeanE /= ncells;
1454 fHCellM->Fill(cellMeanE);
1455 fHCellM2->Fill(cellMeanE*ncells/24/48/nsm);
1456 for (Int_t i=0; i<nsm; ++i)
1457 fHCellMult[i]->Fill(cellModCount[i]);
286b47a5 1458}
1459
1460//________________________________________________________________________
1461void AliAnalysisTaskEMCALPi0PbPb::FillClusHists()
1462{
90d5b88b 1463 // Fill histograms related to cluster properties.
fa443410 1464
90d5b88b 1465 TObjArray *clusters = fEsdClusters;
1466 if (!clusters)
1467 clusters = fAodClusters;
296ea9b4 1468 if (!clusters)
1469 return;
90d5b88b 1470
296ea9b4 1471 Double_t vertex[3] = {0};
1472 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1473
1474 Int_t nclus = clusters->GetEntries();
1475 for(Int_t i = 0; i<nclus; ++i) {
1476 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1477 if (!clus)
1478 continue;
1479 if (!clus->IsEMCAL())
1480 continue;
90d5b88b 1481 TLorentzVector clusterVec;
296ea9b4 1482 clus->GetMomentum(clusterVec,vertex);
3a952328 1483 Double_t maxAxis = clus->GetTOF(); //sigma
1484 Double_t clusterEcc = clus->Chi2(); //eccentricity
296ea9b4 1485 fHClustEccentricity->Fill(clusterEcc);
1486 fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi());
1487 fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt());
1488 fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E());
1489 fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis);
1490 fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E());
f5e0f1e2 1491 fHClustEnergyNCell->Fill(clus->E(),clus->GetNCells());
788ca675 1492 }
1493}
b3ee6797 1494
38727e64 1495//________________________________________________________________________
1496void AliAnalysisTaskEMCALPi0PbPb::CalcMcInfo()
1497{
1498 // Get Mc truth particle information.
1499
cfd7d5b2 1500 if (!fMcParts)
38727e64 1501 return;
1502
807016ea 1503 fMcParts->Clear();
1504
1505 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
1506 Double_t etamin = -0.7;
1507 Double_t etamax = +0.7;
1508 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
1509 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
1510
56fd6cb2 1511 if (fAodEv) {
807016ea 1512 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1513 am->LoadBranch(AliAODMCParticle::StdBranchName());
56fd6cb2 1514 TClonesArray *tca = dynamic_cast<TClonesArray*>(fAodEv->FindListObject(AliAODMCParticle::StdBranchName()));
1515 if (!tca)
1516 return;
1517
1518 Int_t nents = tca->GetEntries();
807016ea 1519 for(int it=0; it<nents; ++it) {
56fd6cb2 1520 AliAODMCParticle *part = static_cast<AliAODMCParticle*>(tca->At(it));
807016ea 1521 part->Print();
56fd6cb2 1522
b5ba4932 1523 // pion or eta meson or direct photon
56fd6cb2 1524 if(part->GetPdgCode() == 111) {
1525 } else if(part->GetPdgCode() == 221) {
b5ba4932 1526 } else if(part->GetPdgCode() == 22 ) {
1527 } else
56fd6cb2 1528 continue;
1529
1530 // primary particle
1531 Double_t dR = TMath::Sqrt((part->Xv()*part->Xv())+(part->Yv()*part->Yv()));
1532 if(dR > 1.0)
1533 continue;
56fd6cb2 1534
807016ea 1535 // kinematic cuts
1536 Double_t pt = part->Pt() ;
1537 if (pt<0.5)
1538 continue;
1539 Double_t eta = part->Eta();
1540 if (eta<etamin||eta>etamax)
1541 continue;
1542 Double_t phi = part->Phi();
1543 if (phi<phimin||phi>phimax)
1544 continue;
1545
1546 ProcessDaughters(part, it, tca);
56fd6cb2 1547 }
807016ea 1548 return;
1549 }
38727e64 1550
1551 AliMCEvent *mcEvent = MCEvent();
1552 if (!mcEvent)
1553 return;
1554
1555 const AliVVertex *evtVtx = mcEvent->GetPrimaryVertex();
1556 if (!evtVtx)
1557 return;
1558
1559 mcEvent->PreReadAll();
38727e64 1560
1561 Int_t nTracks = mcEvent->GetNumberOfPrimaries();
1562 for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
1563 AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
1564 if (!mcP)
1565 continue;
1566
b5ba4932 1567 // pion or eta meson or direct photon
807016ea 1568 if(mcP->PdgCode() == 111) {
1569 } else if(mcP->PdgCode() == 221) {
b5ba4932 1570 } else if(mcP->PdgCode() == 22 ) {
38727e64 1571 } else
1572 continue;
1573
1574 // primary particle
807016ea 1575 Double_t dR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
1576 (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
38727e64 1577 if(dR > 1.0)
1578 continue;
1579
38727e64 1580 // kinematic cuts
807016ea 1581 Double_t pt = mcP->Pt() ;
1582 if (pt<0.5)
38727e64 1583 continue;
807016ea 1584 Double_t eta = mcP->Eta();
1585 if (eta<etamin||eta>etamax)
38727e64 1586 continue;
807016ea 1587 Double_t phi = mcP->Phi();
1588 if (phi<phimin||phi>phimax)
1589 continue;
1590
1591 ProcessDaughters(mcP, iTrack, mcEvent);
38727e64 1592 }
1593}
1594
788ca675 1595//________________________________________________________________________
1596void AliAnalysisTaskEMCALPi0PbPb::FillNtuple()
1597{
1598 // Fill ntuple.
1599
1600 if (!fNtuple)
1601 return;
1602
1603 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1604 if (fAodEv) {
1605 fHeader->fRun = fAodEv->GetRunNumber();
1606 fHeader->fOrbit = fAodEv->GetHeader()->GetOrbitNumber();
1607 fHeader->fPeriod = fAodEv->GetHeader()->GetPeriodNumber();
1608 fHeader->fBx = fAodEv->GetHeader()->GetBunchCrossNumber();
1609 fHeader->fL0 = fAodEv->GetHeader()->GetL0TriggerInputs();
1610 fHeader->fL1 = fAodEv->GetHeader()->GetL1TriggerInputs();
1611 fHeader->fL2 = fAodEv->GetHeader()->GetL2TriggerInputs();
1612 fHeader->fTrClassMask = fAodEv->GetHeader()->GetTriggerMask();
1613 fHeader->fTrCluster = fAodEv->GetHeader()->GetTriggerCluster();
1614 fHeader->fOffTriggers = fAodEv->GetHeader()->GetOfflineTrigger();
1615 fHeader->fFiredTriggers = fAodEv->GetHeader()->GetFiredTriggerClasses();
1616 } else {
1617 fHeader->fRun = fEsdEv->GetRunNumber();
1618 fHeader->fOrbit = fEsdEv->GetHeader()->GetOrbitNumber();
1619 fHeader->fPeriod = fEsdEv->GetESDRun()->GetPeriodNumber();
1620 fHeader->fBx = fEsdEv->GetHeader()->GetBunchCrossNumber();
1621 fHeader->fL0 = fEsdEv->GetHeader()->GetL0TriggerInputs();
1622 fHeader->fL1 = fEsdEv->GetHeader()->GetL1TriggerInputs();
1623 fHeader->fL2 = fEsdEv->GetHeader()->GetL2TriggerInputs();
1624 fHeader->fTrClassMask = fEsdEv->GetHeader()->GetTriggerMask();
1625 fHeader->fTrCluster = fEsdEv->GetHeader()->GetTriggerCluster();
1626 fHeader->fOffTriggers = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
1627 fHeader->fFiredTriggers = fEsdEv->GetFiredTriggerClasses();
5fe1ca23 1628 Float_t v0CorrR = 0;
1629 fHeader->fV0 = AliESDUtils::GetCorrV0(fEsdEv,v0CorrR);
1630 const AliMultiplicity *mult = fEsdEv->GetMultiplicity();
1631 if (mult)
1632 fHeader->fCl1 = mult->GetNumberOfITSClusters(1);
1633 fHeader->fTr = AliESDtrackCuts::GetReferenceMultiplicity(fEsdEv,1);
c0d2e1f2 1634 AliTriggerAnalysis trAn; /// Trigger Analysis
1635 Bool_t v0B = trAn.IsOfflineTriggerFired(fEsdEv, AliTriggerAnalysis::kV0C);
1636 Bool_t v0A = trAn.IsOfflineTriggerFired(fEsdEv, AliTriggerAnalysis::kV0A);
469b2bff 1637 fHeader->fV0And = v0A && v0B;
1638 fHeader->fIsHT = (fHeader->fOffTriggers & AliVEvent::kEMC1) || (fHeader->fOffTriggers & AliVEvent::kEMC7);
1639 am->LoadBranch("SPDPileupVertices");
1640 am->LoadBranch("TrkPileupVertices");
1641 fHeader->fIsPileup = fEsdEv->IsPileupFromSPD(3,0.8);
1642 fHeader->fIsPileup2 = fEsdEv->IsPileupFromSPD(3,0.4);
1643 fHeader->fIsPileup4 = fEsdEv->IsPileupFromSPD(3,0.2);
1644 fHeader->fIsPileup8 = fEsdEv->IsPileupFromSPD(3,0.1);
1645 fHeader->fNSpdVertices = fEsdEv->GetNumberOfPileupVerticesSPD();
1646 fHeader->fNTpcVertices = fEsdEv->GetNumberOfPileupVerticesTracks();
788ca675 1647 }
2ee7bda4 1648
788ca675 1649 AliCentrality *cent = InputEvent()->GetCentrality();
1650 fHeader->fV0Cent = cent->GetCentralityPercentileUnchecked("V0M");
1651 fHeader->fCl1Cent = cent->GetCentralityPercentileUnchecked("CL1");
1652 fHeader->fTrCent = cent->GetCentralityPercentileUnchecked("TRK");
1653 fHeader->fCqual = cent->GetQuality();
b6c599fe 1654
1655 AliEventplane *ep = InputEvent()->GetEventplane();
1656 if (ep) {
1657 if (ep->GetQVector())
1658 fHeader->fPsi = ep->GetQVector()->Phi()/2. ;
f2b8fca6 1659 else
1660 fHeader->fPsi = -1;
b6c599fe 1661 if (ep->GetQsub1()&&ep->GetQsub2())
1662 fHeader->fPsiRes = ep->GetQsub1()->Phi()/2.-ep->GetQsub2()->Phi()/2.;
1663 else
1664 fHeader->fPsiRes = 0;
1665 }
1666
788ca675 1667 Double_t val = 0;
1668 TString trgclasses(fHeader->fFiredTriggers);
1669 for (Int_t j = 0; j<fTrClassNamesArr->GetEntries(); ++j) {
1670 const char *name = fTrClassNamesArr->At(j)->GetName();
2ee7bda4 1671 TRegexp regexp(name);
1672 if (trgclasses.Contains(regexp))
788ca675 1673 val += TMath::Power(2,j);
1674 }
1675 fHeader->fTcls = (UInt_t)val;
1676
5fe1ca23 1677 fHeader->fNSelTr = fSelTracks->GetEntries();
1678 fHeader->fNSelPrimTr = fSelPrimTracks->GetEntries();
f5e0f1e2 1679 fHeader->fNSelPrimTr1 = 0;
1680 fHeader->fNSelPrimTr2 = 0;
1681 for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); iTracks++){
1682 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1683 if(track->Pt()>1)
1684 ++fHeader->fNSelPrimTr1;
1685 if(track->Pt()>2)
1686 ++fHeader->fNSelPrimTr2;
1687 }
5fe1ca23 1688
1689 fHeader->fNCells = 0;
c0d2e1f2 1690 fHeader->fNCells0 = 0;
1691 fHeader->fNCells01 = 0;
1692 fHeader->fNCells03 = 0;
5fe1ca23 1693 fHeader->fNCells1 = 0;
1694 fHeader->fNCells2 = 0;
1695 fHeader->fNCells5 = 0;
1696 fHeader->fMaxCellE = 0;
1697
1698 AliVCaloCells *cells = fEsdCells;
1699 if (!cells)
1700 cells = fAodCells;
1701
1702 if (cells) {
1703 Int_t ncells = cells->GetNumberOfCells();
1704 for(Int_t j=0; j<ncells; ++j) {
1705 Double_t cellen = cells->GetAmplitude(j);
c0d2e1f2 1706 if (cellen>0.045)
1707 ++fHeader->fNCells0;
1708 if (cellen>0.1)
1709 ++fHeader->fNCells01;
1710 if (cellen>0.3)
1711 ++fHeader->fNCells03;
5fe1ca23 1712 if (cellen>1)
1713 ++fHeader->fNCells1;
1714 if (cellen>2)
1715 ++fHeader->fNCells2;
1716 if (cellen>5)
1717 ++fHeader->fNCells5;
1718 if (cellen>fHeader->fMaxCellE)
1719 fHeader->fMaxCellE = cellen;
1720 }
1721 fHeader->fNCells = ncells;
1722 }
1723
1724 fHeader->fNClus = 0;
1725 fHeader->fNClus1 = 0;
1726 fHeader->fNClus2 = 0;
1727 fHeader->fNClus5 = 0;
1728 fHeader->fMaxClusE = 0;
1729
1730 TObjArray *clusters = fEsdClusters;
1731 if (!clusters)
1732 clusters = fAodClusters;
1733
1734 if (clusters) {
1735 Int_t nclus = clusters->GetEntries();
1736 for(Int_t j=0; j<nclus; ++j) {
1737 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(j));
6d4faab0 1738 if (!clus->IsEMCAL())
1739 continue;
5fe1ca23 1740 Double_t clusen = clus->E();
1741 if (clusen>1)
1742 ++fHeader->fNClus1;
1743 if (clusen>2)
1744 ++fHeader->fNClus2;
1745 if (clusen>5)
1746 ++fHeader->fNClus5;
1747 if (clusen>fHeader->fMaxClusE)
1748 fHeader->fMaxClusE = clusen;
1749 }
1750 fHeader->fNClus = nclus;
1751 }
1752
2ee7bda4 1753 fHeader->fMaxTrE = 0;
1754 if (fTriggers) {
1755 Int_t ntrig = fTriggers->GetEntries();
1756 for (Int_t j = 0; j<ntrig; ++j) {
1757 AliStaTrigger *sta = static_cast<AliStaTrigger*>(fTriggers->At(j));
1758 if (!sta)
1759 continue;
1760 if (sta->fE>fHeader->fMaxTrE)
1761 fHeader->fMaxTrE = sta->fE;
1762 }
1763 }
1764
1765 // count cells above 100 MeV on super modules
1766 fHeader->fNcSM0 = GetNCells(0, 0.100);
1767 fHeader->fNcSM1 = GetNCells(1, 0.100);
1768 fHeader->fNcSM2 = GetNCells(2, 0.100);
1769 fHeader->fNcSM3 = GetNCells(3, 0.100);
1770 fHeader->fNcSM4 = GetNCells(4, 0.100);
1771 fHeader->fNcSM5 = GetNCells(5, 0.100);
1772 fHeader->fNcSM6 = GetNCells(6, 0.100);
1773 fHeader->fNcSM7 = GetNCells(7, 0.100);
1774 fHeader->fNcSM8 = GetNCells(8, 0.100);
1775 fHeader->fNcSM9 = GetNCells(9, 0.100);
1776
788ca675 1777 if (fAodEv) {
1778 am->LoadBranch("vertices");
1779 AliAODVertex *pv = fAodEv->GetPrimaryVertex();
1780 FillVertex(fPrimVert, pv);
1781 AliAODVertex *sv = fAodEv->GetPrimaryVertexSPD();
1782 FillVertex(fSpdVert, sv);
1783 } else {
1784 am->LoadBranch("PrimaryVertex.");
1785 const AliESDVertex *pv = fEsdEv->GetPrimaryVertexTracks();
1786 FillVertex(fPrimVert, pv);
1787 am->LoadBranch("SPDVertex.");
1788 const AliESDVertex *sv = fEsdEv->GetPrimaryVertexSPD();
1789 FillVertex(fSpdVert, sv);
1790 am->LoadBranch("TPCVertex.");
1791 const AliESDVertex *tv = fEsdEv->GetPrimaryVertexTPC();
1792 FillVertex(fTpcVert, tv);
fa443410 1793 }
788ca675 1794
788ca675 1795 fNtuple->Fill();
286b47a5 1796}
1797
1798//________________________________________________________________________
1799void AliAnalysisTaskEMCALPi0PbPb::FillPionHists()
1800{
1801 // Fill histograms related to pions.
286b47a5 1802
296ea9b4 1803 Double_t vertex[3] = {0};
fa443410 1804 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1805
90d5b88b 1806 TObjArray *clusters = fEsdClusters;
1807 if (!clusters)
1808 clusters = fAodClusters;
fa443410 1809
90d5b88b 1810 if (clusters) {
1811 TLorentzVector clusterVec1;
1812 TLorentzVector clusterVec2;
1813 TLorentzVector pionVec;
1814
1815 Int_t nclus = clusters->GetEntries();
fa443410 1816 for (Int_t i = 0; i<nclus; ++i) {
90d5b88b 1817 AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i));
fa443410 1818 if (!clus1)
1819 continue;
1820 if (!clus1->IsEMCAL())
1821 continue;
3c661da5 1822 if (clus1->E()<fMinE)
6eb6260e 1823 continue;
a49742b5 1824 if (clus1->GetNCells()<fNminCells)
1825 continue;
f224d35b 1826 if (GetMaxCellEnergy(clus1)/clus1->E()<fMinErat)
1827 continue;
3c661da5 1828 if (clus1->Chi2()<fMinEcc) // eccentricity cut
f224d35b 1829 continue;
fa443410 1830 clus1->GetMomentum(clusterVec1,vertex);
1831 for (Int_t j = i+1; j<nclus; ++j) {
90d5b88b 1832 AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j));
fa443410 1833 if (!clus2)
1834 continue;
1835 if (!clus2->IsEMCAL())
1836 continue;
3c661da5 1837 if (clus2->E()<fMinE)
6eb6260e 1838 continue;
a49742b5 1839 if (clus2->GetNCells()<fNminCells)
1840 continue;
f224d35b 1841 if (GetMaxCellEnergy(clus2)/clus2->E()<fMinErat)
1842 continue;
3c661da5 1843 if (clus2->Chi2()<fMinEcc) // eccentricity cut
f224d35b 1844 continue;
fa443410 1845 clus2->GetMomentum(clusterVec2,vertex);
1846 pionVec = clusterVec1 + clusterVec2;
1847 Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E();
6eb6260e 1848 Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect());
d595acbb 1849 if (pionZgg < fAsymMax) {
1850 fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi());
1851 fHPionMggPt->Fill(pionVec.M(),pionVec.Pt());
1852 fHPionMggAsym->Fill(pionVec.M(),pionZgg);
6eb6260e 1853 fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle);
d595acbb 1854 Int_t bin = fPtRanges->FindBin(pionVec.Pt());
1855 fHPionInvMasses[bin]->Fill(pionVec.M());
1856 }
fa443410 1857 }
1858 }
90d5b88b 1859 }
fa443410 1860}
d595acbb 1861
38727e64 1862//________________________________________________________________________
1863void AliAnalysisTaskEMCALPi0PbPb::FillMcHists()
1864{
1865 // Fill additional MC information histograms.
1866
cfd7d5b2 1867 if (!fMcParts)
38727e64 1868 return;
1869
1870 // check if aod or esd mc mode and the fill histos
1871}
1872
6eb6260e 1873//________________________________________________________________________
323834f0 1874void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists()
6eb6260e 1875{
788ca675 1876 // Fill other histograms.
b5ba4932 1877
1878 for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); ++iTracks){
1879 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1880 if(!track)
1881 continue;
1882 fHPrimTrackPt->Fill(track->Pt());
1883 fHPrimTrackEta->Fill(track->Eta());
1884 fHPrimTrackPhi->Fill(track->Phi());
1885 }
788ca675 1886}
1887
f5e0f1e2 1888//________________________________________________________________________
1889void AliAnalysisTaskEMCALPi0PbPb::FillTrackHists()
1890{
1891 // Fill track histograms.
1892
1893 if (fSelPrimTracks) {
1894 for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); iTracks++) {
1895 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1896 if(!track)
1897 continue;
1898 fHPrimTrackPt->Fill(track->Pt());
1899 fHPrimTrackEta->Fill(track->Eta());
1900 fHPrimTrackPhi->Fill(track->Phi());
1901 }
1902 }
1903}
1904
788ca675 1905//__________________________________________________________________________________________________
1906void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliESDVertex *esdv)
1907{
1908 // Fill vertex from ESD vertex info.
1909
1910 v->fVx = esdv->GetX();
1911 v->fVy = esdv->GetY();
1912 v->fVz = esdv->GetZ();
1913 v->fVc = esdv->GetNContributors();
1914 v->fDisp = esdv->GetDispersion();
1915 v->fZres = esdv->GetZRes();
1916 v->fChi2 = esdv->GetChi2();
1917 v->fSt = esdv->GetStatus();
1918 v->fIs3D = esdv->IsFromVertexer3D();
1919 v->fIsZ = esdv->IsFromVertexerZ();
1920}
1921
1922//__________________________________________________________________________________________________
1923void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliAODVertex *aodv)
1924{
1925 // Fill vertex from AOD vertex info.
1926
1927 v->fVx = aodv->GetX();
1928 v->fVy = aodv->GetY();
1929 v->fVz = aodv->GetZ();
1930 v->fVc = aodv->GetNContributors();
1931 v->fChi2 = aodv->GetChi2();
6eb6260e 1932}
1933
d595acbb 1934//________________________________________________________________________
296ea9b4 1935Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
1936{
1937 // Compute isolation based on cell content.
1938
1939 AliVCaloCells *cells = fEsdCells;
1940 if (!cells)
1941 cells = fAodCells;
1942 if (!cells)
1943 return 0;
1944
1945 Double_t cellIsolation = 0;
1946 Double_t rad2 = radius*radius;
1947 Int_t ncells = cells->GetNumberOfCells();
1948 for (Int_t i = 0; i<ncells; ++i) {
1949 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
0fbe8d4f 1950 Float_t eta=-1, phi=-1;
296ea9b4 1951 fGeom->EtaPhiFromIndex(absID,eta,phi);
0fbe8d4f 1952 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
296ea9b4 1953 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1954 if(dist>rad2)
1955 continue;
0fbe8d4f 1956 Double_t cellE = cells->GetAmplitude(i);
2ee7bda4 1957 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
1958 Double_t cellEt = cellE*sin(theta);
1959 cellIsolation += cellEt;
296ea9b4 1960 }
1961 return cellIsolation;
1962}
2ee7bda4 1963
5fc7508c 1964//________________________________________________________________________
1965Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsoNxM(Double_t cEta, Double_t cPhi, Int_t N, Int_t M) const
1966{
1967 // Compute isolation based on cell content, in a NxM rectangle.
1968
1969 AliVCaloCells *cells = fEsdCells;
1970 if (!cells)
1971 cells = fAodCells;
1972 if (!cells)
1973 return 0;
1974
1975 Double_t cellIsolation = 0;
1976 Int_t ncells = cells->GetNumberOfCells();
1977 for (Int_t i = 0; i<ncells; ++i) {
1978 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
1979 Float_t eta=-1, phi=-1;
1980 fGeom->EtaPhiFromIndex(absID,eta,phi);
1981 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
1982 Double_t etadiff = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1983 if(TMath::Abs(etadiff)/0.014>N)
1984 continue;
1985 if(TMath::Abs(phidiff)/0.014>M)
1986 continue;
1987 Double_t cellE = cells->GetAmplitude(i);
2ee7bda4 1988 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
1989 Double_t cellEt = cellE*sin(theta);
1990 cellIsolation += cellEt;
5fc7508c 1991 }
1992 return cellIsolation;
1993}
1994
296ea9b4 1995//________________________________________________________________________
0fbe8d4f 1996Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellEnergy(const AliVCluster *cluster) const
1997{
1998 // Get maximum energy of attached cell.
1999
2000 Double_t ret = 0;
2001 Int_t ncells = cluster->GetNCells();
2002 if (fEsdCells) {
2003 for (Int_t i=0; i<ncells; i++) {
2004 Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2005 ret += e;
2006 }
2007 } else {
2008 for (Int_t i=0; i<ncells; i++) {
2009 Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2010 ret += e;
2011 }
2012 }
2013 return ret;
2014}
2015
2016//________________________________________________________________________
2017Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
d595acbb 2018{
90d5b88b 2019 // Get maximum energy of attached cell.
2020
788ca675 2021 id = -1;
90d5b88b 2022 Double_t maxe = 0;
90d5b88b 2023 Int_t ncells = cluster->GetNCells();
2024 if (fEsdCells) {
2025 for (Int_t i=0; i<ncells; i++) {
27422847 2026 Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
f0897c18 2027 if (e>maxe) {
90d5b88b 2028 maxe = e;
788ca675 2029 id = cluster->GetCellAbsId(i);
f0897c18 2030 }
90d5b88b 2031 }
2032 } else {
2033 for (Int_t i=0; i<ncells; i++) {
27422847 2034 Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
90d5b88b 2035 if (e>maxe)
2036 maxe = e;
788ca675 2037 id = cluster->GetCellAbsId(i);
90d5b88b 2038 }
2039 }
6eb6260e 2040 return maxe;
d595acbb 2041}
2042
469b2bff 2043//________________________________________________________________________
1f41ed3d 2044Double_t AliAnalysisTaskEMCALPi0PbPb::GetSecondMaxCellEnergy(AliVCluster *clus, Short_t &id) const
469b2bff 2045{
2046 // Get second maximum cell.
2047
2048 AliVCaloCells *cells = fEsdCells;
2049 if (!cells)
2050 cells = fAodCells;
2051 if (!cells)
2052 return -1;
2053
2054 Double_t secondEmax=0, firstEmax=0;
2055 Double_t cellen;
2056 for(Int_t iCell=0;iCell<clus->GetNCells();iCell++){
2057 Int_t absId = clus->GetCellAbsId(iCell);
2058 cellen = cells->GetCellAmplitude(absId);
2059 if(cellen > firstEmax)
2060 firstEmax = cellen;
2061 }
2062 for(Int_t iCell=0;iCell<clus->GetNCells();iCell++){
2063 Int_t absId = clus->GetCellAbsId(iCell);
2064 cellen = cells->GetCellAmplitude(absId);
1f41ed3d 2065 if(cellen < firstEmax && cellen > secondEmax) {
469b2bff 2066 secondEmax = cellen;
1f41ed3d 2067 id = absId;
2068 }
469b2bff 2069 }
2070 return secondEmax;
2071}
2072
d595acbb 2073//________________________________________________________________________
0fbe8d4f 2074void AliAnalysisTaskEMCALPi0PbPb::GetSigma(const AliVCluster *c, Double_t& sigmaMax, Double_t &sigmaMin) const
d595acbb 2075{
90d5b88b 2076 // Calculate the (E) weighted variance along the longer (eigen) axis.
2077
6bf90832 2078 sigmaMax = 0; // cluster variance along its longer axis
2079 sigmaMin = 0; // cluster variance along its shorter axis
2080 Double_t Ec = c->E(); // cluster energy
296ea9b4 2081 if(Ec<=0)
2082 return;
6bf90832 2083 Double_t Xc = 0 ; // cluster first moment along X
2084 Double_t Yc = 0 ; // cluster first moment along Y
2085 Double_t Sxx = 0 ; // cluster second central moment along X (variance_X^2)
2086 Double_t Sxy = 0 ; // cluster second central moment along Y (variance_Y^2)
2087 Double_t Syy = 0 ; // cluster covariance^2
90d5b88b 2088
2089 AliVCaloCells *cells = fEsdCells;
2090 if (!cells)
2091 cells = fAodCells;
2092
6eb6260e 2093 if (!cells)
2094 return;
2095
6bf90832 2096 Int_t ncells = c->GetNCells();
6eb6260e 2097 if (ncells==1)
2098 return;
2099
6eb6260e 2100 for(Int_t j=0; j<ncells; ++j) {
6bf90832 2101 Int_t id = TMath::Abs(c->GetCellAbsId(j));
6eb6260e 2102 Double_t cellen = cells->GetCellAmplitude(id);
3a952328 2103 TVector3 pos;
2104 fGeom->GetGlobal(id,pos);
6eb6260e 2105 Xc += cellen*pos.X();
2106 Yc += cellen*pos.Y();
2107 Sxx += cellen*pos.X()*pos.X();
2108 Syy += cellen*pos.Y()*pos.Y();
2109 Sxy += cellen*pos.X()*pos.Y();
2110 }
2111 Xc /= Ec;
2112 Yc /= Ec;
2113 Sxx /= Ec;
2114 Syy /= Ec;
2115 Sxy /= Ec;
2116 Sxx -= Xc*Xc;
2117 Syy -= Yc*Yc;
2118 Sxy -= Xc*Yc;
f0897c18 2119 Sxx = TMath::Abs(Sxx);
2120 Syy = TMath::Abs(Syy);
296ea9b4 2121 sigmaMax = (Sxx + Syy + TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
2122 sigmaMax = TMath::Sqrt(TMath::Abs(sigmaMax));
2123 sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
2124 sigmaMin = TMath::Sqrt(TMath::Abs(sigmaMin));
d595acbb 2125}
2ee7bda4 2126
5fc7508c 2127//________________________________________________________________________
2128void AliAnalysisTaskEMCALPi0PbPb::GetSigmaEtaEta(const AliVCluster *c, Double_t& sEtaEta, Double_t &sPhiPhi) const
2129{
2130 // Calculate the (E) weighted variance along the pseudorapidity.
2131
2132 sEtaEta = 0;
2133 sPhiPhi = 0;
2134
2135 Double_t Ec = c->E(); // cluster energy
2136 if(Ec<=0)
2137 return;
2138
2139 const Int_t ncells = c->GetNCells();
2140
2141 Double_t EtaC = 0; // cluster first moment along eta
2142 Double_t PhiC = 0; // cluster first moment along phi
2143 Double_t Setaeta = 0; // cluster second central moment along eta
2144 Double_t Sphiphi = 0; // cluster second central moment along phi
2145 Double_t w[ncells]; // weight max(0,4.5*log(E_i/Ec))
2146 Double_t sumw = 0;
2147 Int_t id[ncells];
2148
2149 AliVCaloCells *cells = fEsdCells;
2150 if (!cells)
2151 cells = fAodCells;
2152
2153 if (!cells)
2154 return;
2155
2156 if (ncells==1)
2157 return;
2158
2159 for(Int_t j=0; j<ncells; ++j) {
2160 id[j] = TMath::Abs(c->GetCellAbsId(j));
2161 Double_t cellen = cells->GetCellAmplitude(id[j]);
2162 w[j] = TMath::Max(0., 4.5+TMath::Log(cellen/Ec));
2163 TVector3 pos;
2164 fGeom->GetGlobal(id[j],pos);
2165 EtaC += w[j]*pos.Eta();
2166 PhiC += w[j]*pos.Phi();
2167 sumw += w[j];
2168 }
2169 EtaC /= sumw;
2170 PhiC /= sumw;
2171
2172 for(Int_t j=0; j<ncells; ++j) {
2173 TVector3 pos;
2174 fGeom->GetGlobal(id[j],pos);
2175 Setaeta = w[j]*(pos.Eta() - EtaC)*(pos.Eta() - EtaC);
2176 Sphiphi = w[j]*(pos.Phi() - PhiC)*(pos.Phi() - PhiC);
2177 }
2178 Setaeta /= sumw;
2179 sEtaEta = TMath::Sqrt(Setaeta);
2180 Sphiphi /= sumw;
2181 sPhiPhi = TMath::Sqrt(Sphiphi);
2182}
2183
6bf90832 2184//________________________________________________________________________
0fbe8d4f 2185Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(const AliVCluster *c, Double_t emin) const
6bf90832 2186{
2187 // Calculate number of attached cells above emin.
2188
6bf90832 2189 AliVCaloCells *cells = fEsdCells;
2190 if (!cells)
2191 cells = fAodCells;
6bf90832 2192 if (!cells)
296ea9b4 2193 return 0;
6bf90832 2194
296ea9b4 2195 Int_t n = 0;
6bf90832 2196 Int_t ncells = c->GetNCells();
2197 for(Int_t j=0; j<ncells; ++j) {
2198 Int_t id = TMath::Abs(c->GetCellAbsId(j));
2199 Double_t cellen = cells->GetCellAmplitude(id);
2200 if (cellen>=emin)
2201 ++n;
2202 }
2203 return n;
2204}
2205
2ee7bda4 2206//________________________________________________________________________
2207Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(Int_t sm, Double_t emin) const
2208{
2209 // Calculate number of cells per SM above emin.
2210
2211 AliVCaloCells *cells = fEsdCells;
2212 if (!cells)
2213 cells = fAodCells;
2214 if (!cells)
2215 return 0;
2216
2217 Int_t n = 0;
2218 Int_t ncells = cells->GetNumberOfCells();
2219 for(Int_t j=0; j<ncells; ++j) {
2220 Int_t id = TMath::Abs(cells->GetCellNumber(j));
2221 Double_t cellen = cells->GetCellAmplitude(id);
2222 if (cellen<emin)
2223 continue;
2224 Int_t fsm = fGeom->GetSuperModuleNumber(id);
2225 if (fsm != sm)
2226 continue;
2227 ++n;
2228 }
2229 return n;
2230}
2231
296ea9b4 2232//________________________________________________________________________
b6c599fe 2233Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
296ea9b4 2234{
2235 // Compute isolation based on tracks.
2236
2237 Double_t trkIsolation = 0;
2238 Double_t rad2 = radius*radius;
b6c599fe 2239 Int_t ntrks = fSelPrimTracks->GetEntries();
296ea9b4 2240 for(Int_t j = 0; j<ntrks; ++j) {
1f41ed3d 2241 AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
296ea9b4 2242 if (!track)
2243 continue;
b6c599fe 2244 if (track->Pt()<pt)
2245 continue;
296ea9b4 2246 Float_t eta = track->Eta();
2247 Float_t phi = track->Phi();
0fbe8d4f 2248 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
296ea9b4 2249 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
2250 if(dist>rad2)
2251 continue;
2252 trkIsolation += track->Pt();
2253 }
2254 return trkIsolation;
2255}
2ee7bda4 2256
5fc7508c 2257//________________________________________________________________________
2258Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsoStrip(Double_t cEta, Double_t cPhi, Double_t dEta, Double_t dPhi, Double_t pt) const
2259{
2260 // Compute isolation based on tracks.
2261
2262 Double_t trkIsolation = 0;
2263 Int_t ntrks = fSelPrimTracks->GetEntries();
2264 for(Int_t j = 0; j<ntrks; ++j) {
1f41ed3d 2265 AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
5fc7508c 2266 if (!track)
2267 continue;
2268 if (track->Pt()<pt)
2269 continue;
2270 Float_t eta = track->Eta();
2271 Float_t phi = track->Phi();
2272 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
2273 Double_t etadiff = (eta-cEta);
2274 if(TMath::Abs(etadiff)>dEta || TMath::Abs(phidiff)>dPhi)
2275 continue;
2276 trkIsolation += track->Pt();
2277 }
2278 return trkIsolation;
2279}
2280
0fbe8d4f 2281//________________________________________________________________________
2282Bool_t AliAnalysisTaskEMCALPi0PbPb::IsShared(const AliVCluster *c) const
2283{
2284 // Returns if cluster shared across super modules.
2285
1f41ed3d 2286 if (!c)
0fbe8d4f 2287 return 0;
2288
2289 Int_t n = -1;
2290 Int_t ncells = c->GetNCells();
2291 for(Int_t j=0; j<ncells; ++j) {
2292 Int_t id = TMath::Abs(c->GetCellAbsId(j));
2293 Int_t got = id / (24*48);
2294 if (n==-1) {
2295 n = got;
2296 continue;
2297 }
2298 if (got!=n)
2299 return 1;
2300 }
2301 return 0;
2302}
3a952328 2303
2ee7bda4 2304//________________________________________________________________________
2305Bool_t AliAnalysisTaskEMCALPi0PbPb::IsIdPartOfCluster(const AliVCluster *c, Short_t id) const
2306{
2307 // Returns if id is part of cluster.
2308
2309 AliVCaloCells *cells = fEsdCells;
2310 if (!cells)
2311 cells = fAodCells;
2312 if (!cells)
2313 return 0;
2314
2315 Int_t ncells = c->GetNCells();
2316 for(Int_t j=0; j<ncells; ++j) {
2317 Int_t cid = TMath::Abs(c->GetCellAbsId(j));
2318 if (cid == id)
2319 return 1;
2320 }
2321 return 0;
2322}
2323
38727e64 2324//________________________________________________________________________
807016ea 2325void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(const AliVParticle *p, const TObjArray *arr, Int_t level) const
38727e64 2326{
56fd6cb2 2327 // Print recursively daughter information.
2328
2329 if (!p || !arr)
2330 return;
2331
807016ea 2332 const AliAODMCParticle *amc = dynamic_cast<const AliAODMCParticle*>(p);
2333 if (!amc)
2334 return;
2335 for (Int_t i=0; i<level; ++i) printf(" ");
2336 amc->Print();
2337
2338 Int_t n = amc->GetNDaughters();
2339 for (Int_t i=0; i<n; ++i) {
2340 Int_t d = amc->GetDaughter(i);
2341 const AliVParticle *dmc = static_cast<const AliVParticle*>(arr->At(d));
2342 PrintDaughters(dmc,arr,level+1);
2343 }
2344}
2345
2346//________________________________________________________________________
2347void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(const AliMCParticle *p, const AliMCEvent *arr, Int_t level) const
2348{
2349 // Print recursively daughter information.
2350
2351 if (!p || !arr)
2352 return;
2353
2354 for (Int_t i=0; i<level; ++i) printf(" ");
2355 Int_t d1 = p->GetFirstDaughter();
2356 Int_t d2 = p->GetLastDaughter();
2357 printf("pid=%d: %.2f %.2f %.2f (%.2f %.2f %.2f); nd=%d,%d\n",
2358 p->PdgCode(),p->Px(),p->Py(),p->Pz(),p->Xv(),p->Yv(),p->Zv(),d1,d2);
2359 if (d1<0)
2360 return;
2361 if (d2<0)
2362 d2=d1;
2363 for (Int_t i=d1;i<=d2;++i) {
2364 const AliMCParticle *dmc = static_cast<const AliMCParticle *>(arr->GetTrack(i));
6d4faab0 2365 PrintDaughters(dmc,arr,level+1);
56fd6cb2 2366 }
38727e64 2367}
2368
2369//________________________________________________________________________
2370void AliAnalysisTaskEMCALPi0PbPb::PrintTrackRefs(AliMCParticle *p) const
2371{
56fd6cb2 2372 // Print track ref array.
2373
2374 if (!p)
2375 return;
2376
2377 Int_t n = p->GetNumberOfTrackReferences();
2378 for (Int_t i=0; i<n; ++i) {
2379 AliTrackReference *ref = p->GetTrackReference(i);
bc8a45bd 2380 if (!ref)
2381 continue;
56fd6cb2 2382 ref->SetUserId(ref->DetectorId());
2383 ref->Print();
2384 }
38727e64 2385}
2386
807016ea 2387//________________________________________________________________________
2388void AliAnalysisTaskEMCALPi0PbPb::ProcessDaughters(AliVParticle *p, Int_t index, const TObjArray *arr)
2389{
2390 // Process and create daughters.
2391
2392 if (!p || !arr)
2393 return;
2394
2395 AliAODMCParticle *amc = dynamic_cast<AliAODMCParticle*>(p);
2396 if (!amc)
2397 return;
2398
2399 //amc->Print();
2400
2401 Int_t nparts = arr->GetEntries();
2402 Int_t nents = fMcParts->GetEntries();
2403
2404 AliStaPart *newp = static_cast<AliStaPart*>(fMcParts->New(nents));
2405 newp->fPt = amc->Pt();
2406 newp->fEta = amc->Eta();
2407 newp->fPhi = amc->Phi();
2408 if (amc->Xv() != 0 || amc->Yv() != 0 || amc->Zv() != 0) {
2409 TVector3 vec(amc->Xv(),amc->Yv(),amc->Zv());
2410 newp->fVR = vec.Perp();
2411 newp->fVEta = vec.Eta();
2412 newp->fVPhi = vec.Phi();
2413 }
2414 newp->fPid = amc->PdgCode();
8c56d760 2415 newp->fLab = nents;
807016ea 2416 Int_t moi = amc->GetMother();
2417 if (moi>=0&&moi<nparts) {
2418 const AliAODMCParticle *mmc = static_cast<const AliAODMCParticle*>(arr->At(moi));
2419 moi = mmc->GetUniqueID();
2420 }
2421 newp->fMo = moi;
2422 p->SetUniqueID(nents);
8c56d760 2423
2424 // TODO: Determine which detector was hit
2425 //newp->fDet = ???
2426
807016ea 2427 Int_t n = amc->GetNDaughters();
2428 for (Int_t i=0; i<n; ++i) {
2429 Int_t d = amc->GetDaughter(i);
2430 if (d<=index || d>=nparts)
2431 continue;
2432 AliVParticle *dmc = static_cast<AliVParticle*>(arr->At(d));
2433 ProcessDaughters(dmc,d,arr);
2434 }
2435}
2436
2437//________________________________________________________________________
2438void AliAnalysisTaskEMCALPi0PbPb::ProcessDaughters(AliMCParticle *p, Int_t index, const AliMCEvent *arr)
2439{
2440 // Process and create daughters.
2441
2442 if (!p || !arr)
2443 return;
2444
2445 Int_t d1 = p->GetFirstDaughter();
2446 Int_t d2 = p->GetLastDaughter();
2447 if (0) {
2448 printf("%d pid=%d: %.3f %.3f %.3f (%.2f %.2f %.2f); nd=%d,%d, mo=%d\n",
2449 index,p->PdgCode(),p->Px(),p->Py(),p->Pz(),p->Xv(),p->Yv(),p->Zv(),d1,d2, p->GetMother());
bc8a45bd 2450 PrintTrackRefs(p);
807016ea 2451 }
2452 Int_t nents = fMcParts->GetEntries();
2453
2454 AliStaPart *newp = static_cast<AliStaPart*>(fMcParts->New(nents));
2455 newp->fPt = p->Pt();
2456 newp->fEta = p->Eta();
2457 newp->fPhi = p->Phi();
2458 if (p->Xv() != 0 || p->Yv() != 0 || p->Zv() != 0) {
2459 TVector3 vec(p->Xv(),p->Yv(),p->Zv());
2460 newp->fVR = vec.Perp();
2461 newp->fVEta = vec.Eta();
2462 newp->fVPhi = vec.Phi();
2463 }
2464 newp->fPid = p->PdgCode();
8c56d760 2465 newp->fLab = nents;
807016ea 2466 Int_t moi = p->GetMother();
2467 if (moi>=0) {
2468 const AliMCParticle *mmc = static_cast<const AliMCParticle *>(arr->GetTrack(moi));
2469 moi = mmc->GetUniqueID();
2470 }
2471 newp->fMo = moi;
2472 p->SetUniqueID(nents);
2473
2474 Int_t nref = p->GetNumberOfTrackReferences();
2475 if (nref>0) {
2476 AliTrackReference *ref = p->GetTrackReference(nref-1);
2477 if (ref) {
2478 newp->fDet = ref->DetectorId();
2479 }
2480 }
2481
2482 if (d1<0)
2483 return;
2484 if (d2<0)
2485 d2=d1;
2486 for (Int_t i=d1;i<=d2;++i) {
2487 AliMCParticle *dmc = static_cast<AliMCParticle *>(arr->GetTrack(i));
2488 if (dmc->P()<0.01)
2489 continue;
2490 ProcessDaughters(dmc,i,arr);
2491 }
2492}
6dbf6b27 2493
2494//__________________________________________________________________________________________________
2495void AliStaCluster::GetMom(TLorentzVector& p, Double_t *vertex)
2496{
2497 // Calculate momentum.
2498
2499 TVector3 pos;
2500 pos.SetPtEtaPhi(fR,fEta,fPhi);
2501
a186e444 2502 if(vertex){ //calculate direction relative to vertex
6dbf6b27 2503 pos -= vertex;
2504 }
2505
2506 Double_t r = pos.Mag();
2507 p.SetPxPyPzE(fE*pos.x()/r, fE*pos.y()/r, fE*pos.z()/r, fE);
2508}
2509
2510//__________________________________________________________________________________________________
2511void AliStaCluster::GetMom(TLorentzVector& p, AliStaVertex *vertex)
2512{
2513 // Calculate momentum.
2514
2515 Double_t v[3] = {0,0,0};
2516 if (vertex) {
2517 v[0] = vertex->fVx;
2518 v[1] = vertex->fVy;
2519 v[2] = vertex->fVz;
2520 }
2521 GetMom(p, v);
2522}