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