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