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