]> 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
in internal trigger selection, no GetTriggerMask unless ( kNoSelection == fInternalTr...
[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) {
60d77596 1359 Short_t id=-1;
1360 Int_t mclabel = -1;
77e93dc2 1361 Double_t amp=0,time=0, efrac = 0;
1362 if (!cells->GetCell(i, id, amp, time,mclabel,efrac))
323834f0 1363 continue;
1364 amp -= cellMeanE;
1365 if (amp<0.001)
1366 amp = 0;
1367 cells->SetCell(i, id, amp, time);
1368 }
1369
1370 TObjArray *clusters = fEsdClusters;
1371 if (!clusters)
1372 clusters = fAodClusters;
323834f0 1373 if (!clusters)
1374 return;
1375
1376 Int_t nclus = clusters->GetEntries();
1377 for (Int_t i = 0; i<nclus; ++i) {
1378 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1379 if (!clus->IsEMCAL())
1380 continue;
1381 Int_t nc = clus->GetNCells();
1382 Double_t clusE = 0;
1383 UShort_t ids[100] = {0};
1384 Double_t fra[100] = {0};
1385 for (Int_t j = 0; j<nc; ++j) {
1386 Short_t id = TMath::Abs(clus->GetCellAbsId(j));
1387 Double_t cen = cells->GetCellAmplitude(id);
1388 clusE += cen;
1389 if (cen>0) {
1390 ids[nc] = id;
1391 ++nc;
1392 }
1393 }
1394 if (clusE<=0) {
1395 clusters->RemoveAt(i);
1396 continue;
1397 }
1398
1399 for (Int_t j = 0; j<nc; ++j) {
1400 Short_t id = ids[j];
1401 Double_t cen = cells->GetCellAmplitude(id);
1402 fra[j] = cen/clusE;
1403 }
1404 clus->SetE(clusE);
1405 AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus);
1406 if (aodclus) {
1407 aodclus->Clear("");
1408 aodclus->SetNCells(nc);
1409 aodclus->SetCellsAmplitudeFraction(fra);
1410 aodclus->SetCellsAbsId(ids);
1411 }
1412 }
1413 clusters->Compress();
1414}
1415
286b47a5 1416//________________________________________________________________________
1417void AliAnalysisTaskEMCALPi0PbPb::FillCellHists()
1418{
1419 // Fill histograms related to cell properties.
d595acbb 1420
90d5b88b 1421 AliVCaloCells *cells = fEsdCells;
1422 if (!cells)
1423 cells = fAodCells;
1424
296ea9b4 1425 if (!cells)
1426 return;
90d5b88b 1427
296ea9b4 1428 Int_t cellModCount[12] = {0};
1429 Double_t cellMaxE = 0;
1430 Double_t cellMeanE = 0;
1431 Int_t ncells = cells->GetNumberOfCells();
1432 if (ncells==0)
1433 return;
1434
1435 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
1436
1437 for (Int_t i = 0; i<ncells; ++i) {
1438 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
1439 Double_t cellE = cells->GetAmplitude(i);
1440 fHCellE->Fill(cellE);
1441 if (cellE>cellMaxE)
1442 cellMaxE = cellE;
1443 cellMeanE += cellE;
1444
1445 Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1;
1446 Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta);
1447 if (!ret) {
1448 AliError(Form("Could not get cell index for %d", absID));
1449 continue;
1450 }
1451 ++cellModCount[iSM];
1452 Int_t iPhi=-1, iEta=-1;
1453 fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta);
1454 fHColuRow[iSM]->Fill(iEta,iPhi,1);
1455 fHColuRowE[iSM]->Fill(iEta,iPhi,cellE);
1456 fHCellFreqNoCut[iSM]->Fill(absID);
2e4d8148 1457 if (cellE > 0.1) fHCellFreqCut100M[iSM]->Fill(absID);
1458 if (cellE > 0.3) fHCellFreqCut300M[iSM]->Fill(absID);
296ea9b4 1459 if (fHCellCheckE && fHCellCheckE[absID])
1460 fHCellCheckE[absID]->Fill(cellE);
2e4d8148 1461 fHCellFreqE[iSM]->Fill(absID, cellE);
296ea9b4 1462 }
1463 fHCellH->Fill(cellMaxE);
1464 cellMeanE /= ncells;
1465 fHCellM->Fill(cellMeanE);
1466 fHCellM2->Fill(cellMeanE*ncells/24/48/nsm);
1467 for (Int_t i=0; i<nsm; ++i)
1468 fHCellMult[i]->Fill(cellModCount[i]);
286b47a5 1469}
1470
1471//________________________________________________________________________
1472void AliAnalysisTaskEMCALPi0PbPb::FillClusHists()
1473{
90d5b88b 1474 // Fill histograms related to cluster properties.
fa443410 1475
90d5b88b 1476 TObjArray *clusters = fEsdClusters;
1477 if (!clusters)
1478 clusters = fAodClusters;
296ea9b4 1479 if (!clusters)
1480 return;
90d5b88b 1481
296ea9b4 1482 Double_t vertex[3] = {0};
1483 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1484
1485 Int_t nclus = clusters->GetEntries();
1486 for(Int_t i = 0; i<nclus; ++i) {
1487 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1488 if (!clus)
1489 continue;
1490 if (!clus->IsEMCAL())
1491 continue;
90d5b88b 1492 TLorentzVector clusterVec;
296ea9b4 1493 clus->GetMomentum(clusterVec,vertex);
3a952328 1494 Double_t maxAxis = clus->GetTOF(); //sigma
1495 Double_t clusterEcc = clus->Chi2(); //eccentricity
296ea9b4 1496 fHClustEccentricity->Fill(clusterEcc);
1497 fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi());
1498 fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt());
1499 fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E());
1500 fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis);
1501 fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E());
f5e0f1e2 1502 fHClustEnergyNCell->Fill(clus->E(),clus->GetNCells());
788ca675 1503 }
1504}
b3ee6797 1505
38727e64 1506//________________________________________________________________________
1507void AliAnalysisTaskEMCALPi0PbPb::CalcMcInfo()
1508{
1509 // Get Mc truth particle information.
1510
cfd7d5b2 1511 if (!fMcParts)
38727e64 1512 return;
1513
807016ea 1514 fMcParts->Clear();
1515
1516 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
1517 Double_t etamin = -0.7;
1518 Double_t etamax = +0.7;
1519 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
1520 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
1521
56fd6cb2 1522 if (fAodEv) {
807016ea 1523 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1524 am->LoadBranch(AliAODMCParticle::StdBranchName());
56fd6cb2 1525 TClonesArray *tca = dynamic_cast<TClonesArray*>(fAodEv->FindListObject(AliAODMCParticle::StdBranchName()));
1526 if (!tca)
1527 return;
1528
1529 Int_t nents = tca->GetEntries();
807016ea 1530 for(int it=0; it<nents; ++it) {
56fd6cb2 1531 AliAODMCParticle *part = static_cast<AliAODMCParticle*>(tca->At(it));
807016ea 1532 part->Print();
56fd6cb2 1533
b5ba4932 1534 // pion or eta meson or direct photon
56fd6cb2 1535 if(part->GetPdgCode() == 111) {
1536 } else if(part->GetPdgCode() == 221) {
b5ba4932 1537 } else if(part->GetPdgCode() == 22 ) {
1538 } else
56fd6cb2 1539 continue;
1540
1541 // primary particle
1542 Double_t dR = TMath::Sqrt((part->Xv()*part->Xv())+(part->Yv()*part->Yv()));
1543 if(dR > 1.0)
1544 continue;
56fd6cb2 1545
807016ea 1546 // kinematic cuts
1547 Double_t pt = part->Pt() ;
1548 if (pt<0.5)
1549 continue;
1550 Double_t eta = part->Eta();
1551 if (eta<etamin||eta>etamax)
1552 continue;
1553 Double_t phi = part->Phi();
1554 if (phi<phimin||phi>phimax)
1555 continue;
1556
1557 ProcessDaughters(part, it, tca);
56fd6cb2 1558 }
807016ea 1559 return;
1560 }
38727e64 1561
1562 AliMCEvent *mcEvent = MCEvent();
1563 if (!mcEvent)
1564 return;
1565
1566 const AliVVertex *evtVtx = mcEvent->GetPrimaryVertex();
1567 if (!evtVtx)
1568 return;
1569
1570 mcEvent->PreReadAll();
38727e64 1571
1572 Int_t nTracks = mcEvent->GetNumberOfPrimaries();
1573 for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
1574 AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
1575 if (!mcP)
1576 continue;
1577
b5ba4932 1578 // pion or eta meson or direct photon
807016ea 1579 if(mcP->PdgCode() == 111) {
1580 } else if(mcP->PdgCode() == 221) {
b5ba4932 1581 } else if(mcP->PdgCode() == 22 ) {
38727e64 1582 } else
1583 continue;
1584
1585 // primary particle
807016ea 1586 Double_t dR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
1587 (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
38727e64 1588 if(dR > 1.0)
1589 continue;
1590
38727e64 1591 // kinematic cuts
807016ea 1592 Double_t pt = mcP->Pt() ;
1593 if (pt<0.5)
38727e64 1594 continue;
807016ea 1595 Double_t eta = mcP->Eta();
1596 if (eta<etamin||eta>etamax)
38727e64 1597 continue;
807016ea 1598 Double_t phi = mcP->Phi();
1599 if (phi<phimin||phi>phimax)
1600 continue;
1601
1602 ProcessDaughters(mcP, iTrack, mcEvent);
38727e64 1603 }
1604}
1605
788ca675 1606//________________________________________________________________________
1607void AliAnalysisTaskEMCALPi0PbPb::FillNtuple()
1608{
1609 // Fill ntuple.
1610
1611 if (!fNtuple)
1612 return;
1613
1614 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1615 if (fAodEv) {
1616 fHeader->fRun = fAodEv->GetRunNumber();
1617 fHeader->fOrbit = fAodEv->GetHeader()->GetOrbitNumber();
1618 fHeader->fPeriod = fAodEv->GetHeader()->GetPeriodNumber();
1619 fHeader->fBx = fAodEv->GetHeader()->GetBunchCrossNumber();
1620 fHeader->fL0 = fAodEv->GetHeader()->GetL0TriggerInputs();
1621 fHeader->fL1 = fAodEv->GetHeader()->GetL1TriggerInputs();
1622 fHeader->fL2 = fAodEv->GetHeader()->GetL2TriggerInputs();
1623 fHeader->fTrClassMask = fAodEv->GetHeader()->GetTriggerMask();
1624 fHeader->fTrCluster = fAodEv->GetHeader()->GetTriggerCluster();
1625 fHeader->fOffTriggers = fAodEv->GetHeader()->GetOfflineTrigger();
1626 fHeader->fFiredTriggers = fAodEv->GetHeader()->GetFiredTriggerClasses();
1627 } else {
1628 fHeader->fRun = fEsdEv->GetRunNumber();
1629 fHeader->fOrbit = fEsdEv->GetHeader()->GetOrbitNumber();
1630 fHeader->fPeriod = fEsdEv->GetESDRun()->GetPeriodNumber();
1631 fHeader->fBx = fEsdEv->GetHeader()->GetBunchCrossNumber();
1632 fHeader->fL0 = fEsdEv->GetHeader()->GetL0TriggerInputs();
1633 fHeader->fL1 = fEsdEv->GetHeader()->GetL1TriggerInputs();
1634 fHeader->fL2 = fEsdEv->GetHeader()->GetL2TriggerInputs();
1635 fHeader->fTrClassMask = fEsdEv->GetHeader()->GetTriggerMask();
1636 fHeader->fTrCluster = fEsdEv->GetHeader()->GetTriggerCluster();
1637 fHeader->fOffTriggers = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
1638 fHeader->fFiredTriggers = fEsdEv->GetFiredTriggerClasses();
5fe1ca23 1639 Float_t v0CorrR = 0;
1640 fHeader->fV0 = AliESDUtils::GetCorrV0(fEsdEv,v0CorrR);
1641 const AliMultiplicity *mult = fEsdEv->GetMultiplicity();
1642 if (mult)
1643 fHeader->fCl1 = mult->GetNumberOfITSClusters(1);
1644 fHeader->fTr = AliESDtrackCuts::GetReferenceMultiplicity(fEsdEv,1);
c0d2e1f2 1645 AliTriggerAnalysis trAn; /// Trigger Analysis
1646 Bool_t v0B = trAn.IsOfflineTriggerFired(fEsdEv, AliTriggerAnalysis::kV0C);
1647 Bool_t v0A = trAn.IsOfflineTriggerFired(fEsdEv, AliTriggerAnalysis::kV0A);
469b2bff 1648 fHeader->fV0And = v0A && v0B;
1649 fHeader->fIsHT = (fHeader->fOffTriggers & AliVEvent::kEMC1) || (fHeader->fOffTriggers & AliVEvent::kEMC7);
1650 am->LoadBranch("SPDPileupVertices");
1651 am->LoadBranch("TrkPileupVertices");
1652 fHeader->fIsPileup = fEsdEv->IsPileupFromSPD(3,0.8);
1653 fHeader->fIsPileup2 = fEsdEv->IsPileupFromSPD(3,0.4);
1654 fHeader->fIsPileup4 = fEsdEv->IsPileupFromSPD(3,0.2);
1655 fHeader->fIsPileup8 = fEsdEv->IsPileupFromSPD(3,0.1);
1656 fHeader->fNSpdVertices = fEsdEv->GetNumberOfPileupVerticesSPD();
1657 fHeader->fNTpcVertices = fEsdEv->GetNumberOfPileupVerticesTracks();
788ca675 1658 }
2ee7bda4 1659
788ca675 1660 AliCentrality *cent = InputEvent()->GetCentrality();
1661 fHeader->fV0Cent = cent->GetCentralityPercentileUnchecked("V0M");
1662 fHeader->fCl1Cent = cent->GetCentralityPercentileUnchecked("CL1");
1663 fHeader->fTrCent = cent->GetCentralityPercentileUnchecked("TRK");
1664 fHeader->fCqual = cent->GetQuality();
b6c599fe 1665
1666 AliEventplane *ep = InputEvent()->GetEventplane();
1667 if (ep) {
1668 if (ep->GetQVector())
1669 fHeader->fPsi = ep->GetQVector()->Phi()/2. ;
f2b8fca6 1670 else
1671 fHeader->fPsi = -1;
b6c599fe 1672 if (ep->GetQsub1()&&ep->GetQsub2())
1673 fHeader->fPsiRes = ep->GetQsub1()->Phi()/2.-ep->GetQsub2()->Phi()/2.;
1674 else
1675 fHeader->fPsiRes = 0;
1676 }
1677
788ca675 1678 Double_t val = 0;
1679 TString trgclasses(fHeader->fFiredTriggers);
1680 for (Int_t j = 0; j<fTrClassNamesArr->GetEntries(); ++j) {
1681 const char *name = fTrClassNamesArr->At(j)->GetName();
2ee7bda4 1682 TRegexp regexp(name);
1683 if (trgclasses.Contains(regexp))
788ca675 1684 val += TMath::Power(2,j);
1685 }
1686 fHeader->fTcls = (UInt_t)val;
1687
5fe1ca23 1688 fHeader->fNSelTr = fSelTracks->GetEntries();
1689 fHeader->fNSelPrimTr = fSelPrimTracks->GetEntries();
f5e0f1e2 1690 fHeader->fNSelPrimTr1 = 0;
1691 fHeader->fNSelPrimTr2 = 0;
1692 for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); iTracks++){
1693 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1694 if(track->Pt()>1)
1695 ++fHeader->fNSelPrimTr1;
1696 if(track->Pt()>2)
1697 ++fHeader->fNSelPrimTr2;
1698 }
5fe1ca23 1699
1700 fHeader->fNCells = 0;
c0d2e1f2 1701 fHeader->fNCells0 = 0;
1702 fHeader->fNCells01 = 0;
1703 fHeader->fNCells03 = 0;
5fe1ca23 1704 fHeader->fNCells1 = 0;
1705 fHeader->fNCells2 = 0;
1706 fHeader->fNCells5 = 0;
1707 fHeader->fMaxCellE = 0;
1708
1709 AliVCaloCells *cells = fEsdCells;
1710 if (!cells)
1711 cells = fAodCells;
1712
1713 if (cells) {
1714 Int_t ncells = cells->GetNumberOfCells();
1715 for(Int_t j=0; j<ncells; ++j) {
1716 Double_t cellen = cells->GetAmplitude(j);
c0d2e1f2 1717 if (cellen>0.045)
1718 ++fHeader->fNCells0;
1719 if (cellen>0.1)
1720 ++fHeader->fNCells01;
1721 if (cellen>0.3)
1722 ++fHeader->fNCells03;
5fe1ca23 1723 if (cellen>1)
1724 ++fHeader->fNCells1;
1725 if (cellen>2)
1726 ++fHeader->fNCells2;
1727 if (cellen>5)
1728 ++fHeader->fNCells5;
1729 if (cellen>fHeader->fMaxCellE)
1730 fHeader->fMaxCellE = cellen;
1731 }
1732 fHeader->fNCells = ncells;
1733 }
1734
1735 fHeader->fNClus = 0;
1736 fHeader->fNClus1 = 0;
1737 fHeader->fNClus2 = 0;
1738 fHeader->fNClus5 = 0;
1739 fHeader->fMaxClusE = 0;
1740
1741 TObjArray *clusters = fEsdClusters;
1742 if (!clusters)
1743 clusters = fAodClusters;
1744
1745 if (clusters) {
1746 Int_t nclus = clusters->GetEntries();
1747 for(Int_t j=0; j<nclus; ++j) {
1748 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(j));
6d4faab0 1749 if (!clus->IsEMCAL())
1750 continue;
5fe1ca23 1751 Double_t clusen = clus->E();
1752 if (clusen>1)
1753 ++fHeader->fNClus1;
1754 if (clusen>2)
1755 ++fHeader->fNClus2;
1756 if (clusen>5)
1757 ++fHeader->fNClus5;
1758 if (clusen>fHeader->fMaxClusE)
1759 fHeader->fMaxClusE = clusen;
1760 }
1761 fHeader->fNClus = nclus;
1762 }
1763
2ee7bda4 1764 fHeader->fMaxTrE = 0;
1765 if (fTriggers) {
1766 Int_t ntrig = fTriggers->GetEntries();
1767 for (Int_t j = 0; j<ntrig; ++j) {
1768 AliStaTrigger *sta = static_cast<AliStaTrigger*>(fTriggers->At(j));
1769 if (!sta)
1770 continue;
1771 if (sta->fE>fHeader->fMaxTrE)
1772 fHeader->fMaxTrE = sta->fE;
1773 }
1774 }
1775
1776 // count cells above 100 MeV on super modules
1777 fHeader->fNcSM0 = GetNCells(0, 0.100);
1778 fHeader->fNcSM1 = GetNCells(1, 0.100);
1779 fHeader->fNcSM2 = GetNCells(2, 0.100);
1780 fHeader->fNcSM3 = GetNCells(3, 0.100);
1781 fHeader->fNcSM4 = GetNCells(4, 0.100);
1782 fHeader->fNcSM5 = GetNCells(5, 0.100);
1783 fHeader->fNcSM6 = GetNCells(6, 0.100);
1784 fHeader->fNcSM7 = GetNCells(7, 0.100);
1785 fHeader->fNcSM8 = GetNCells(8, 0.100);
1786 fHeader->fNcSM9 = GetNCells(9, 0.100);
1787
788ca675 1788 if (fAodEv) {
1789 am->LoadBranch("vertices");
1790 AliAODVertex *pv = fAodEv->GetPrimaryVertex();
1791 FillVertex(fPrimVert, pv);
1792 AliAODVertex *sv = fAodEv->GetPrimaryVertexSPD();
1793 FillVertex(fSpdVert, sv);
1794 } else {
1795 am->LoadBranch("PrimaryVertex.");
1796 const AliESDVertex *pv = fEsdEv->GetPrimaryVertexTracks();
1797 FillVertex(fPrimVert, pv);
1798 am->LoadBranch("SPDVertex.");
1799 const AliESDVertex *sv = fEsdEv->GetPrimaryVertexSPD();
1800 FillVertex(fSpdVert, sv);
1801 am->LoadBranch("TPCVertex.");
1802 const AliESDVertex *tv = fEsdEv->GetPrimaryVertexTPC();
1803 FillVertex(fTpcVert, tv);
fa443410 1804 }
788ca675 1805
788ca675 1806 fNtuple->Fill();
286b47a5 1807}
1808
1809//________________________________________________________________________
1810void AliAnalysisTaskEMCALPi0PbPb::FillPionHists()
1811{
1812 // Fill histograms related to pions.
286b47a5 1813
296ea9b4 1814 Double_t vertex[3] = {0};
fa443410 1815 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1816
90d5b88b 1817 TObjArray *clusters = fEsdClusters;
1818 if (!clusters)
1819 clusters = fAodClusters;
fa443410 1820
90d5b88b 1821 if (clusters) {
1822 TLorentzVector clusterVec1;
1823 TLorentzVector clusterVec2;
1824 TLorentzVector pionVec;
1825
1826 Int_t nclus = clusters->GetEntries();
fa443410 1827 for (Int_t i = 0; i<nclus; ++i) {
90d5b88b 1828 AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i));
fa443410 1829 if (!clus1)
1830 continue;
1831 if (!clus1->IsEMCAL())
1832 continue;
3c661da5 1833 if (clus1->E()<fMinE)
6eb6260e 1834 continue;
a49742b5 1835 if (clus1->GetNCells()<fNminCells)
1836 continue;
f224d35b 1837 if (GetMaxCellEnergy(clus1)/clus1->E()<fMinErat)
1838 continue;
3c661da5 1839 if (clus1->Chi2()<fMinEcc) // eccentricity cut
f224d35b 1840 continue;
fa443410 1841 clus1->GetMomentum(clusterVec1,vertex);
1842 for (Int_t j = i+1; j<nclus; ++j) {
90d5b88b 1843 AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j));
fa443410 1844 if (!clus2)
1845 continue;
1846 if (!clus2->IsEMCAL())
1847 continue;
3c661da5 1848 if (clus2->E()<fMinE)
6eb6260e 1849 continue;
a49742b5 1850 if (clus2->GetNCells()<fNminCells)
1851 continue;
f224d35b 1852 if (GetMaxCellEnergy(clus2)/clus2->E()<fMinErat)
1853 continue;
3c661da5 1854 if (clus2->Chi2()<fMinEcc) // eccentricity cut
f224d35b 1855 continue;
fa443410 1856 clus2->GetMomentum(clusterVec2,vertex);
1857 pionVec = clusterVec1 + clusterVec2;
1858 Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E();
6eb6260e 1859 Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect());
d595acbb 1860 if (pionZgg < fAsymMax) {
1861 fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi());
1862 fHPionMggPt->Fill(pionVec.M(),pionVec.Pt());
1863 fHPionMggAsym->Fill(pionVec.M(),pionZgg);
6eb6260e 1864 fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle);
d595acbb 1865 Int_t bin = fPtRanges->FindBin(pionVec.Pt());
1866 fHPionInvMasses[bin]->Fill(pionVec.M());
1867 }
fa443410 1868 }
1869 }
90d5b88b 1870 }
fa443410 1871}
d595acbb 1872
38727e64 1873//________________________________________________________________________
1874void AliAnalysisTaskEMCALPi0PbPb::FillMcHists()
1875{
1876 // Fill additional MC information histograms.
1877
cfd7d5b2 1878 if (!fMcParts)
38727e64 1879 return;
1880
1881 // check if aod or esd mc mode and the fill histos
1882}
1883
6eb6260e 1884//________________________________________________________________________
323834f0 1885void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists()
6eb6260e 1886{
788ca675 1887 // Fill other histograms.
b5ba4932 1888
1889 for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); ++iTracks){
1890 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1891 if(!track)
1892 continue;
1893 fHPrimTrackPt->Fill(track->Pt());
1894 fHPrimTrackEta->Fill(track->Eta());
1895 fHPrimTrackPhi->Fill(track->Phi());
1896 }
788ca675 1897}
1898
f5e0f1e2 1899//________________________________________________________________________
1900void AliAnalysisTaskEMCALPi0PbPb::FillTrackHists()
1901{
1902 // Fill track histograms.
1903
1904 if (fSelPrimTracks) {
1905 for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); iTracks++) {
1906 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1907 if(!track)
1908 continue;
1909 fHPrimTrackPt->Fill(track->Pt());
1910 fHPrimTrackEta->Fill(track->Eta());
1911 fHPrimTrackPhi->Fill(track->Phi());
1912 }
1913 }
1914}
1915
788ca675 1916//__________________________________________________________________________________________________
1917void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliESDVertex *esdv)
1918{
1919 // Fill vertex from ESD vertex info.
1920
1921 v->fVx = esdv->GetX();
1922 v->fVy = esdv->GetY();
1923 v->fVz = esdv->GetZ();
1924 v->fVc = esdv->GetNContributors();
1925 v->fDisp = esdv->GetDispersion();
1926 v->fZres = esdv->GetZRes();
1927 v->fChi2 = esdv->GetChi2();
1928 v->fSt = esdv->GetStatus();
1929 v->fIs3D = esdv->IsFromVertexer3D();
1930 v->fIsZ = esdv->IsFromVertexerZ();
1931}
1932
1933//__________________________________________________________________________________________________
1934void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliAODVertex *aodv)
1935{
1936 // Fill vertex from AOD vertex info.
1937
1938 v->fVx = aodv->GetX();
1939 v->fVy = aodv->GetY();
1940 v->fVz = aodv->GetZ();
1941 v->fVc = aodv->GetNContributors();
1942 v->fChi2 = aodv->GetChi2();
6eb6260e 1943}
1944
d595acbb 1945//________________________________________________________________________
296ea9b4 1946Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
1947{
1948 // Compute isolation based on cell content.
1949
1950 AliVCaloCells *cells = fEsdCells;
1951 if (!cells)
1952 cells = fAodCells;
1953 if (!cells)
1954 return 0;
1955
1956 Double_t cellIsolation = 0;
1957 Double_t rad2 = radius*radius;
1958 Int_t ncells = cells->GetNumberOfCells();
1959 for (Int_t i = 0; i<ncells; ++i) {
1960 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
0fbe8d4f 1961 Float_t eta=-1, phi=-1;
296ea9b4 1962 fGeom->EtaPhiFromIndex(absID,eta,phi);
0fbe8d4f 1963 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
296ea9b4 1964 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1965 if(dist>rad2)
1966 continue;
0fbe8d4f 1967 Double_t cellE = cells->GetAmplitude(i);
2ee7bda4 1968 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
1969 Double_t cellEt = cellE*sin(theta);
1970 cellIsolation += cellEt;
296ea9b4 1971 }
1972 return cellIsolation;
1973}
2ee7bda4 1974
5fc7508c 1975//________________________________________________________________________
1976Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsoNxM(Double_t cEta, Double_t cPhi, Int_t N, Int_t M) const
1977{
1978 // Compute isolation based on cell content, in a NxM rectangle.
1979
1980 AliVCaloCells *cells = fEsdCells;
1981 if (!cells)
1982 cells = fAodCells;
1983 if (!cells)
1984 return 0;
1985
1986 Double_t cellIsolation = 0;
1987 Int_t ncells = cells->GetNumberOfCells();
1988 for (Int_t i = 0; i<ncells; ++i) {
1989 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
1990 Float_t eta=-1, phi=-1;
1991 fGeom->EtaPhiFromIndex(absID,eta,phi);
1992 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
1993 Double_t etadiff = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1994 if(TMath::Abs(etadiff)/0.014>N)
1995 continue;
1996 if(TMath::Abs(phidiff)/0.014>M)
1997 continue;
1998 Double_t cellE = cells->GetAmplitude(i);
2ee7bda4 1999 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
2000 Double_t cellEt = cellE*sin(theta);
2001 cellIsolation += cellEt;
5fc7508c 2002 }
2003 return cellIsolation;
2004}
2005
296ea9b4 2006//________________________________________________________________________
0fbe8d4f 2007Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellEnergy(const AliVCluster *cluster) const
2008{
2009 // Get maximum energy of attached cell.
2010
2011 Double_t ret = 0;
2012 Int_t ncells = cluster->GetNCells();
2013 if (fEsdCells) {
2014 for (Int_t i=0; i<ncells; i++) {
2015 Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2016 ret += e;
2017 }
2018 } else {
2019 for (Int_t i=0; i<ncells; i++) {
2020 Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2021 ret += e;
2022 }
2023 }
2024 return ret;
2025}
2026
2027//________________________________________________________________________
2028Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
d595acbb 2029{
90d5b88b 2030 // Get maximum energy of attached cell.
2031
788ca675 2032 id = -1;
90d5b88b 2033 Double_t maxe = 0;
90d5b88b 2034 Int_t ncells = cluster->GetNCells();
2035 if (fEsdCells) {
2036 for (Int_t i=0; i<ncells; i++) {
27422847 2037 Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
f0897c18 2038 if (e>maxe) {
90d5b88b 2039 maxe = e;
788ca675 2040 id = cluster->GetCellAbsId(i);
f0897c18 2041 }
90d5b88b 2042 }
2043 } else {
2044 for (Int_t i=0; i<ncells; i++) {
27422847 2045 Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
90d5b88b 2046 if (e>maxe)
2047 maxe = e;
788ca675 2048 id = cluster->GetCellAbsId(i);
90d5b88b 2049 }
2050 }
6eb6260e 2051 return maxe;
d595acbb 2052}
2053
469b2bff 2054//________________________________________________________________________
1f41ed3d 2055Double_t AliAnalysisTaskEMCALPi0PbPb::GetSecondMaxCellEnergy(AliVCluster *clus, Short_t &id) const
469b2bff 2056{
2057 // Get second maximum cell.
2058
2059 AliVCaloCells *cells = fEsdCells;
2060 if (!cells)
2061 cells = fAodCells;
2062 if (!cells)
2063 return -1;
2064
2065 Double_t secondEmax=0, firstEmax=0;
2066 Double_t cellen;
2067 for(Int_t iCell=0;iCell<clus->GetNCells();iCell++){
2068 Int_t absId = clus->GetCellAbsId(iCell);
2069 cellen = cells->GetCellAmplitude(absId);
2070 if(cellen > firstEmax)
2071 firstEmax = cellen;
2072 }
2073 for(Int_t iCell=0;iCell<clus->GetNCells();iCell++){
2074 Int_t absId = clus->GetCellAbsId(iCell);
2075 cellen = cells->GetCellAmplitude(absId);
1f41ed3d 2076 if(cellen < firstEmax && cellen > secondEmax) {
469b2bff 2077 secondEmax = cellen;
1f41ed3d 2078 id = absId;
2079 }
469b2bff 2080 }
2081 return secondEmax;
2082}
2083
d595acbb 2084//________________________________________________________________________
0fbe8d4f 2085void AliAnalysisTaskEMCALPi0PbPb::GetSigma(const AliVCluster *c, Double_t& sigmaMax, Double_t &sigmaMin) const
d595acbb 2086{
90d5b88b 2087 // Calculate the (E) weighted variance along the longer (eigen) axis.
2088
6bf90832 2089 sigmaMax = 0; // cluster variance along its longer axis
2090 sigmaMin = 0; // cluster variance along its shorter axis
2091 Double_t Ec = c->E(); // cluster energy
296ea9b4 2092 if(Ec<=0)
2093 return;
6bf90832 2094 Double_t Xc = 0 ; // cluster first moment along X
2095 Double_t Yc = 0 ; // cluster first moment along Y
2096 Double_t Sxx = 0 ; // cluster second central moment along X (variance_X^2)
2097 Double_t Sxy = 0 ; // cluster second central moment along Y (variance_Y^2)
2098 Double_t Syy = 0 ; // cluster covariance^2
90d5b88b 2099
2100 AliVCaloCells *cells = fEsdCells;
2101 if (!cells)
2102 cells = fAodCells;
2103
6eb6260e 2104 if (!cells)
2105 return;
2106
6bf90832 2107 Int_t ncells = c->GetNCells();
6eb6260e 2108 if (ncells==1)
2109 return;
2110
6eb6260e 2111 for(Int_t j=0; j<ncells; ++j) {
6bf90832 2112 Int_t id = TMath::Abs(c->GetCellAbsId(j));
6eb6260e 2113 Double_t cellen = cells->GetCellAmplitude(id);
3a952328 2114 TVector3 pos;
2115 fGeom->GetGlobal(id,pos);
6eb6260e 2116 Xc += cellen*pos.X();
2117 Yc += cellen*pos.Y();
2118 Sxx += cellen*pos.X()*pos.X();
2119 Syy += cellen*pos.Y()*pos.Y();
2120 Sxy += cellen*pos.X()*pos.Y();
2121 }
2122 Xc /= Ec;
2123 Yc /= Ec;
2124 Sxx /= Ec;
2125 Syy /= Ec;
2126 Sxy /= Ec;
2127 Sxx -= Xc*Xc;
2128 Syy -= Yc*Yc;
2129 Sxy -= Xc*Yc;
f0897c18 2130 Sxx = TMath::Abs(Sxx);
2131 Syy = TMath::Abs(Syy);
296ea9b4 2132 sigmaMax = (Sxx + Syy + TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
2133 sigmaMax = TMath::Sqrt(TMath::Abs(sigmaMax));
2134 sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
2135 sigmaMin = TMath::Sqrt(TMath::Abs(sigmaMin));
d595acbb 2136}
2ee7bda4 2137
5fc7508c 2138//________________________________________________________________________
2139void AliAnalysisTaskEMCALPi0PbPb::GetSigmaEtaEta(const AliVCluster *c, Double_t& sEtaEta, Double_t &sPhiPhi) const
2140{
2141 // Calculate the (E) weighted variance along the pseudorapidity.
2142
2143 sEtaEta = 0;
2144 sPhiPhi = 0;
2145
2146 Double_t Ec = c->E(); // cluster energy
2147 if(Ec<=0)
2148 return;
2149
2150 const Int_t ncells = c->GetNCells();
2151
2152 Double_t EtaC = 0; // cluster first moment along eta
2153 Double_t PhiC = 0; // cluster first moment along phi
2154 Double_t Setaeta = 0; // cluster second central moment along eta
2155 Double_t Sphiphi = 0; // cluster second central moment along phi
2156 Double_t w[ncells]; // weight max(0,4.5*log(E_i/Ec))
2157 Double_t sumw = 0;
2158 Int_t id[ncells];
2159
2160 AliVCaloCells *cells = fEsdCells;
2161 if (!cells)
2162 cells = fAodCells;
2163
2164 if (!cells)
2165 return;
2166
2167 if (ncells==1)
2168 return;
2169
2170 for(Int_t j=0; j<ncells; ++j) {
2171 id[j] = TMath::Abs(c->GetCellAbsId(j));
2172 Double_t cellen = cells->GetCellAmplitude(id[j]);
2173 w[j] = TMath::Max(0., 4.5+TMath::Log(cellen/Ec));
2174 TVector3 pos;
2175 fGeom->GetGlobal(id[j],pos);
2176 EtaC += w[j]*pos.Eta();
2177 PhiC += w[j]*pos.Phi();
2178 sumw += w[j];
2179 }
2180 EtaC /= sumw;
2181 PhiC /= sumw;
2182
2183 for(Int_t j=0; j<ncells; ++j) {
2184 TVector3 pos;
2185 fGeom->GetGlobal(id[j],pos);
2186 Setaeta = w[j]*(pos.Eta() - EtaC)*(pos.Eta() - EtaC);
2187 Sphiphi = w[j]*(pos.Phi() - PhiC)*(pos.Phi() - PhiC);
2188 }
2189 Setaeta /= sumw;
2190 sEtaEta = TMath::Sqrt(Setaeta);
2191 Sphiphi /= sumw;
2192 sPhiPhi = TMath::Sqrt(Sphiphi);
2193}
2194
6bf90832 2195//________________________________________________________________________
0fbe8d4f 2196Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(const AliVCluster *c, Double_t emin) const
6bf90832 2197{
2198 // Calculate number of attached cells above emin.
2199
6bf90832 2200 AliVCaloCells *cells = fEsdCells;
2201 if (!cells)
2202 cells = fAodCells;
6bf90832 2203 if (!cells)
296ea9b4 2204 return 0;
6bf90832 2205
296ea9b4 2206 Int_t n = 0;
6bf90832 2207 Int_t ncells = c->GetNCells();
2208 for(Int_t j=0; j<ncells; ++j) {
2209 Int_t id = TMath::Abs(c->GetCellAbsId(j));
2210 Double_t cellen = cells->GetCellAmplitude(id);
2211 if (cellen>=emin)
2212 ++n;
2213 }
2214 return n;
2215}
2216
2ee7bda4 2217//________________________________________________________________________
2218Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(Int_t sm, Double_t emin) const
2219{
2220 // Calculate number of cells per SM above emin.
2221
2222 AliVCaloCells *cells = fEsdCells;
2223 if (!cells)
2224 cells = fAodCells;
2225 if (!cells)
2226 return 0;
2227
2228 Int_t n = 0;
2229 Int_t ncells = cells->GetNumberOfCells();
2230 for(Int_t j=0; j<ncells; ++j) {
2231 Int_t id = TMath::Abs(cells->GetCellNumber(j));
2232 Double_t cellen = cells->GetCellAmplitude(id);
2233 if (cellen<emin)
2234 continue;
2235 Int_t fsm = fGeom->GetSuperModuleNumber(id);
2236 if (fsm != sm)
2237 continue;
2238 ++n;
2239 }
2240 return n;
2241}
2242
296ea9b4 2243//________________________________________________________________________
b6c599fe 2244Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
296ea9b4 2245{
2246 // Compute isolation based on tracks.
2247
2248 Double_t trkIsolation = 0;
2249 Double_t rad2 = radius*radius;
b6c599fe 2250 Int_t ntrks = fSelPrimTracks->GetEntries();
296ea9b4 2251 for(Int_t j = 0; j<ntrks; ++j) {
1f41ed3d 2252 AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
296ea9b4 2253 if (!track)
2254 continue;
b6c599fe 2255 if (track->Pt()<pt)
2256 continue;
296ea9b4 2257 Float_t eta = track->Eta();
2258 Float_t phi = track->Phi();
0fbe8d4f 2259 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
296ea9b4 2260 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
2261 if(dist>rad2)
2262 continue;
2263 trkIsolation += track->Pt();
2264 }
2265 return trkIsolation;
2266}
2ee7bda4 2267
5fc7508c 2268//________________________________________________________________________
2269Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsoStrip(Double_t cEta, Double_t cPhi, Double_t dEta, Double_t dPhi, Double_t pt) const
2270{
2271 // Compute isolation based on tracks.
2272
2273 Double_t trkIsolation = 0;
2274 Int_t ntrks = fSelPrimTracks->GetEntries();
2275 for(Int_t j = 0; j<ntrks; ++j) {
1f41ed3d 2276 AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
5fc7508c 2277 if (!track)
2278 continue;
2279 if (track->Pt()<pt)
2280 continue;
2281 Float_t eta = track->Eta();
2282 Float_t phi = track->Phi();
2283 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
2284 Double_t etadiff = (eta-cEta);
2285 if(TMath::Abs(etadiff)>dEta || TMath::Abs(phidiff)>dPhi)
2286 continue;
2287 trkIsolation += track->Pt();
2288 }
2289 return trkIsolation;
2290}
2291
0fbe8d4f 2292//________________________________________________________________________
2293Bool_t AliAnalysisTaskEMCALPi0PbPb::IsShared(const AliVCluster *c) const
2294{
2295 // Returns if cluster shared across super modules.
2296
1f41ed3d 2297 if (!c)
0fbe8d4f 2298 return 0;
2299
2300 Int_t n = -1;
2301 Int_t ncells = c->GetNCells();
2302 for(Int_t j=0; j<ncells; ++j) {
2303 Int_t id = TMath::Abs(c->GetCellAbsId(j));
2304 Int_t got = id / (24*48);
2305 if (n==-1) {
2306 n = got;
2307 continue;
2308 }
2309 if (got!=n)
2310 return 1;
2311 }
2312 return 0;
2313}
3a952328 2314
2ee7bda4 2315//________________________________________________________________________
2316Bool_t AliAnalysisTaskEMCALPi0PbPb::IsIdPartOfCluster(const AliVCluster *c, Short_t id) const
2317{
2318 // Returns if id is part of cluster.
2319
2320 AliVCaloCells *cells = fEsdCells;
2321 if (!cells)
2322 cells = fAodCells;
2323 if (!cells)
2324 return 0;
2325
2326 Int_t ncells = c->GetNCells();
2327 for(Int_t j=0; j<ncells; ++j) {
2328 Int_t cid = TMath::Abs(c->GetCellAbsId(j));
2329 if (cid == id)
2330 return 1;
2331 }
2332 return 0;
2333}
2334
38727e64 2335//________________________________________________________________________
807016ea 2336void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(const AliVParticle *p, const TObjArray *arr, Int_t level) const
38727e64 2337{
56fd6cb2 2338 // Print recursively daughter information.
2339
2340 if (!p || !arr)
2341 return;
2342
807016ea 2343 const AliAODMCParticle *amc = dynamic_cast<const AliAODMCParticle*>(p);
2344 if (!amc)
2345 return;
2346 for (Int_t i=0; i<level; ++i) printf(" ");
2347 amc->Print();
2348
2349 Int_t n = amc->GetNDaughters();
2350 for (Int_t i=0; i<n; ++i) {
2351 Int_t d = amc->GetDaughter(i);
2352 const AliVParticle *dmc = static_cast<const AliVParticle*>(arr->At(d));
2353 PrintDaughters(dmc,arr,level+1);
2354 }
2355}
2356
2357//________________________________________________________________________
2358void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(const AliMCParticle *p, const AliMCEvent *arr, Int_t level) const
2359{
2360 // Print recursively daughter information.
2361
2362 if (!p || !arr)
2363 return;
2364
2365 for (Int_t i=0; i<level; ++i) printf(" ");
2366 Int_t d1 = p->GetFirstDaughter();
2367 Int_t d2 = p->GetLastDaughter();
2368 printf("pid=%d: %.2f %.2f %.2f (%.2f %.2f %.2f); nd=%d,%d\n",
2369 p->PdgCode(),p->Px(),p->Py(),p->Pz(),p->Xv(),p->Yv(),p->Zv(),d1,d2);
2370 if (d1<0)
2371 return;
2372 if (d2<0)
2373 d2=d1;
2374 for (Int_t i=d1;i<=d2;++i) {
2375 const AliMCParticle *dmc = static_cast<const AliMCParticle *>(arr->GetTrack(i));
6d4faab0 2376 PrintDaughters(dmc,arr,level+1);
56fd6cb2 2377 }
38727e64 2378}
2379
2380//________________________________________________________________________
2381void AliAnalysisTaskEMCALPi0PbPb::PrintTrackRefs(AliMCParticle *p) const
2382{
56fd6cb2 2383 // Print track ref array.
2384
2385 if (!p)
2386 return;
2387
2388 Int_t n = p->GetNumberOfTrackReferences();
2389 for (Int_t i=0; i<n; ++i) {
2390 AliTrackReference *ref = p->GetTrackReference(i);
bc8a45bd 2391 if (!ref)
2392 continue;
56fd6cb2 2393 ref->SetUserId(ref->DetectorId());
2394 ref->Print();
2395 }
38727e64 2396}
2397
807016ea 2398//________________________________________________________________________
2399void AliAnalysisTaskEMCALPi0PbPb::ProcessDaughters(AliVParticle *p, Int_t index, const TObjArray *arr)
2400{
2401 // Process and create daughters.
2402
2403 if (!p || !arr)
2404 return;
2405
2406 AliAODMCParticle *amc = dynamic_cast<AliAODMCParticle*>(p);
2407 if (!amc)
2408 return;
2409
2410 //amc->Print();
2411
2412 Int_t nparts = arr->GetEntries();
2413 Int_t nents = fMcParts->GetEntries();
2414
2415 AliStaPart *newp = static_cast<AliStaPart*>(fMcParts->New(nents));
2416 newp->fPt = amc->Pt();
2417 newp->fEta = amc->Eta();
2418 newp->fPhi = amc->Phi();
2419 if (amc->Xv() != 0 || amc->Yv() != 0 || amc->Zv() != 0) {
2420 TVector3 vec(amc->Xv(),amc->Yv(),amc->Zv());
2421 newp->fVR = vec.Perp();
2422 newp->fVEta = vec.Eta();
2423 newp->fVPhi = vec.Phi();
2424 }
2425 newp->fPid = amc->PdgCode();
8c56d760 2426 newp->fLab = nents;
807016ea 2427 Int_t moi = amc->GetMother();
2428 if (moi>=0&&moi<nparts) {
2429 const AliAODMCParticle *mmc = static_cast<const AliAODMCParticle*>(arr->At(moi));
2430 moi = mmc->GetUniqueID();
2431 }
2432 newp->fMo = moi;
2433 p->SetUniqueID(nents);
8c56d760 2434
2435 // TODO: Determine which detector was hit
2436 //newp->fDet = ???
2437
807016ea 2438 Int_t n = amc->GetNDaughters();
2439 for (Int_t i=0; i<n; ++i) {
2440 Int_t d = amc->GetDaughter(i);
2441 if (d<=index || d>=nparts)
2442 continue;
2443 AliVParticle *dmc = static_cast<AliVParticle*>(arr->At(d));
2444 ProcessDaughters(dmc,d,arr);
2445 }
2446}
2447
2448//________________________________________________________________________
2449void AliAnalysisTaskEMCALPi0PbPb::ProcessDaughters(AliMCParticle *p, Int_t index, const AliMCEvent *arr)
2450{
2451 // Process and create daughters.
2452
2453 if (!p || !arr)
2454 return;
2455
2456 Int_t d1 = p->GetFirstDaughter();
2457 Int_t d2 = p->GetLastDaughter();
2458 if (0) {
2459 printf("%d pid=%d: %.3f %.3f %.3f (%.2f %.2f %.2f); nd=%d,%d, mo=%d\n",
2460 index,p->PdgCode(),p->Px(),p->Py(),p->Pz(),p->Xv(),p->Yv(),p->Zv(),d1,d2, p->GetMother());
bc8a45bd 2461 PrintTrackRefs(p);
807016ea 2462 }
2463 Int_t nents = fMcParts->GetEntries();
2464
2465 AliStaPart *newp = static_cast<AliStaPart*>(fMcParts->New(nents));
2466 newp->fPt = p->Pt();
2467 newp->fEta = p->Eta();
2468 newp->fPhi = p->Phi();
2469 if (p->Xv() != 0 || p->Yv() != 0 || p->Zv() != 0) {
2470 TVector3 vec(p->Xv(),p->Yv(),p->Zv());
2471 newp->fVR = vec.Perp();
2472 newp->fVEta = vec.Eta();
2473 newp->fVPhi = vec.Phi();
2474 }
2475 newp->fPid = p->PdgCode();
8c56d760 2476 newp->fLab = nents;
807016ea 2477 Int_t moi = p->GetMother();
2478 if (moi>=0) {
2479 const AliMCParticle *mmc = static_cast<const AliMCParticle *>(arr->GetTrack(moi));
2480 moi = mmc->GetUniqueID();
2481 }
2482 newp->fMo = moi;
2483 p->SetUniqueID(nents);
2484
2485 Int_t nref = p->GetNumberOfTrackReferences();
2486 if (nref>0) {
2487 AliTrackReference *ref = p->GetTrackReference(nref-1);
2488 if (ref) {
2489 newp->fDet = ref->DetectorId();
2490 }
2491 }
2492
2493 if (d1<0)
2494 return;
2495 if (d2<0)
2496 d2=d1;
2497 for (Int_t i=d1;i<=d2;++i) {
2498 AliMCParticle *dmc = static_cast<AliMCParticle *>(arr->GetTrack(i));
2499 if (dmc->P()<0.01)
2500 continue;
2501 ProcessDaughters(dmc,i,arr);
2502 }
2503}
6dbf6b27 2504
2505//__________________________________________________________________________________________________
2506void AliStaCluster::GetMom(TLorentzVector& p, Double_t *vertex)
2507{
2508 // Calculate momentum.
2509
2510 TVector3 pos;
2511 pos.SetPtEtaPhi(fR,fEta,fPhi);
2512
a186e444 2513 if(vertex){ //calculate direction relative to vertex
6dbf6b27 2514 pos -= vertex;
2515 }
2516
2517 Double_t r = pos.Mag();
2518 p.SetPxPyPzE(fE*pos.x()/r, fE*pos.y()/r, fE*pos.z()/r, fE);
2519}
2520
2521//__________________________________________________________________________________________________
2522void AliStaCluster::GetMom(TLorentzVector& p, AliStaVertex *vertex)
2523{
2524 // Calculate momentum.
2525
2526 Double_t v[3] = {0,0,0};
2527 if (vertex) {
2528 v[0] = vertex->fVx;
2529 v[1] = vertex->fVy;
2530 v[2] = vertex->fVz;
2531 }
2532 GetMom(p, v);
2533}