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