]> 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
Merge branch 'feature-movesplit'
[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");
0a918d8d 606 offtrigger = ((AliVAODHeader*)fAodEv->GetHeader())->GetOfflineTrigger();
44310bce 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);
0a918d8d 629 else {
630 AliAODHeader * aodheader = dynamic_cast<AliAODHeader*>(fAodEv->GetHeader());
631 if(!aodheader) AliFatal("Not a standard AOD");
632 geom = aodheader->GetEMCALMatrix(i);
633 }
9809ed8c 634 if (!geom)
635 continue;
636 geom->Print();
637 fGeom->SetMisalMatrix(geom,i);
27c2e3d9 638 }
639 fIsGeoMatsSet = kTRUE;
640 }
641
642 if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
643 if (fEsdEv) {
644 const AliESDRun *erun = fEsdEv->GetESDRun();
645 AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(),
646 erun->GetCurrentDip(),
647 AliMagF::kConvLHC,
648 kFALSE,
649 erun->GetBeamEnergy(),
650 erun->GetBeamType());
651 TGeoGlobalMagField::Instance()->SetField(field);
652 } else {
653 Double_t pol = -1; //polarity
654 Double_t be = -1; //beam energy
655 AliMagF::BeamType_t btype = AliMagF::kBeamTypepp;
656 Int_t runno = fAodEv->GetRunNumber();
657 if (runno>=136851 && runno<138275) {
658 pol = -1;
659 be = 2760;
660 btype = AliMagF::kBeamTypeAA;
661 } else if (runno>=138275 && runno<=139517) {
662 pol = +1;
663 be = 2760;
664 btype = AliMagF::kBeamTypeAA;
665 } else {
666 AliError(Form("Do not know the bfield parameters for run %d! Using defaults!!!", runno));
667 }
668 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", pol, pol, AliMagF::k5kG, btype, be));
669 }
670 }
671
b3ee6797 672 Int_t cut = 1;
673 fHCuts->Fill(cut++);
674
675 TString trgclasses;
676 AliESDHeader *h = dynamic_cast<AliESDHeader*>(InputEvent()->GetHeader());
677 if (h) {
678 trgclasses = fEsdEv->GetFiredTriggerClasses();
679 } else {
680 AliAODHeader *h2 = dynamic_cast<AliAODHeader*>(InputEvent()->GetHeader());
681 if (h2)
682 trgclasses = h2->GetFiredTriggerClasses();
683 }
684 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
685 const char *name = fTrClassNamesArr->At(i)->GetName();
0ceca1f8 686 TRegexp regexp(name);
687 if (trgclasses.Contains(regexp))
b3ee6797 688 fHTclsBeforeCuts->Fill(1+i);
689 }
690
807016ea 691 if (fDoPSel && offtrigger==0)
b3ee6797 692 return;
693
fa443410 694 fHCuts->Fill(cut++);
286b47a5 695
717fe7de 696 const AliCentrality *centP = InputEvent()->GetCentrality();
697 Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar);
fa443410 698 fHCent->Fill(cent);
717fe7de 699 if (cent<fCentFrom||cent>fCentTo)
286b47a5 700 return;
43cfaa06 701
fa443410 702 fHCuts->Fill(cut++);
286b47a5 703
43cfaa06 704 if (fUseQualFlag) {
705 if (centP->GetQuality()>0)
706 return;
707 }
286b47a5 708
fa443410 709 fHCentQual->Fill(cent);
710 fHCuts->Fill(cut++);
717fe7de 711
43cfaa06 712 if (fEsdEv) {
fa443410 713 am->LoadBranch("PrimaryVertex.");
714 am->LoadBranch("SPDVertex.");
715 am->LoadBranch("TPCVertex.");
43cfaa06 716 } else {
717 fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
718 am->LoadBranch("vertices");
3c661da5 719 if (!fAodEv) return;
43cfaa06 720 }
721
722 const AliVVertex *vertex = InputEvent()->GetPrimaryVertex();
717fe7de 723 if (!vertex)
724 return;
725
fa443410 726 fHVertexZ->Fill(vertex->GetZ());
286b47a5 727
717fe7de 728 if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax)
729 return;
286b47a5 730
fa443410 731 fHCuts->Fill(cut++);
76332037 732 fHVertexZ2->Fill(vertex->GetZ());
717fe7de 733
b3ee6797 734 // count number of accepted events
735 ++fNEvs;
736
737 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
738 const char *name = fTrClassNamesArr->At(i)->GetName();
0ceca1f8 739 TRegexp regexp(name);
740 if (trgclasses.Contains(regexp))
b3ee6797 741 fHTclsAfterCuts->Fill(1+i);
742 }
743
43cfaa06 744 fRecPoints = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
c2fe5f0e 745 fDigits = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
b6c599fe 746 fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set or if clusters are attached
43cfaa06 747 fEsdCells = 0; // will be set if ESD input used
b6c599fe 748 fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set or if clusters are attached
43cfaa06 749 // or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode
750 fAodCells = 0; // will be set if AOD input used
717fe7de 751
43cfaa06 752 // deal with special output from AliAnalysisTaskEMCALClusterizeFast first
2ee7bda4 753 Bool_t overwrite = 0;
43cfaa06 754 Bool_t clusattached = 0;
755 Bool_t recalibrated = 0;
756 if (1 && !fClusName.IsNull()) {
757 AliAnalysisTaskEMCALClusterizeFast *cltask = 0;
758 TObjArray *ts = am->GetTasks();
759 cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName));
760 if (cltask && cltask->GetClusters()) {
2ee7bda4 761 fRecPoints = cltask->GetClusters();
762 fDigits = cltask->GetDigits();
43cfaa06 763 clusattached = cltask->GetAttachClusters();
2ee7bda4 764 overwrite = cltask->GetOverwrite();
43cfaa06 765 if (cltask->GetCalibData()!=0)
766 recalibrated = kTRUE;
767 }
768 }
38727e64 769 if (1 && !fClusName.IsNull()) {
770 TList *l = 0;
771 if (AODEvent())
772 l = AODEvent()->GetList();
773 else if (fAodEv)
774 l = fAodEv->GetList();
717fe7de 775 if (l) {
38727e64 776 fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
717fe7de 777 }
ea3fd2d5 778 }
717fe7de 779
780 if (fEsdEv) { // ESD input mode
43cfaa06 781 if (1 && (!fRecPoints||clusattached)) {
2ee7bda4 782 if (!clusattached && !overwrite)
43cfaa06 783 am->LoadBranch("CaloClusters");
717fe7de 784 TList *l = fEsdEv->GetList();
2ee7bda4 785 if (clusattached) {
786 fEsdClusters = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
1829ebfe 787 } else {
38727e64 788 fEsdClusters = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
ea3fd2d5 789 }
790 }
43cfaa06 791 if (1) {
792 if (!recalibrated)
793 am->LoadBranch("EMCALCells.");
717fe7de 794 fEsdCells = fEsdEv->GetEMCALCells();
795 }
796 } else if (fAodEv) { // AOD input mode
43cfaa06 797 if (1 && (!fAodClusters || clusattached)) {
798 if (!clusattached)
799 am->LoadBranch("caloClusters");
717fe7de 800 TList *l = fAodEv->GetList();
717fe7de 801 if (l) {
38727e64 802 fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
ea3fd2d5 803 }
804 }
43cfaa06 805 if (1) {
806 if (!recalibrated)
807 am->LoadBranch("emcalCells");
717fe7de 808 fAodCells = fAodEv->GetEMCALCells();
809 }
810 } else {
811 AliFatal("Impossible to not have either pointer to ESD or AOD event");
812 }
ea3fd2d5 813
43cfaa06 814 if (1) {
815 AliDebug(2,Form("fRecPoints set: %p", fRecPoints));
c2fe5f0e 816 AliDebug(2,Form("fDigits set: %p", fDigits));
43cfaa06 817 AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters));
818 AliDebug(2,Form("fEsdCells set: %p", fEsdCells));
819 AliDebug(2,Form("fAodClusters set: %p", fAodClusters));
820 AliDebug(2,Form("fAodCells set: %p", fAodCells));
821 }
822
a49742b5 823 if (fDoAfterburner)
824 ClusterAfterburner();
6eb6260e 825
38727e64 826 CalcMcInfo();
3a952328 827 CalcCaloTriggers();
b6c599fe 828 CalcPrimTracks();
3a952328 829 CalcTracks();
296ea9b4 830 CalcClusterProps();
831
286b47a5 832 FillCellHists();
3a952328 833 if (!fTrainMode) {
834 FillClusHists();
835 FillPionHists();
836 FillOtherHists();
837 }
38727e64 838 FillMcHists();
788ca675 839 FillNtuple();
ea3fd2d5 840
3a952328 841 if (fTrainMode) {
842 fSelTracks->Clear();
843 fSelPrimTracks->Clear();
cfd7d5b2 844 if (fMcParts)
807016ea 845 fMcParts->Clear();
4ea96211 846 if (fTriggers)
847 fTriggers->Clear();
848 if (fClusters)
849 fClusters->Clear();
3a952328 850 }
de4cee41 851
717fe7de 852 PostData(1, fOutput);
ea3fd2d5 853}
854
855//________________________________________________________________________
856void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *)
857{
717fe7de 858 // Terminate called at the end of analysis.
f5d4ab70 859
860 if (fNtuple) {
4ea96211 861 TFile *f = OpenFile(1);
a9c1389d 862 TDirectory::TContext context(f);
f5d4ab70 863 if (f)
864 fNtuple->Write();
865 }
d9f26424 866
bc8a45bd 867 AliInfo(Form("%s: Accepted %lld events ", GetName(), fNEvs));
ea3fd2d5 868}
ea3fd2d5 869
296ea9b4 870//________________________________________________________________________
3a952328 871void AliAnalysisTaskEMCALPi0PbPb::CalcCaloTriggers()
296ea9b4 872{
3a952328 873 // Calculate triggers
296ea9b4 874
38727e64 875 if (fAodEv)
876 return; // information not available in AOD
877
4ea96211 878 if (!fTriggers)
879 return;
880
3a952328 881 fTriggers->Clear();
296ea9b4 882
2ee7bda4 883 if (fTrigName.Length()<=0)
3a952328 884 return;
885
2ee7bda4 886 TClonesArray *arr = dynamic_cast<TClonesArray*>(fEsdEv->FindListObject(fTrigName));
887 if (!arr) {
888 AliError(Form("Could not get array with name %s", fTrigName.Data()));
3a952328 889 return;
3a952328 890 }
891
2ee7bda4 892 Int_t nNumberOfCaloClusters = arr->GetEntries();
893 for(Int_t j = 0, ntrigs = 0; j < nNumberOfCaloClusters; ++j) {
894 AliVCluster *cl = dynamic_cast<AliVCluster*>(arr->At(j));
895 if (!cl)
3a952328 896 continue;
2ee7bda4 897 if (!cl->IsEMCAL())
3a952328 898 continue;
2ee7bda4 899 if (cl->E()<1)
c47674e0 900 continue;
3a952328 901 AliStaTrigger *trignew = static_cast<AliStaTrigger*>(fTriggers->New(ntrigs++));
2ee7bda4 902 Float_t pos[3] = {0,0,0};
903 cl->GetPosition(pos);
904 TVector3 vpos(pos);
905 trignew->fE = cl->E();
906 trignew->fEta = vpos.Eta();
907 trignew->fPhi = vpos.Phi();
908 Short_t id = -1;
909 GetMaxCellEnergy(cl, id);
910 trignew->fIdMax = id;
3a952328 911 }
912}
913
914//________________________________________________________________________
915void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
916{
917 // Calculate cluster properties
918
4ea96211 919 if (!fClusters)
920 return;
921
3a952328 922 fClusters->Clear();
923
924 AliVCaloCells *cells = fEsdCells;
925 if (!cells)
926 cells = fAodCells;
927 if (!cells)
928 return;
929
930 TObjArray *clusters = fEsdClusters;
931 if (!clusters)
932 clusters = fAodClusters;
933 if (!clusters)
934 return;
935
a186e444 936 const Int_t ncells = cells->GetNumberOfCells();
937 const Int_t nclus = clusters->GetEntries();
938 const Int_t ntrks = fSelTracks->GetEntries();
3a952328 939
940 std::map<Short_t,Short_t> map;
a186e444 941 for (Short_t pos=0; pos<ncells; ++pos) {
3a952328 942 Short_t id = cells->GetCellNumber(pos);
a186e444 943 if (id>24*48*fGeom->GetNumberOfSuperModules()) {
944 AliError(Form("Id of cell found to be too large: %d", id));
945 continue;
946 }
947 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigits->At(pos));
948 if (digit->GetId() != id) {
949 AliError(Form("Ids should be equal: %d %d", id, digit->GetId()));
950 return;
951 }
952 map[id] = pos;
3a952328 953 }
954
8c56d760 955 TObjArray filtMcParts;
cfd7d5b2 956 if (fMcParts) {
8c56d760 957 Int_t nmc = fMcParts->GetEntries();
958 for (Int_t i=0; i<nmc; ++i) {
959 AliStaPart *pa = static_cast<AliStaPart*>(fMcParts->At(i));
960 if (pa->OnEmcal())
961 filtMcParts.Add(pa);
962 }
963 }
964
fe38566a 965 Double_t vertex[3] = {0,0,0};
966 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
967
3a952328 968 for(Int_t i=0, ncl=0; i<nclus; ++i) {
969 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
b5ba4932 970
3a952328 971 if (!clus)
972 continue;
973 if (!clus->IsEMCAL())
974 continue;
975 if (clus->E()<fMinE)
976 continue;
977
fe38566a 978 Float_t clsPos[3] = {0,0,0};
3a952328 979 clus->GetPosition(clsPos);
980 TVector3 clsVec(clsPos);
3a952328 981 TLorentzVector clusterVec;
982 clus->GetMomentum(clusterVec,vertex);
983 Double_t clsEta = clusterVec.Eta();
984
985 AliStaCluster *cl = static_cast<AliStaCluster*>(fClusters->New(ncl++));
986 cl->fE = clus->E();
987 cl->fR = clsVec.Perp();
988 cl->fEta = clsVec.Eta();
989 cl->fPhi = clsVec.Phi();
990 cl->fN = clus->GetNCells();
991 cl->fN1 = GetNCells(clus,0.100);
992 cl->fN3 = GetNCells(clus,0.300);
993 Short_t id = -1;
994 Double_t emax = GetMaxCellEnergy(clus, id);
995 cl->fIdMax = id;
2ee7bda4 996 cl->fSM = fGeom->GetSuperModuleNumber(id);
3a952328 997 cl->fEmax = emax;
1f41ed3d 998 Short_t id2 = -1;
999 cl->fE2max = GetSecondMaxCellEnergy(clus,id2);
2ee7bda4 1000 cl->fTmax = cells->GetCellTime(id);
3a952328 1001 if (clus->GetDistanceToBadChannel()<10000)
1002 cl->fDbc = clus->GetDistanceToBadChannel();
1003 if (!TMath::IsNaN(clus->GetDispersion()))
1004 cl->fDisp = clus->GetDispersion();
1005 if (!TMath::IsNaN(clus->GetM20()))
0fbe8d4f 1006 cl->fM20 = clus->GetM20();
3a952328 1007 if (!TMath::IsNaN(clus->GetM02()))
0fbe8d4f 1008 cl->fM02 = clus->GetM02();
2ee7bda4 1009 Double_t maxAxis = -1, minAxis = -1;
3a952328 1010 GetSigma(clus,maxAxis,minAxis);
2ee7bda4 1011 clus->SetTOF(maxAxis); // store sigma in TOF for later plotting
3a952328 1012 cl->fSig = maxAxis;
2ee7bda4 1013 Double_t sEtaEta = -1;
1014 Double_t sPhiPhi = -1;
5fc7508c 1015 GetSigmaEtaEta(clus, sEtaEta, sPhiPhi);
1016 cl->fSigEtaEta = sEtaEta;
1017 cl->fSigPhiPhi = sPhiPhi;
2ee7bda4 1018 Double_t clusterEcc = -1;
3a952328 1019 if (maxAxis > 0)
1020 clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
2ee7bda4 1021 clus->SetChi2(clusterEcc); // store ecc in chi2 for later plotting
1022 cl->fEcc = clusterEcc;
1023 cl->fTrIso = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
1024 cl->fTrIso1 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 1);
1025 cl->fTrIso2 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 2);
5fc7508c 1026 cl->fTrIsoD1 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist-0.1);
1027 cl->fTrIso1D1 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist-0.1, 1);
1028 cl->fTrIso2D1 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist-0.1, 2);
1029 cl->fTrIsoD3 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.1);
1030 cl->fTrIso1D3 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.1, 1);
1031 cl->fTrIso2D3 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.1, 2);
1f41ed3d 1032 cl->fTrIsoD4 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.2);
1033 cl->fTrIso1D4 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.2, 1);
1034 cl->fTrIso2D4 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.2, 2);
5fc7508c 1035 cl->fTrIsoStrip = GetTrackIsoStrip(clusterVec.Eta(), clusterVec.Phi());
2ee7bda4 1036 cl->fCeCore = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.05);
1037 cl->fCeIso = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
1038 cl->fCeIso1 = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.10);
1039 cl->fCeIso3 = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.30);
1f41ed3d 1040 cl->fCeIso4 = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.40);
a186e444 1041 cl->fCeIso3x3 = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 3, 3);
2ee7bda4 1042 cl->fCeIso4x4 = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 4, 4);
1043 cl->fCeIso5x5 = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 5, 5);
1044 cl->fCeIso3x22 = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 3, 22);
1045 cl->fIsShared = IsShared(clus);
1046 cl->fTrigId = -1;
1047 cl->fTrigE = 0;
e0e1022c 1048
2ee7bda4 1049 if (fTriggers) {
1050 Int_t ntrig = fTriggers->GetEntries();
1051 for (Int_t j = 0; j<ntrig; ++j) {
1052 AliStaTrigger *sta = static_cast<AliStaTrigger*>(fTriggers->At(j));
1053 if (!sta)
38727e64 1054 continue;
2ee7bda4 1055 Short_t idmax = sta->fIdMax;
1056 Bool_t inc = IsIdPartOfCluster(clus, idmax);
1057 if (inc) {
1058 cl->fTrigId = j;
1059 cl->fTrigE = sta->fE;
1060 break;
1061 }
38727e64 1062 }
3a952328 1063 }
1064
1065 // track matching
a186e444 1066 if (fDoTrMatGeom) {
1067 Double_t mind2 = 1e10;
1068 for(Int_t j = 0; j<ntrks; ++j) {
1069 AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
1070 if (!track)
1071 continue;
1072 if (TMath::Abs(clsEta-track->Eta())>0.5)
1073 continue;
1074 TVector3 vec(clsPos);
1075 Float_t tmpEta=-1, tmpPhi=-1, tmpR2=1e10;
1076 Double_t dedx = -1;
1077 if (fEsdEv) {
1078 AliESDtrack *esdTrack = static_cast<AliESDtrack*>(track);
1079 dedx = esdTrack->GetTPCsignal();
1080 }
3a952328 1081 if (TMath::Abs(clsEta-track->Eta())>fIsoDist)
1082 continue;
e0e1022c 1083 AliExternalTrackParam tParam;
1084 tParam.CopyFromVTrack(track);
1085 if (!fReco->ExtrapolateTrackToCluster(&tParam, clus, tmpEta, tmpPhi))
3a952328 1086 continue;
e0e1022c 1087 tmpR2 = tmpEta*tmpEta + tmpPhi*tmpPhi;
a186e444 1088 Double_t d2 = tmpR2;
1089 if (mind2>d2) {
1090 mind2=d2;
1091 cl->fTrDz = tmpEta;
1092 cl->fTrDr = tmpPhi;
1093 cl->fTrEp = clus->E()/track->P();
1094 cl->fTrDedx = dedx;
1095 cl->fIsTrackM = 1;
1096 }
b6c599fe 1097 }
a186e444 1098 if (cl->fIsTrackM) {
1099 fHMatchDr->Fill(cl->fTrDr);
1100 fHMatchDz->Fill(cl->fTrDz);
1101 fHMatchEp->Fill(cl->fTrEp);
3a952328 1102 }
1103 }
8c56d760 1104
1105 //mc matching
cfd7d5b2 1106 if (fMcParts) {
8c56d760 1107 Int_t nmc = filtMcParts.GetEntries();
1108 Double_t diffR2 = 1e9;
1109 AliStaPart *msta = 0;
1110 for (Int_t j=0; j<nmc; ++j) {
1111 AliStaPart *pa = static_cast<AliStaPart*>(filtMcParts.At(j));
1112 Double_t t1=clsVec.Eta()-pa->fVEta;
1113 Double_t t2=TVector2::Phi_mpi_pi(clsVec.Phi()-pa->fVPhi);
1114 Double_t tmp = t1*t1+t2*t2;
1115 if (tmp<diffR2) {
1116 diffR2 = tmp;
1117 msta = pa;
1118 }
1119 }
c2fe5f0e 1120 if (diffR2<10 && msta!=0) {
1121 cl->fMcLabel = msta->fLab;
1122 }
1123 }
1124
1125 cl->fEmbE = 0;
1126 if (fDigits && fEmbedMode) {
1127 for(Int_t j=0; j<cl->fN; ++j) {
1128 Short_t cid = TMath::Abs(clus->GetCellAbsId(j));
1129 Short_t pos = -1;
1130 std::map<Short_t,Short_t>::iterator it = map.find(cid);
1131 if (it!=map.end())
1132 pos = it->second;
1133 if (pos<0)
1134 continue;
1135 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigits->At(pos));
1136 if (!digit)
1137 continue;
c2fe5f0e 1138 if (digit->GetType()<-1) {
1139 cl->fEmbE += digit->GetChi2();
1140 }
1141 }
8c56d760 1142 }
b6c599fe 1143 }
1144}
1145
1146//________________________________________________________________________
1147void AliAnalysisTaskEMCALPi0PbPb::CalcPrimTracks()
1148{
a2de5ca1 1149 // Calculate track properties for primary tracks.
b6c599fe 1150
1151 fSelPrimTracks->Clear();
1152
1153 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1154 TClonesArray *tracks = 0;
e0e1022c 1155 Bool_t doConstrain = kFALSE;
1156 TString trkname(fPrimTracksName);
b6c599fe 1157 if (fEsdEv) {
e0e1022c 1158 if (trkname.IsNull()) {
1159 trkname = "Tracks";
1160 am->LoadBranch("Tracks");
1161 fSelPrimTracks->SetOwner(kTRUE);
1162 doConstrain = kTRUE;
1163 }
b6c599fe 1164 TList *l = fEsdEv->GetList();
e0e1022c 1165 tracks = dynamic_cast<TClonesArray*>(l->FindObject(trkname));
b6c599fe 1166 } else {
e0e1022c 1167 trkname = "tracks";
b6c599fe 1168 am->LoadBranch("tracks");
1169 TList *l = fAodEv->GetList();
1170 tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1171 }
1172
e0e1022c 1173 if (!tracks) {
1174 AliError(Form("Pointer to tracks %s == 0", trkname.Data() ));
296ea9b4 1175 return;
e0e1022c 1176 }
296ea9b4 1177
0ec74551 1178 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
b6c599fe 1179 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-fIsoDist*1.25;
1180 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+fIsoDist*1.25;
0ec74551 1181
296ea9b4 1182 if (fEsdEv) {
0ec74551 1183 am->LoadBranch("PrimaryVertex.");
1184 const AliESDVertex *vtx = fEsdEv->GetPrimaryVertexSPD();
1185 am->LoadBranch("SPDVertex.");
1186 const AliESDVertex *vtxSPD = fEsdEv->GetPrimaryVertexSPD();
1187 am->LoadBranch("Tracks");
1188 const Int_t Ntracks = tracks->GetEntries();
1189 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1190 AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1191 if (!track) {
1192 AliWarning(Form("Could not receive track %d\n", iTracks));
1193 continue;
1194 }
e0e1022c 1195 if (fPrimTrCuts && !fPrimTrCuts->IsSelected(track))
296ea9b4 1196 continue;
1197 Double_t eta = track->Eta();
1198 if (eta<-1||eta>1)
1199 continue;
0ec74551 1200 if (track->Phi()<phimin||track->Phi()>phimax)
1201 continue;
e0e1022c 1202 if (!doConstrain) {
1203 fSelPrimTracks->Add(track);
1204 continue;
1205 }
0ec74551 1206 AliESDtrack copyt(*track);
1207 Double_t bfield[3];
1208 copyt.GetBxByBz(bfield);
1209 AliExternalTrackParam tParam;
1210 Bool_t relate = copyt.RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&tParam);
1211 if (!relate)
1212 continue;
1213 copyt.Set(tParam.GetX(),tParam.GetAlpha(),tParam.GetParameter(),tParam.GetCovariance());
1214 Double_t p[3] = { 0. };
1215 copyt.GetPxPyPz(p);
1216 Double_t pos[3] = { 0. };
1217 copyt.GetXYZ(pos);
1218 Double_t covTr[21] = { 0. };
1219 copyt.GetCovarianceXYZPxPyPz(covTr);
76b98553 1220 // Double_t pid[10] = { 0. };
1221 // copyt.GetESDpid(pid);
0ec74551 1222 AliAODTrack *aTrack = new AliAODTrack(copyt.GetID(),
1223 copyt.GetLabel(),
1224 p,
1225 kTRUE,
1226 pos,
1227 kFALSE,
1228 covTr,
1229 (Short_t)copyt.GetSign(),
1230 copyt.GetITSClusterMap(),
76b98553 1231 //pid,
0ec74551 1232 0,/*fPrimaryVertex,*/
1233 kTRUE, // check if this is right
1234 vtx->UsesTrack(copyt.GetID()));
76b98553 1235 aTrack->SetPIDForTracking(copyt.GetPIDForTracking());
0ec74551 1236 aTrack->SetTPCClusterMap(copyt.GetTPCClusterMap());
1237 aTrack->SetTPCSharedMap (copyt.GetTPCSharedMap());
1238 Float_t ndf = copyt.GetTPCNcls() + 1 - 5 ;
1239 if(ndf>0)
1240 aTrack->SetChi2perNDF(copyt.GetTPCchi2()/ndf);
1241 else
1242 aTrack->SetChi2perNDF(-1);
1243 aTrack->SetFlags(copyt.GetStatus());
1244 aTrack->SetTPCPointsF(copyt.GetTPCNclsF());
b6c599fe 1245 fSelPrimTracks->Add(aTrack);
296ea9b4 1246 }
1247 } else {
1248 Int_t ntracks = tracks->GetEntries();
1249 for (Int_t i=0; i<ntracks; ++i) {
1250 AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
1251 if (!track)
1252 continue;
1253 Double_t eta = track->Eta();
1254 if (eta<-1||eta>1)
1255 continue;
0ec74551 1256 if (track->Phi()<phimin||track->Phi()>phimax)
1257 continue;
b6c599fe 1258 if(track->GetTPCNcls()<fMinNClusPerTr)
296ea9b4 1259 continue;
b6c599fe 1260 //todo: Learn how to set/filter AODs for prim/sec tracks
1261 fSelPrimTracks->Add(track);
296ea9b4 1262 }
1263 }
1264}
1265
1266//________________________________________________________________________
3a952328 1267void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
296ea9b4 1268{
3a952328 1269 // Calculate track properties (including secondaries).
b6c599fe 1270
3a952328 1271 fSelTracks->Clear();
296ea9b4 1272
3a952328 1273 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1274 TClonesArray *tracks = 0;
1275 if (fEsdEv) {
1276 am->LoadBranch("Tracks");
1277 TList *l = fEsdEv->GetList();
1278 tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
1279 } else {
1280 am->LoadBranch("tracks");
1281 TList *l = fAodEv->GetList();
1282 tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1283 }
296ea9b4 1284
3a952328 1285 if (!tracks)
1286 return;
296ea9b4 1287
3a952328 1288 if (fEsdEv) {
1289 const Int_t Ntracks = tracks->GetEntries();
1290 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1291 AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1292 if (!track) {
1293 AliWarning(Form("Could not receive track %d\n", iTracks));
1294 continue;
1295 }
1296 if (fTrCuts && !fTrCuts->IsSelected(track))
1297 continue;
1298 Double_t eta = track->Eta();
1299 if (eta<-1||eta>1)
1300 continue;
1301 fSelTracks->Add(track);
1302 }
1303 } else {
1304 Int_t ntracks = tracks->GetEntries();
1305 for (Int_t i=0; i<ntracks; ++i) {
1306 AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
296ea9b4 1307 if (!track)
1308 continue;
3a952328 1309 Double_t eta = track->Eta();
1310 if (eta<-1||eta>1)
c4236a58 1311 continue;
3a952328 1312 if(track->GetTPCNcls()<fMinNClusPerTr)
b6c599fe 1313 continue;
2e4d8148 1314
3a952328 1315 if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL
e0e1022c 1316 AliExternalTrackParam tParam;
1317 tParam.CopyFromVTrack(track);
a186e444 1318 if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 430, 0.139, 1, kTRUE)) {
3a952328 1319 Double_t trkPos[3];
1320 tParam.GetXYZ(trkPos);
1321 track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]);
1322 }
b6c599fe 1323 }
3a952328 1324 fSelTracks->Add(track);
c4236a58 1325 }
296ea9b4 1326 }
1327}
1328
323834f0 1329//________________________________________________________________________
3a952328 1330void AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner()
323834f0 1331{
296ea9b4 1332 // Run custer reconstruction afterburner.
323834f0 1333
1334 AliVCaloCells *cells = fEsdCells;
1335 if (!cells)
1336 cells = fAodCells;
1337
1338 if (!cells)
1339 return;
1340
1341 Int_t ncells = cells->GetNumberOfCells();
1342 if (ncells<=0)
1343 return;
1344
1345 Double_t cellMeanE = 0, cellSigE = 0;
1346 for (Int_t i = 0; i<ncells; ++i) {
1347 Double_t cellE = cells->GetAmplitude(i);
1348 cellMeanE += cellE;
1349 cellSigE += cellE*cellE;
1350 }
1351 cellMeanE /= ncells;
1352 cellSigE /= ncells;
1353 cellSigE -= cellMeanE*cellMeanE;
1354 if (cellSigE<0)
1355 cellSigE = 0;
1356 cellSigE = TMath::Sqrt(cellSigE / ncells);
1357
1358 Double_t subE = cellMeanE - 7*cellSigE;
1359 if (subE<0)
1360 return;
1361
1362 for (Int_t i = 0; i<ncells; ++i) {
60d77596 1363 Short_t id=-1;
1364 Int_t mclabel = -1;
77e93dc2 1365 Double_t amp=0,time=0, efrac = 0;
1366 if (!cells->GetCell(i, id, amp, time,mclabel,efrac))
323834f0 1367 continue;
1368 amp -= cellMeanE;
1369 if (amp<0.001)
1370 amp = 0;
1371 cells->SetCell(i, id, amp, time);
1372 }
1373
1374 TObjArray *clusters = fEsdClusters;
1375 if (!clusters)
1376 clusters = fAodClusters;
323834f0 1377 if (!clusters)
1378 return;
1379
1380 Int_t nclus = clusters->GetEntries();
1381 for (Int_t i = 0; i<nclus; ++i) {
1382 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1383 if (!clus->IsEMCAL())
1384 continue;
1385 Int_t nc = clus->GetNCells();
1386 Double_t clusE = 0;
1387 UShort_t ids[100] = {0};
1388 Double_t fra[100] = {0};
1389 for (Int_t j = 0; j<nc; ++j) {
1390 Short_t id = TMath::Abs(clus->GetCellAbsId(j));
1391 Double_t cen = cells->GetCellAmplitude(id);
1392 clusE += cen;
1393 if (cen>0) {
1394 ids[nc] = id;
1395 ++nc;
1396 }
1397 }
1398 if (clusE<=0) {
1399 clusters->RemoveAt(i);
1400 continue;
1401 }
1402
1403 for (Int_t j = 0; j<nc; ++j) {
1404 Short_t id = ids[j];
1405 Double_t cen = cells->GetCellAmplitude(id);
1406 fra[j] = cen/clusE;
1407 }
1408 clus->SetE(clusE);
1409 AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus);
1410 if (aodclus) {
1411 aodclus->Clear("");
1412 aodclus->SetNCells(nc);
1413 aodclus->SetCellsAmplitudeFraction(fra);
1414 aodclus->SetCellsAbsId(ids);
1415 }
1416 }
1417 clusters->Compress();
1418}
1419
286b47a5 1420//________________________________________________________________________
1421void AliAnalysisTaskEMCALPi0PbPb::FillCellHists()
1422{
1423 // Fill histograms related to cell properties.
d595acbb 1424
90d5b88b 1425 AliVCaloCells *cells = fEsdCells;
1426 if (!cells)
1427 cells = fAodCells;
1428
296ea9b4 1429 if (!cells)
1430 return;
90d5b88b 1431
296ea9b4 1432 Int_t cellModCount[12] = {0};
1433 Double_t cellMaxE = 0;
1434 Double_t cellMeanE = 0;
1435 Int_t ncells = cells->GetNumberOfCells();
1436 if (ncells==0)
1437 return;
1438
1439 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
1440
1441 for (Int_t i = 0; i<ncells; ++i) {
1442 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
1443 Double_t cellE = cells->GetAmplitude(i);
1444 fHCellE->Fill(cellE);
1445 if (cellE>cellMaxE)
1446 cellMaxE = cellE;
1447 cellMeanE += cellE;
1448
1449 Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1;
1450 Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta);
1451 if (!ret) {
1452 AliError(Form("Could not get cell index for %d", absID));
1453 continue;
1454 }
1455 ++cellModCount[iSM];
1456 Int_t iPhi=-1, iEta=-1;
1457 fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta);
1458 fHColuRow[iSM]->Fill(iEta,iPhi,1);
1459 fHColuRowE[iSM]->Fill(iEta,iPhi,cellE);
1460 fHCellFreqNoCut[iSM]->Fill(absID);
2e4d8148 1461 if (cellE > 0.1) fHCellFreqCut100M[iSM]->Fill(absID);
1462 if (cellE > 0.3) fHCellFreqCut300M[iSM]->Fill(absID);
296ea9b4 1463 if (fHCellCheckE && fHCellCheckE[absID])
1464 fHCellCheckE[absID]->Fill(cellE);
2e4d8148 1465 fHCellFreqE[iSM]->Fill(absID, cellE);
296ea9b4 1466 }
1467 fHCellH->Fill(cellMaxE);
1468 cellMeanE /= ncells;
1469 fHCellM->Fill(cellMeanE);
1470 fHCellM2->Fill(cellMeanE*ncells/24/48/nsm);
1471 for (Int_t i=0; i<nsm; ++i)
1472 fHCellMult[i]->Fill(cellModCount[i]);
286b47a5 1473}
1474
1475//________________________________________________________________________
1476void AliAnalysisTaskEMCALPi0PbPb::FillClusHists()
1477{
90d5b88b 1478 // Fill histograms related to cluster properties.
fa443410 1479
90d5b88b 1480 TObjArray *clusters = fEsdClusters;
1481 if (!clusters)
1482 clusters = fAodClusters;
296ea9b4 1483 if (!clusters)
1484 return;
90d5b88b 1485
296ea9b4 1486 Double_t vertex[3] = {0};
1487 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1488
1489 Int_t nclus = clusters->GetEntries();
1490 for(Int_t i = 0; i<nclus; ++i) {
1491 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1492 if (!clus)
1493 continue;
1494 if (!clus->IsEMCAL())
1495 continue;
90d5b88b 1496 TLorentzVector clusterVec;
296ea9b4 1497 clus->GetMomentum(clusterVec,vertex);
3a952328 1498 Double_t maxAxis = clus->GetTOF(); //sigma
1499 Double_t clusterEcc = clus->Chi2(); //eccentricity
296ea9b4 1500 fHClustEccentricity->Fill(clusterEcc);
1501 fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi());
1502 fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt());
1503 fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E());
1504 fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis);
1505 fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E());
f5e0f1e2 1506 fHClustEnergyNCell->Fill(clus->E(),clus->GetNCells());
788ca675 1507 }
1508}
b3ee6797 1509
38727e64 1510//________________________________________________________________________
1511void AliAnalysisTaskEMCALPi0PbPb::CalcMcInfo()
1512{
1513 // Get Mc truth particle information.
1514
cfd7d5b2 1515 if (!fMcParts)
38727e64 1516 return;
1517
807016ea 1518 fMcParts->Clear();
1519
1520 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
1521 Double_t etamin = -0.7;
1522 Double_t etamax = +0.7;
1523 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
1524 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
1525
56fd6cb2 1526 if (fAodEv) {
807016ea 1527 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1528 am->LoadBranch(AliAODMCParticle::StdBranchName());
56fd6cb2 1529 TClonesArray *tca = dynamic_cast<TClonesArray*>(fAodEv->FindListObject(AliAODMCParticle::StdBranchName()));
1530 if (!tca)
1531 return;
1532
1533 Int_t nents = tca->GetEntries();
807016ea 1534 for(int it=0; it<nents; ++it) {
56fd6cb2 1535 AliAODMCParticle *part = static_cast<AliAODMCParticle*>(tca->At(it));
807016ea 1536 part->Print();
56fd6cb2 1537
b5ba4932 1538 // pion or eta meson or direct photon
56fd6cb2 1539 if(part->GetPdgCode() == 111) {
1540 } else if(part->GetPdgCode() == 221) {
b5ba4932 1541 } else if(part->GetPdgCode() == 22 ) {
1542 } else
56fd6cb2 1543 continue;
1544
1545 // primary particle
1546 Double_t dR = TMath::Sqrt((part->Xv()*part->Xv())+(part->Yv()*part->Yv()));
1547 if(dR > 1.0)
1548 continue;
56fd6cb2 1549
807016ea 1550 // kinematic cuts
1551 Double_t pt = part->Pt() ;
1552 if (pt<0.5)
1553 continue;
1554 Double_t eta = part->Eta();
1555 if (eta<etamin||eta>etamax)
1556 continue;
1557 Double_t phi = part->Phi();
1558 if (phi<phimin||phi>phimax)
1559 continue;
1560
1561 ProcessDaughters(part, it, tca);
56fd6cb2 1562 }
807016ea 1563 return;
1564 }
38727e64 1565
1566 AliMCEvent *mcEvent = MCEvent();
1567 if (!mcEvent)
1568 return;
1569
1570 const AliVVertex *evtVtx = mcEvent->GetPrimaryVertex();
1571 if (!evtVtx)
1572 return;
1573
1574 mcEvent->PreReadAll();
38727e64 1575
1576 Int_t nTracks = mcEvent->GetNumberOfPrimaries();
1577 for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
1578 AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
1579 if (!mcP)
1580 continue;
1581
b5ba4932 1582 // pion or eta meson or direct photon
807016ea 1583 if(mcP->PdgCode() == 111) {
1584 } else if(mcP->PdgCode() == 221) {
b5ba4932 1585 } else if(mcP->PdgCode() == 22 ) {
38727e64 1586 } else
1587 continue;
1588
1589 // primary particle
807016ea 1590 Double_t dR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
1591 (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
38727e64 1592 if(dR > 1.0)
1593 continue;
1594
38727e64 1595 // kinematic cuts
807016ea 1596 Double_t pt = mcP->Pt() ;
1597 if (pt<0.5)
38727e64 1598 continue;
807016ea 1599 Double_t eta = mcP->Eta();
1600 if (eta<etamin||eta>etamax)
38727e64 1601 continue;
807016ea 1602 Double_t phi = mcP->Phi();
1603 if (phi<phimin||phi>phimax)
1604 continue;
1605
1606 ProcessDaughters(mcP, iTrack, mcEvent);
38727e64 1607 }
1608}
1609
788ca675 1610//________________________________________________________________________
1611void AliAnalysisTaskEMCALPi0PbPb::FillNtuple()
1612{
1613 // Fill ntuple.
1614
1615 if (!fNtuple)
1616 return;
1617
1618 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1619 if (fAodEv) {
0a918d8d 1620 AliAODHeader * aodheader = dynamic_cast<AliAODHeader*>(fAodEv->GetHeader());
1621 if(!aodheader) AliFatal("Not a standard AOD");
1622
788ca675 1623 fHeader->fRun = fAodEv->GetRunNumber();
0a918d8d 1624 fHeader->fOrbit = aodheader->GetOrbitNumber();
1625 fHeader->fPeriod = aodheader->GetPeriodNumber();
1626 fHeader->fBx = aodheader->GetBunchCrossNumber();
1627 fHeader->fL0 = aodheader->GetL0TriggerInputs();
1628 fHeader->fL1 = aodheader->GetL1TriggerInputs();
1629 fHeader->fL2 = aodheader->GetL2TriggerInputs();
1630 fHeader->fTrClassMask = aodheader->GetTriggerMask();
1631 fHeader->fTrCluster = aodheader->GetTriggerCluster();
1632 fHeader->fOffTriggers = aodheader->GetOfflineTrigger();
1633 fHeader->fFiredTriggers = aodheader->GetFiredTriggerClasses();
788ca675 1634 } else {
1635 fHeader->fRun = fEsdEv->GetRunNumber();
1636 fHeader->fOrbit = fEsdEv->GetHeader()->GetOrbitNumber();
1637 fHeader->fPeriod = fEsdEv->GetESDRun()->GetPeriodNumber();
1638 fHeader->fBx = fEsdEv->GetHeader()->GetBunchCrossNumber();
1639 fHeader->fL0 = fEsdEv->GetHeader()->GetL0TriggerInputs();
1640 fHeader->fL1 = fEsdEv->GetHeader()->GetL1TriggerInputs();
1641 fHeader->fL2 = fEsdEv->GetHeader()->GetL2TriggerInputs();
1642 fHeader->fTrClassMask = fEsdEv->GetHeader()->GetTriggerMask();
1643 fHeader->fTrCluster = fEsdEv->GetHeader()->GetTriggerCluster();
1644 fHeader->fOffTriggers = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
1645 fHeader->fFiredTriggers = fEsdEv->GetFiredTriggerClasses();
5fe1ca23 1646 Float_t v0CorrR = 0;
1647 fHeader->fV0 = AliESDUtils::GetCorrV0(fEsdEv,v0CorrR);
1648 const AliMultiplicity *mult = fEsdEv->GetMultiplicity();
1649 if (mult)
1650 fHeader->fCl1 = mult->GetNumberOfITSClusters(1);
1651 fHeader->fTr = AliESDtrackCuts::GetReferenceMultiplicity(fEsdEv,1);
c0d2e1f2 1652 AliTriggerAnalysis trAn; /// Trigger Analysis
1653 Bool_t v0B = trAn.IsOfflineTriggerFired(fEsdEv, AliTriggerAnalysis::kV0C);
1654 Bool_t v0A = trAn.IsOfflineTriggerFired(fEsdEv, AliTriggerAnalysis::kV0A);
469b2bff 1655 fHeader->fV0And = v0A && v0B;
1656 fHeader->fIsHT = (fHeader->fOffTriggers & AliVEvent::kEMC1) || (fHeader->fOffTriggers & AliVEvent::kEMC7);
1657 am->LoadBranch("SPDPileupVertices");
1658 am->LoadBranch("TrkPileupVertices");
1659 fHeader->fIsPileup = fEsdEv->IsPileupFromSPD(3,0.8);
1660 fHeader->fIsPileup2 = fEsdEv->IsPileupFromSPD(3,0.4);
1661 fHeader->fIsPileup4 = fEsdEv->IsPileupFromSPD(3,0.2);
1662 fHeader->fIsPileup8 = fEsdEv->IsPileupFromSPD(3,0.1);
1663 fHeader->fNSpdVertices = fEsdEv->GetNumberOfPileupVerticesSPD();
1664 fHeader->fNTpcVertices = fEsdEv->GetNumberOfPileupVerticesTracks();
788ca675 1665 }
2ee7bda4 1666
788ca675 1667 AliCentrality *cent = InputEvent()->GetCentrality();
1668 fHeader->fV0Cent = cent->GetCentralityPercentileUnchecked("V0M");
1669 fHeader->fCl1Cent = cent->GetCentralityPercentileUnchecked("CL1");
1670 fHeader->fTrCent = cent->GetCentralityPercentileUnchecked("TRK");
1671 fHeader->fCqual = cent->GetQuality();
b6c599fe 1672
1673 AliEventplane *ep = InputEvent()->GetEventplane();
1674 if (ep) {
1675 if (ep->GetQVector())
1676 fHeader->fPsi = ep->GetQVector()->Phi()/2. ;
f2b8fca6 1677 else
1678 fHeader->fPsi = -1;
b6c599fe 1679 if (ep->GetQsub1()&&ep->GetQsub2())
1680 fHeader->fPsiRes = ep->GetQsub1()->Phi()/2.-ep->GetQsub2()->Phi()/2.;
1681 else
1682 fHeader->fPsiRes = 0;
1683 }
1684
788ca675 1685 Double_t val = 0;
1686 TString trgclasses(fHeader->fFiredTriggers);
1687 for (Int_t j = 0; j<fTrClassNamesArr->GetEntries(); ++j) {
1688 const char *name = fTrClassNamesArr->At(j)->GetName();
2ee7bda4 1689 TRegexp regexp(name);
1690 if (trgclasses.Contains(regexp))
788ca675 1691 val += TMath::Power(2,j);
1692 }
1693 fHeader->fTcls = (UInt_t)val;
1694
5fe1ca23 1695 fHeader->fNSelTr = fSelTracks->GetEntries();
1696 fHeader->fNSelPrimTr = fSelPrimTracks->GetEntries();
f5e0f1e2 1697 fHeader->fNSelPrimTr1 = 0;
1698 fHeader->fNSelPrimTr2 = 0;
1699 for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); iTracks++){
1700 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1701 if(track->Pt()>1)
1702 ++fHeader->fNSelPrimTr1;
1703 if(track->Pt()>2)
1704 ++fHeader->fNSelPrimTr2;
1705 }
5fe1ca23 1706
1707 fHeader->fNCells = 0;
c0d2e1f2 1708 fHeader->fNCells0 = 0;
1709 fHeader->fNCells01 = 0;
1710 fHeader->fNCells03 = 0;
5fe1ca23 1711 fHeader->fNCells1 = 0;
1712 fHeader->fNCells2 = 0;
1713 fHeader->fNCells5 = 0;
1714 fHeader->fMaxCellE = 0;
1715
1716 AliVCaloCells *cells = fEsdCells;
1717 if (!cells)
1718 cells = fAodCells;
1719
1720 if (cells) {
1721 Int_t ncells = cells->GetNumberOfCells();
1722 for(Int_t j=0; j<ncells; ++j) {
1723 Double_t cellen = cells->GetAmplitude(j);
c0d2e1f2 1724 if (cellen>0.045)
1725 ++fHeader->fNCells0;
1726 if (cellen>0.1)
1727 ++fHeader->fNCells01;
1728 if (cellen>0.3)
1729 ++fHeader->fNCells03;
5fe1ca23 1730 if (cellen>1)
1731 ++fHeader->fNCells1;
1732 if (cellen>2)
1733 ++fHeader->fNCells2;
1734 if (cellen>5)
1735 ++fHeader->fNCells5;
1736 if (cellen>fHeader->fMaxCellE)
1737 fHeader->fMaxCellE = cellen;
1738 }
1739 fHeader->fNCells = ncells;
1740 }
1741
1742 fHeader->fNClus = 0;
1743 fHeader->fNClus1 = 0;
1744 fHeader->fNClus2 = 0;
1745 fHeader->fNClus5 = 0;
1746 fHeader->fMaxClusE = 0;
1747
1748 TObjArray *clusters = fEsdClusters;
1749 if (!clusters)
1750 clusters = fAodClusters;
1751
1752 if (clusters) {
1753 Int_t nclus = clusters->GetEntries();
1754 for(Int_t j=0; j<nclus; ++j) {
1755 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(j));
6d4faab0 1756 if (!clus->IsEMCAL())
1757 continue;
5fe1ca23 1758 Double_t clusen = clus->E();
1759 if (clusen>1)
1760 ++fHeader->fNClus1;
1761 if (clusen>2)
1762 ++fHeader->fNClus2;
1763 if (clusen>5)
1764 ++fHeader->fNClus5;
1765 if (clusen>fHeader->fMaxClusE)
1766 fHeader->fMaxClusE = clusen;
1767 }
1768 fHeader->fNClus = nclus;
1769 }
1770
2ee7bda4 1771 fHeader->fMaxTrE = 0;
1772 if (fTriggers) {
1773 Int_t ntrig = fTriggers->GetEntries();
1774 for (Int_t j = 0; j<ntrig; ++j) {
1775 AliStaTrigger *sta = static_cast<AliStaTrigger*>(fTriggers->At(j));
1776 if (!sta)
1777 continue;
1778 if (sta->fE>fHeader->fMaxTrE)
1779 fHeader->fMaxTrE = sta->fE;
1780 }
1781 }
1782
1783 // count cells above 100 MeV on super modules
1784 fHeader->fNcSM0 = GetNCells(0, 0.100);
1785 fHeader->fNcSM1 = GetNCells(1, 0.100);
1786 fHeader->fNcSM2 = GetNCells(2, 0.100);
1787 fHeader->fNcSM3 = GetNCells(3, 0.100);
1788 fHeader->fNcSM4 = GetNCells(4, 0.100);
1789 fHeader->fNcSM5 = GetNCells(5, 0.100);
1790 fHeader->fNcSM6 = GetNCells(6, 0.100);
1791 fHeader->fNcSM7 = GetNCells(7, 0.100);
1792 fHeader->fNcSM8 = GetNCells(8, 0.100);
1793 fHeader->fNcSM9 = GetNCells(9, 0.100);
1794
788ca675 1795 if (fAodEv) {
1796 am->LoadBranch("vertices");
1797 AliAODVertex *pv = fAodEv->GetPrimaryVertex();
1798 FillVertex(fPrimVert, pv);
1799 AliAODVertex *sv = fAodEv->GetPrimaryVertexSPD();
1800 FillVertex(fSpdVert, sv);
1801 } else {
1802 am->LoadBranch("PrimaryVertex.");
1803 const AliESDVertex *pv = fEsdEv->GetPrimaryVertexTracks();
1804 FillVertex(fPrimVert, pv);
1805 am->LoadBranch("SPDVertex.");
1806 const AliESDVertex *sv = fEsdEv->GetPrimaryVertexSPD();
1807 FillVertex(fSpdVert, sv);
1808 am->LoadBranch("TPCVertex.");
1809 const AliESDVertex *tv = fEsdEv->GetPrimaryVertexTPC();
1810 FillVertex(fTpcVert, tv);
fa443410 1811 }
788ca675 1812
788ca675 1813 fNtuple->Fill();
286b47a5 1814}
1815
1816//________________________________________________________________________
1817void AliAnalysisTaskEMCALPi0PbPb::FillPionHists()
1818{
1819 // Fill histograms related to pions.
286b47a5 1820
296ea9b4 1821 Double_t vertex[3] = {0};
fa443410 1822 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1823
90d5b88b 1824 TObjArray *clusters = fEsdClusters;
1825 if (!clusters)
1826 clusters = fAodClusters;
fa443410 1827
90d5b88b 1828 if (clusters) {
1829 TLorentzVector clusterVec1;
1830 TLorentzVector clusterVec2;
1831 TLorentzVector pionVec;
1832
1833 Int_t nclus = clusters->GetEntries();
fa443410 1834 for (Int_t i = 0; i<nclus; ++i) {
90d5b88b 1835 AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i));
fa443410 1836 if (!clus1)
1837 continue;
1838 if (!clus1->IsEMCAL())
1839 continue;
3c661da5 1840 if (clus1->E()<fMinE)
6eb6260e 1841 continue;
a49742b5 1842 if (clus1->GetNCells()<fNminCells)
1843 continue;
f224d35b 1844 if (GetMaxCellEnergy(clus1)/clus1->E()<fMinErat)
1845 continue;
3c661da5 1846 if (clus1->Chi2()<fMinEcc) // eccentricity cut
f224d35b 1847 continue;
fa443410 1848 clus1->GetMomentum(clusterVec1,vertex);
1849 for (Int_t j = i+1; j<nclus; ++j) {
90d5b88b 1850 AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j));
fa443410 1851 if (!clus2)
1852 continue;
1853 if (!clus2->IsEMCAL())
1854 continue;
3c661da5 1855 if (clus2->E()<fMinE)
6eb6260e 1856 continue;
a49742b5 1857 if (clus2->GetNCells()<fNminCells)
1858 continue;
f224d35b 1859 if (GetMaxCellEnergy(clus2)/clus2->E()<fMinErat)
1860 continue;
3c661da5 1861 if (clus2->Chi2()<fMinEcc) // eccentricity cut
f224d35b 1862 continue;
fa443410 1863 clus2->GetMomentum(clusterVec2,vertex);
1864 pionVec = clusterVec1 + clusterVec2;
1865 Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E();
6eb6260e 1866 Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect());
d595acbb 1867 if (pionZgg < fAsymMax) {
1868 fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi());
1869 fHPionMggPt->Fill(pionVec.M(),pionVec.Pt());
1870 fHPionMggAsym->Fill(pionVec.M(),pionZgg);
6eb6260e 1871 fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle);
d595acbb 1872 Int_t bin = fPtRanges->FindBin(pionVec.Pt());
1873 fHPionInvMasses[bin]->Fill(pionVec.M());
1874 }
fa443410 1875 }
1876 }
90d5b88b 1877 }
fa443410 1878}
d595acbb 1879
38727e64 1880//________________________________________________________________________
1881void AliAnalysisTaskEMCALPi0PbPb::FillMcHists()
1882{
1883 // Fill additional MC information histograms.
1884
cfd7d5b2 1885 if (!fMcParts)
38727e64 1886 return;
1887
1888 // check if aod or esd mc mode and the fill histos
1889}
1890
6eb6260e 1891//________________________________________________________________________
323834f0 1892void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists()
6eb6260e 1893{
788ca675 1894 // Fill other histograms.
b5ba4932 1895
1896 for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); ++iTracks){
1897 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1898 if(!track)
1899 continue;
1900 fHPrimTrackPt->Fill(track->Pt());
1901 fHPrimTrackEta->Fill(track->Eta());
1902 fHPrimTrackPhi->Fill(track->Phi());
1903 }
788ca675 1904}
1905
f5e0f1e2 1906//________________________________________________________________________
1907void AliAnalysisTaskEMCALPi0PbPb::FillTrackHists()
1908{
1909 // Fill track histograms.
1910
1911 if (fSelPrimTracks) {
1912 for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); iTracks++) {
1913 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1914 if(!track)
1915 continue;
1916 fHPrimTrackPt->Fill(track->Pt());
1917 fHPrimTrackEta->Fill(track->Eta());
1918 fHPrimTrackPhi->Fill(track->Phi());
1919 }
1920 }
1921}
1922
788ca675 1923//__________________________________________________________________________________________________
1924void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliESDVertex *esdv)
1925{
1926 // Fill vertex from ESD vertex info.
1927
1928 v->fVx = esdv->GetX();
1929 v->fVy = esdv->GetY();
1930 v->fVz = esdv->GetZ();
1931 v->fVc = esdv->GetNContributors();
1932 v->fDisp = esdv->GetDispersion();
1933 v->fZres = esdv->GetZRes();
1934 v->fChi2 = esdv->GetChi2();
1935 v->fSt = esdv->GetStatus();
1936 v->fIs3D = esdv->IsFromVertexer3D();
1937 v->fIsZ = esdv->IsFromVertexerZ();
1938}
1939
1940//__________________________________________________________________________________________________
1941void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliAODVertex *aodv)
1942{
1943 // Fill vertex from AOD vertex info.
1944
1945 v->fVx = aodv->GetX();
1946 v->fVy = aodv->GetY();
1947 v->fVz = aodv->GetZ();
1948 v->fVc = aodv->GetNContributors();
1949 v->fChi2 = aodv->GetChi2();
6eb6260e 1950}
1951
d595acbb 1952//________________________________________________________________________
296ea9b4 1953Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
1954{
1955 // Compute isolation based on cell content.
1956
1957 AliVCaloCells *cells = fEsdCells;
1958 if (!cells)
1959 cells = fAodCells;
1960 if (!cells)
1961 return 0;
1962
1963 Double_t cellIsolation = 0;
1964 Double_t rad2 = radius*radius;
1965 Int_t ncells = cells->GetNumberOfCells();
1966 for (Int_t i = 0; i<ncells; ++i) {
1967 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
0fbe8d4f 1968 Float_t eta=-1, phi=-1;
296ea9b4 1969 fGeom->EtaPhiFromIndex(absID,eta,phi);
0fbe8d4f 1970 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
296ea9b4 1971 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1972 if(dist>rad2)
1973 continue;
0fbe8d4f 1974 Double_t cellE = cells->GetAmplitude(i);
2ee7bda4 1975 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
1976 Double_t cellEt = cellE*sin(theta);
1977 cellIsolation += cellEt;
296ea9b4 1978 }
1979 return cellIsolation;
1980}
2ee7bda4 1981
5fc7508c 1982//________________________________________________________________________
1983Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsoNxM(Double_t cEta, Double_t cPhi, Int_t N, Int_t M) const
1984{
1985 // Compute isolation based on cell content, in a NxM rectangle.
1986
1987 AliVCaloCells *cells = fEsdCells;
1988 if (!cells)
1989 cells = fAodCells;
1990 if (!cells)
1991 return 0;
1992
1993 Double_t cellIsolation = 0;
1994 Int_t ncells = cells->GetNumberOfCells();
1995 for (Int_t i = 0; i<ncells; ++i) {
1996 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
1997 Float_t eta=-1, phi=-1;
1998 fGeom->EtaPhiFromIndex(absID,eta,phi);
1999 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
2000 Double_t etadiff = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
2001 if(TMath::Abs(etadiff)/0.014>N)
2002 continue;
2003 if(TMath::Abs(phidiff)/0.014>M)
2004 continue;
2005 Double_t cellE = cells->GetAmplitude(i);
2ee7bda4 2006 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
2007 Double_t cellEt = cellE*sin(theta);
2008 cellIsolation += cellEt;
5fc7508c 2009 }
2010 return cellIsolation;
2011}
2012
296ea9b4 2013//________________________________________________________________________
0fbe8d4f 2014Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellEnergy(const AliVCluster *cluster) const
2015{
2016 // Get maximum energy of attached cell.
2017
2018 Double_t ret = 0;
2019 Int_t ncells = cluster->GetNCells();
2020 if (fEsdCells) {
2021 for (Int_t i=0; i<ncells; i++) {
2022 Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2023 ret += e;
2024 }
2025 } else {
2026 for (Int_t i=0; i<ncells; i++) {
2027 Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2028 ret += e;
2029 }
2030 }
2031 return ret;
2032}
2033
2034//________________________________________________________________________
2035Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
d595acbb 2036{
90d5b88b 2037 // Get maximum energy of attached cell.
2038
788ca675 2039 id = -1;
90d5b88b 2040 Double_t maxe = 0;
90d5b88b 2041 Int_t ncells = cluster->GetNCells();
2042 if (fEsdCells) {
2043 for (Int_t i=0; i<ncells; i++) {
27422847 2044 Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
f0897c18 2045 if (e>maxe) {
90d5b88b 2046 maxe = e;
788ca675 2047 id = cluster->GetCellAbsId(i);
f0897c18 2048 }
90d5b88b 2049 }
2050 } else {
2051 for (Int_t i=0; i<ncells; i++) {
27422847 2052 Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
90d5b88b 2053 if (e>maxe)
2054 maxe = e;
788ca675 2055 id = cluster->GetCellAbsId(i);
90d5b88b 2056 }
2057 }
6eb6260e 2058 return maxe;
d595acbb 2059}
2060
469b2bff 2061//________________________________________________________________________
1f41ed3d 2062Double_t AliAnalysisTaskEMCALPi0PbPb::GetSecondMaxCellEnergy(AliVCluster *clus, Short_t &id) const
469b2bff 2063{
2064 // Get second maximum cell.
2065
2066 AliVCaloCells *cells = fEsdCells;
2067 if (!cells)
2068 cells = fAodCells;
2069 if (!cells)
2070 return -1;
2071
2072 Double_t secondEmax=0, firstEmax=0;
2073 Double_t cellen;
2074 for(Int_t iCell=0;iCell<clus->GetNCells();iCell++){
2075 Int_t absId = clus->GetCellAbsId(iCell);
2076 cellen = cells->GetCellAmplitude(absId);
2077 if(cellen > firstEmax)
2078 firstEmax = cellen;
2079 }
2080 for(Int_t iCell=0;iCell<clus->GetNCells();iCell++){
2081 Int_t absId = clus->GetCellAbsId(iCell);
2082 cellen = cells->GetCellAmplitude(absId);
1f41ed3d 2083 if(cellen < firstEmax && cellen > secondEmax) {
469b2bff 2084 secondEmax = cellen;
1f41ed3d 2085 id = absId;
2086 }
469b2bff 2087 }
2088 return secondEmax;
2089}
2090
d595acbb 2091//________________________________________________________________________
0fbe8d4f 2092void AliAnalysisTaskEMCALPi0PbPb::GetSigma(const AliVCluster *c, Double_t& sigmaMax, Double_t &sigmaMin) const
d595acbb 2093{
90d5b88b 2094 // Calculate the (E) weighted variance along the longer (eigen) axis.
2095
6bf90832 2096 sigmaMax = 0; // cluster variance along its longer axis
2097 sigmaMin = 0; // cluster variance along its shorter axis
2098 Double_t Ec = c->E(); // cluster energy
296ea9b4 2099 if(Ec<=0)
2100 return;
6bf90832 2101 Double_t Xc = 0 ; // cluster first moment along X
2102 Double_t Yc = 0 ; // cluster first moment along Y
2103 Double_t Sxx = 0 ; // cluster second central moment along X (variance_X^2)
2104 Double_t Sxy = 0 ; // cluster second central moment along Y (variance_Y^2)
2105 Double_t Syy = 0 ; // cluster covariance^2
90d5b88b 2106
2107 AliVCaloCells *cells = fEsdCells;
2108 if (!cells)
2109 cells = fAodCells;
2110
6eb6260e 2111 if (!cells)
2112 return;
2113
6bf90832 2114 Int_t ncells = c->GetNCells();
6eb6260e 2115 if (ncells==1)
2116 return;
2117
6eb6260e 2118 for(Int_t j=0; j<ncells; ++j) {
6bf90832 2119 Int_t id = TMath::Abs(c->GetCellAbsId(j));
6eb6260e 2120 Double_t cellen = cells->GetCellAmplitude(id);
3a952328 2121 TVector3 pos;
2122 fGeom->GetGlobal(id,pos);
6eb6260e 2123 Xc += cellen*pos.X();
2124 Yc += cellen*pos.Y();
2125 Sxx += cellen*pos.X()*pos.X();
2126 Syy += cellen*pos.Y()*pos.Y();
2127 Sxy += cellen*pos.X()*pos.Y();
2128 }
2129 Xc /= Ec;
2130 Yc /= Ec;
2131 Sxx /= Ec;
2132 Syy /= Ec;
2133 Sxy /= Ec;
2134 Sxx -= Xc*Xc;
2135 Syy -= Yc*Yc;
2136 Sxy -= Xc*Yc;
f0897c18 2137 Sxx = TMath::Abs(Sxx);
2138 Syy = TMath::Abs(Syy);
296ea9b4 2139 sigmaMax = (Sxx + Syy + TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
2140 sigmaMax = TMath::Sqrt(TMath::Abs(sigmaMax));
2141 sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
2142 sigmaMin = TMath::Sqrt(TMath::Abs(sigmaMin));
d595acbb 2143}
2ee7bda4 2144
5fc7508c 2145//________________________________________________________________________
2146void AliAnalysisTaskEMCALPi0PbPb::GetSigmaEtaEta(const AliVCluster *c, Double_t& sEtaEta, Double_t &sPhiPhi) const
2147{
2148 // Calculate the (E) weighted variance along the pseudorapidity.
2149
2150 sEtaEta = 0;
2151 sPhiPhi = 0;
2152
2153 Double_t Ec = c->E(); // cluster energy
2154 if(Ec<=0)
2155 return;
2156
2157 const Int_t ncells = c->GetNCells();
2158
2159 Double_t EtaC = 0; // cluster first moment along eta
2160 Double_t PhiC = 0; // cluster first moment along phi
2161 Double_t Setaeta = 0; // cluster second central moment along eta
2162 Double_t Sphiphi = 0; // cluster second central moment along phi
2163 Double_t w[ncells]; // weight max(0,4.5*log(E_i/Ec))
2164 Double_t sumw = 0;
2165 Int_t id[ncells];
2166
2167 AliVCaloCells *cells = fEsdCells;
2168 if (!cells)
2169 cells = fAodCells;
2170
2171 if (!cells)
2172 return;
2173
2174 if (ncells==1)
2175 return;
2176
2177 for(Int_t j=0; j<ncells; ++j) {
2178 id[j] = TMath::Abs(c->GetCellAbsId(j));
2179 Double_t cellen = cells->GetCellAmplitude(id[j]);
2180 w[j] = TMath::Max(0., 4.5+TMath::Log(cellen/Ec));
2181 TVector3 pos;
2182 fGeom->GetGlobal(id[j],pos);
2183 EtaC += w[j]*pos.Eta();
2184 PhiC += w[j]*pos.Phi();
2185 sumw += w[j];
2186 }
2187 EtaC /= sumw;
2188 PhiC /= sumw;
2189
2190 for(Int_t j=0; j<ncells; ++j) {
2191 TVector3 pos;
2192 fGeom->GetGlobal(id[j],pos);
2193 Setaeta = w[j]*(pos.Eta() - EtaC)*(pos.Eta() - EtaC);
2194 Sphiphi = w[j]*(pos.Phi() - PhiC)*(pos.Phi() - PhiC);
2195 }
2196 Setaeta /= sumw;
2197 sEtaEta = TMath::Sqrt(Setaeta);
2198 Sphiphi /= sumw;
2199 sPhiPhi = TMath::Sqrt(Sphiphi);
2200}
2201
6bf90832 2202//________________________________________________________________________
0fbe8d4f 2203Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(const AliVCluster *c, Double_t emin) const
6bf90832 2204{
2205 // Calculate number of attached cells above emin.
2206
6bf90832 2207 AliVCaloCells *cells = fEsdCells;
2208 if (!cells)
2209 cells = fAodCells;
6bf90832 2210 if (!cells)
296ea9b4 2211 return 0;
6bf90832 2212
296ea9b4 2213 Int_t n = 0;
6bf90832 2214 Int_t ncells = c->GetNCells();
2215 for(Int_t j=0; j<ncells; ++j) {
2216 Int_t id = TMath::Abs(c->GetCellAbsId(j));
2217 Double_t cellen = cells->GetCellAmplitude(id);
2218 if (cellen>=emin)
2219 ++n;
2220 }
2221 return n;
2222}
2223
2ee7bda4 2224//________________________________________________________________________
2225Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(Int_t sm, Double_t emin) const
2226{
2227 // Calculate number of cells per SM above emin.
2228
2229 AliVCaloCells *cells = fEsdCells;
2230 if (!cells)
2231 cells = fAodCells;
2232 if (!cells)
2233 return 0;
2234
2235 Int_t n = 0;
2236 Int_t ncells = cells->GetNumberOfCells();
2237 for(Int_t j=0; j<ncells; ++j) {
2238 Int_t id = TMath::Abs(cells->GetCellNumber(j));
2239 Double_t cellen = cells->GetCellAmplitude(id);
2240 if (cellen<emin)
2241 continue;
2242 Int_t fsm = fGeom->GetSuperModuleNumber(id);
2243 if (fsm != sm)
2244 continue;
2245 ++n;
2246 }
2247 return n;
2248}
2249
296ea9b4 2250//________________________________________________________________________
b6c599fe 2251Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
296ea9b4 2252{
2253 // Compute isolation based on tracks.
2254
2255 Double_t trkIsolation = 0;
2256 Double_t rad2 = radius*radius;
b6c599fe 2257 Int_t ntrks = fSelPrimTracks->GetEntries();
296ea9b4 2258 for(Int_t j = 0; j<ntrks; ++j) {
1f41ed3d 2259 AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
296ea9b4 2260 if (!track)
2261 continue;
b6c599fe 2262 if (track->Pt()<pt)
2263 continue;
296ea9b4 2264 Float_t eta = track->Eta();
2265 Float_t phi = track->Phi();
0fbe8d4f 2266 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
296ea9b4 2267 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
2268 if(dist>rad2)
2269 continue;
2270 trkIsolation += track->Pt();
2271 }
2272 return trkIsolation;
2273}
2ee7bda4 2274
5fc7508c 2275//________________________________________________________________________
2276Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsoStrip(Double_t cEta, Double_t cPhi, Double_t dEta, Double_t dPhi, Double_t pt) const
2277{
2278 // Compute isolation based on tracks.
2279
2280 Double_t trkIsolation = 0;
2281 Int_t ntrks = fSelPrimTracks->GetEntries();
2282 for(Int_t j = 0; j<ntrks; ++j) {
1f41ed3d 2283 AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
5fc7508c 2284 if (!track)
2285 continue;
2286 if (track->Pt()<pt)
2287 continue;
2288 Float_t eta = track->Eta();
2289 Float_t phi = track->Phi();
2290 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
2291 Double_t etadiff = (eta-cEta);
2292 if(TMath::Abs(etadiff)>dEta || TMath::Abs(phidiff)>dPhi)
2293 continue;
2294 trkIsolation += track->Pt();
2295 }
2296 return trkIsolation;
2297}
2298
0fbe8d4f 2299//________________________________________________________________________
2300Bool_t AliAnalysisTaskEMCALPi0PbPb::IsShared(const AliVCluster *c) const
2301{
2302 // Returns if cluster shared across super modules.
2303
1f41ed3d 2304 if (!c)
0fbe8d4f 2305 return 0;
2306
2307 Int_t n = -1;
2308 Int_t ncells = c->GetNCells();
2309 for(Int_t j=0; j<ncells; ++j) {
2310 Int_t id = TMath::Abs(c->GetCellAbsId(j));
2311 Int_t got = id / (24*48);
2312 if (n==-1) {
2313 n = got;
2314 continue;
2315 }
2316 if (got!=n)
2317 return 1;
2318 }
2319 return 0;
2320}
3a952328 2321
2ee7bda4 2322//________________________________________________________________________
2323Bool_t AliAnalysisTaskEMCALPi0PbPb::IsIdPartOfCluster(const AliVCluster *c, Short_t id) const
2324{
2325 // Returns if id is part of cluster.
2326
2327 AliVCaloCells *cells = fEsdCells;
2328 if (!cells)
2329 cells = fAodCells;
2330 if (!cells)
2331 return 0;
2332
2333 Int_t ncells = c->GetNCells();
2334 for(Int_t j=0; j<ncells; ++j) {
2335 Int_t cid = TMath::Abs(c->GetCellAbsId(j));
2336 if (cid == id)
2337 return 1;
2338 }
2339 return 0;
2340}
2341
38727e64 2342//________________________________________________________________________
807016ea 2343void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(const AliVParticle *p, const TObjArray *arr, Int_t level) const
38727e64 2344{
56fd6cb2 2345 // Print recursively daughter information.
2346
2347 if (!p || !arr)
2348 return;
2349
807016ea 2350 const AliAODMCParticle *amc = dynamic_cast<const AliAODMCParticle*>(p);
2351 if (!amc)
2352 return;
2353 for (Int_t i=0; i<level; ++i) printf(" ");
2354 amc->Print();
2355
2356 Int_t n = amc->GetNDaughters();
2357 for (Int_t i=0; i<n; ++i) {
2358 Int_t d = amc->GetDaughter(i);
2359 const AliVParticle *dmc = static_cast<const AliVParticle*>(arr->At(d));
2360 PrintDaughters(dmc,arr,level+1);
2361 }
2362}
2363
2364//________________________________________________________________________
2365void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(const AliMCParticle *p, const AliMCEvent *arr, Int_t level) const
2366{
2367 // Print recursively daughter information.
2368
2369 if (!p || !arr)
2370 return;
2371
2372 for (Int_t i=0; i<level; ++i) printf(" ");
2373 Int_t d1 = p->GetFirstDaughter();
2374 Int_t d2 = p->GetLastDaughter();
2375 printf("pid=%d: %.2f %.2f %.2f (%.2f %.2f %.2f); nd=%d,%d\n",
2376 p->PdgCode(),p->Px(),p->Py(),p->Pz(),p->Xv(),p->Yv(),p->Zv(),d1,d2);
2377 if (d1<0)
2378 return;
2379 if (d2<0)
2380 d2=d1;
2381 for (Int_t i=d1;i<=d2;++i) {
2382 const AliMCParticle *dmc = static_cast<const AliMCParticle *>(arr->GetTrack(i));
6d4faab0 2383 PrintDaughters(dmc,arr,level+1);
56fd6cb2 2384 }
38727e64 2385}
2386
2387//________________________________________________________________________
2388void AliAnalysisTaskEMCALPi0PbPb::PrintTrackRefs(AliMCParticle *p) const
2389{
56fd6cb2 2390 // Print track ref array.
2391
2392 if (!p)
2393 return;
2394
2395 Int_t n = p->GetNumberOfTrackReferences();
2396 for (Int_t i=0; i<n; ++i) {
2397 AliTrackReference *ref = p->GetTrackReference(i);
bc8a45bd 2398 if (!ref)
2399 continue;
56fd6cb2 2400 ref->SetUserId(ref->DetectorId());
2401 ref->Print();
2402 }
38727e64 2403}
2404
807016ea 2405//________________________________________________________________________
2406void AliAnalysisTaskEMCALPi0PbPb::ProcessDaughters(AliVParticle *p, Int_t index, const TObjArray *arr)
2407{
2408 // Process and create daughters.
2409
2410 if (!p || !arr)
2411 return;
2412
2413 AliAODMCParticle *amc = dynamic_cast<AliAODMCParticle*>(p);
2414 if (!amc)
2415 return;
2416
2417 //amc->Print();
2418
2419 Int_t nparts = arr->GetEntries();
2420 Int_t nents = fMcParts->GetEntries();
2421
2422 AliStaPart *newp = static_cast<AliStaPart*>(fMcParts->New(nents));
2423 newp->fPt = amc->Pt();
2424 newp->fEta = amc->Eta();
2425 newp->fPhi = amc->Phi();
2426 if (amc->Xv() != 0 || amc->Yv() != 0 || amc->Zv() != 0) {
2427 TVector3 vec(amc->Xv(),amc->Yv(),amc->Zv());
2428 newp->fVR = vec.Perp();
2429 newp->fVEta = vec.Eta();
2430 newp->fVPhi = vec.Phi();
2431 }
2432 newp->fPid = amc->PdgCode();
8c56d760 2433 newp->fLab = nents;
807016ea 2434 Int_t moi = amc->GetMother();
2435 if (moi>=0&&moi<nparts) {
2436 const AliAODMCParticle *mmc = static_cast<const AliAODMCParticle*>(arr->At(moi));
2437 moi = mmc->GetUniqueID();
2438 }
2439 newp->fMo = moi;
2440 p->SetUniqueID(nents);
8c56d760 2441
2442 // TODO: Determine which detector was hit
2443 //newp->fDet = ???
2444
807016ea 2445 Int_t n = amc->GetNDaughters();
2446 for (Int_t i=0; i<n; ++i) {
2447 Int_t d = amc->GetDaughter(i);
2448 if (d<=index || d>=nparts)
2449 continue;
2450 AliVParticle *dmc = static_cast<AliVParticle*>(arr->At(d));
2451 ProcessDaughters(dmc,d,arr);
2452 }
2453}
2454
2455//________________________________________________________________________
2456void AliAnalysisTaskEMCALPi0PbPb::ProcessDaughters(AliMCParticle *p, Int_t index, const AliMCEvent *arr)
2457{
2458 // Process and create daughters.
2459
2460 if (!p || !arr)
2461 return;
2462
2463 Int_t d1 = p->GetFirstDaughter();
2464 Int_t d2 = p->GetLastDaughter();
2465 if (0) {
2466 printf("%d pid=%d: %.3f %.3f %.3f (%.2f %.2f %.2f); nd=%d,%d, mo=%d\n",
2467 index,p->PdgCode(),p->Px(),p->Py(),p->Pz(),p->Xv(),p->Yv(),p->Zv(),d1,d2, p->GetMother());
bc8a45bd 2468 PrintTrackRefs(p);
807016ea 2469 }
2470 Int_t nents = fMcParts->GetEntries();
2471
2472 AliStaPart *newp = static_cast<AliStaPart*>(fMcParts->New(nents));
2473 newp->fPt = p->Pt();
2474 newp->fEta = p->Eta();
2475 newp->fPhi = p->Phi();
2476 if (p->Xv() != 0 || p->Yv() != 0 || p->Zv() != 0) {
2477 TVector3 vec(p->Xv(),p->Yv(),p->Zv());
2478 newp->fVR = vec.Perp();
2479 newp->fVEta = vec.Eta();
2480 newp->fVPhi = vec.Phi();
2481 }
2482 newp->fPid = p->PdgCode();
8c56d760 2483 newp->fLab = nents;
807016ea 2484 Int_t moi = p->GetMother();
2485 if (moi>=0) {
2486 const AliMCParticle *mmc = static_cast<const AliMCParticle *>(arr->GetTrack(moi));
2487 moi = mmc->GetUniqueID();
2488 }
2489 newp->fMo = moi;
2490 p->SetUniqueID(nents);
2491
2492 Int_t nref = p->GetNumberOfTrackReferences();
2493 if (nref>0) {
2494 AliTrackReference *ref = p->GetTrackReference(nref-1);
2495 if (ref) {
2496 newp->fDet = ref->DetectorId();
2497 }
2498 }
2499
2500 if (d1<0)
2501 return;
2502 if (d2<0)
2503 d2=d1;
2504 for (Int_t i=d1;i<=d2;++i) {
2505 AliMCParticle *dmc = static_cast<AliMCParticle *>(arr->GetTrack(i));
2506 if (dmc->P()<0.01)
2507 continue;
2508 ProcessDaughters(dmc,i,arr);
2509 }
2510}
6dbf6b27 2511
2512//__________________________________________________________________________________________________
2513void AliStaCluster::GetMom(TLorentzVector& p, Double_t *vertex)
2514{
2515 // Calculate momentum.
2516
2517 TVector3 pos;
2518 pos.SetPtEtaPhi(fR,fEta,fPhi);
2519
a186e444 2520 if(vertex){ //calculate direction relative to vertex
6dbf6b27 2521 pos -= vertex;
2522 }
2523
2524 Double_t r = pos.Mag();
2525 p.SetPxPyPzE(fE*pos.x()/r, fE*pos.y()/r, fE*pos.z()/r, fE);
2526}
2527
2528//__________________________________________________________________________________________________
2529void AliStaCluster::GetMom(TLorentzVector& p, AliStaVertex *vertex)
2530{
2531 // Calculate momentum.
2532
2533 Double_t v[3] = {0,0,0};
2534 if (vertex) {
2535 v[0] = vertex->fVx;
2536 v[1] = vertex->fVy;
2537 v[2] = vertex->fVz;
2538 }
2539 GetMom(p, v);
2540}