]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/CaloCalib/AliAnalysisTaskEMCALPi0PbPb.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
improved handling of aods, added beginning of mc info
[u/mrichter/AliRoot.git] / PWG4 / CaloCalib / AliAnalysisTaskEMCALPi0PbPb.cxx
CommitLineData
ea3fd2d5 1// $Id$
2
ea3fd2d5 3#include "AliAnalysisTaskEMCALPi0PbPb.h"
fa443410 4#include <TAxis.h>
717fe7de 5#include <TChain.h>
6#include <TClonesArray.h>
f5d4ab70 7#include <TFile.h>
296ea9b4 8#include <TGeoGlobalMagField.h>
717fe7de 9#include <TH1F.h>
10#include <TH2F.h>
11#include <TList.h>
12#include <TLorentzVector.h>
f5d4ab70 13#include <TNtuple.h>
296ea9b4 14#include <TProfile.h>
b3ee6797 15#include <TString.h>
296ea9b4 16#include <TVector2.h>
717fe7de 17#include "AliAODEvent.h"
18#include "AliAODVertex.h"
19#include "AliAnalysisManager.h"
20#include "AliAnalysisTaskEMCALClusterizeFast.h"
296ea9b4 21#include "AliCDBManager.h"
ea3fd2d5 22#include "AliCentrality.h"
ea3fd2d5 23#include "AliEMCALGeoUtils.h"
0ec74551 24#include "AliEMCALGeometry.h"
0fbe8d4f 25#include "AliEMCALRecPoint.h"
296ea9b4 26#include "AliEMCALRecoUtils.h"
3a952328 27#include "AliESDCaloTrigger.h"
717fe7de 28#include "AliESDEvent.h"
5fe1ca23 29#include "AliESDUtils.h"
717fe7de 30#include "AliESDVertex.h"
296ea9b4 31#include "AliESDtrack.h"
0ec74551 32#include "AliESDtrackCuts.h"
b6c599fe 33#include "AliEventplane.h"
296ea9b4 34#include "AliGeomManager.h"
b3ee6797 35#include "AliInputEventHandler.h"
43cfaa06 36#include "AliLog.h"
38727e64 37#include "AliMCEvent.h"
38#include "AliMCParticle.h"
296ea9b4 39#include "AliMagF.h"
5fe1ca23 40#include "AliMultiplicity.h"
38727e64 41#include "AliStack.h"
0ec74551 42#include "AliTrackerBase.h"
ea3fd2d5 43
44ClassImp(AliAnalysisTaskEMCALPi0PbPb)
717fe7de 45
ea3fd2d5 46//________________________________________________________________________
47AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name)
48 : AliAnalysisTaskSE(name),
717fe7de 49 fCentVar("V0M"),
50 fCentFrom(0),
51 fCentTo(100),
76332037 52 fVtxZMin(-10),
53 fVtxZMax(+10),
43cfaa06 54 fUseQualFlag(1),
717fe7de 55 fClusName(),
f5d4ab70 56 fDoNtuple(0),
a49742b5 57 fDoAfterburner(0),
f224d35b 58 fAsymMax(1),
a49742b5 59 fNminCells(1),
3c661da5 60 fMinE(0.100),
f224d35b 61 fMinErat(0),
62 fMinEcc(0),
6bf90832 63 fGeoName("EMCAL_FIRSTYEARV1"),
b6c599fe 64 fMinNClusPerTr(50),
296ea9b4 65 fIsoDist(0.2),
b3ee6797 66 fTrClassNames(""),
0ec74551 67 fTrCuts(0),
b6c599fe 68 fPrimTrCuts(0),
69 fDoTrMatGeom(0),
3a952328 70 fTrainMode(0),
71 fMarkCells(),
72 fMinL0Time(-1),
73 fMaxL0Time(1024),
38727e64 74 fMcMode(0),
d595acbb 75 fGeom(0),
296ea9b4 76 fReco(0),
38727e64 77 fIsGeoMatsSet(0),
78 fNEvs(0),
717fe7de 79 fOutput(0),
b3ee6797 80 fTrClassNamesArr(0),
717fe7de 81 fEsdEv(0),
82 fAodEv(0),
43cfaa06 83 fRecPoints(0),
717fe7de 84 fEsdClusters(0),
85 fEsdCells(0),
86 fAodClusters(0),
286b47a5 87 fAodCells(0),
fa443410 88 fPtRanges(0),
296ea9b4 89 fSelTracks(0),
b6c599fe 90 fSelPrimTracks(0),
3a952328 91 fNAmpInTrigger(0),
92 fAmpInTrigger(0),
788ca675 93 fNtuple(0),
94 fHeader(0),
95 fPrimVert(0),
96 fSpdVert(0),
97 fTpcVert(0),
98 fClusters(0),
3a952328 99 fTriggers(0),
fa443410 100 fHCuts(0x0),
101 fHVertexZ(0x0),
76332037 102 fHVertexZ2(0x0),
fa443410 103 fHCent(0x0),
104 fHCentQual(0x0),
b3ee6797 105 fHTclsBeforeCuts(0x0),
106 fHTclsAfterCuts(0x0),
d595acbb 107 fHColuRow(0x0),
108 fHColuRowE(0x0),
109 fHCellMult(0x0),
110 fHCellE(0x0),
111 fHCellH(0x0),
6eb6260e 112 fHCellM(0x0),
a49742b5 113 fHCellM2(0x0),
296ea9b4 114 fHCellFreqNoCut(0x0),
2e4d8148 115 fHCellFreqCut100M(0x0),
116 fHCellFreqCut300M(0x0),
117 fHCellFreqE(0x0),
296ea9b4 118 fHCellCheckE(0x0),
6eb6260e 119 fHClustEccentricity(0),
fa443410 120 fHClustEtaPhi(0x0),
121 fHClustEnergyPt(0x0),
122 fHClustEnergySigma(0x0),
d595acbb 123 fHClustSigmaSigma(0x0),
6eb6260e 124 fHClustNCellEnergyRatio(0x0),
b6c599fe 125 fHMatchDr(0x0),
126 fHMatchDz(0x0),
127 fHMatchEp(0x0),
fa443410 128 fHPionEtaPhi(0x0),
129 fHPionMggPt(0x0),
6eb6260e 130 fHPionMggAsym(0x0),
131 fHPionMggDgg(0x0)
ea3fd2d5 132{
717fe7de 133 // Constructor.
ea3fd2d5 134
d595acbb 135 if (!name)
136 return;
137 SetName(name);
ea3fd2d5 138 DefineInput(0, TChain::Class());
ea3fd2d5 139 DefineOutput(1, TList::Class());
3a952328 140 fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,EMCALCells.,Tracks EMCALTrigger."
296ea9b4 141 "AOD:header,vertices,emcalCells,tracks";
ea3fd2d5 142}
ea3fd2d5 143
717fe7de 144//________________________________________________________________________
145AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb()
146{
147 // Destructor.
ea3fd2d5 148
b3ee6797 149 if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
150 delete fOutput; fOutput = 0;
151 }
fa443410 152 delete fPtRanges; fPtRanges = 0;
d595acbb 153 delete fGeom; fGeom = 0;
296ea9b4 154 delete fReco; fReco = 0;
b3ee6797 155 delete fTrClassNamesArr;
296ea9b4 156 delete fSelTracks;
b6c599fe 157 delete fSelPrimTracks;
3a952328 158 delete [] fAmpInTrigger;
d595acbb 159 delete [] fHColuRow;
160 delete [] fHColuRowE;
161 delete [] fHCellMult;
296ea9b4 162 delete [] fHCellFreqNoCut;
2e4d8148 163 delete [] fHCellFreqCut100M;
164 delete [] fHCellFreqCut300M;
ea3fd2d5 165}
717fe7de 166
ea3fd2d5 167//________________________________________________________________________
168void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects()
169{
717fe7de 170 // Create user objects here.
ea3fd2d5 171
f3582e89 172 cout << "AliAnalysisTaskEMCALPi0PbPb: Input settings" << endl;
173 cout << " fCentVar: " << fCentVar << endl;
174 cout << " fCentFrom: " << fCentFrom << endl;
175 cout << " fCentTo: " << fCentTo << endl;
176 cout << " fVtxZMin: " << fVtxZMin << endl;
177 cout << " fVtxZMax: " << fVtxZMax << endl;
178 cout << " fUseQualFlag: " << fUseQualFlag << endl;
179 cout << " fClusName: \"" << fClusName << "\"" << endl;
180 cout << " fDoNtuple: " << fDoNtuple << endl;
181 cout << " fDoAfterburner: " << fDoAfterburner << endl;
182 cout << " fAsymMax: " << fAsymMax << endl;
183 cout << " fNminCells: " << fNminCells << endl;
184 cout << " fMinE: " << fMinE << endl;
185 cout << " fMinErat: " << fMinErat << endl;
186 cout << " fMinEcc: " << fMinEcc << endl;
187 cout << " fGeoName: \"" << fGeoName << "\"" << endl;
188 cout << " fMinNClusPerTr: " << fMinNClusPerTr << endl;
38727e64 189 cout << " fIsoDist: " << fIsoDist << endl;
f3582e89 190 cout << " fTrClassNames: \"" << fTrClassNames << "\"" << endl;
f3582e89 191 cout << " fTrCuts: " << fTrCuts << endl;
192 cout << " fPrimTrCuts: " << fPrimTrCuts << endl;
193 cout << " fDoTrMatGeom: " << fDoTrMatGeom << endl;
194 cout << " fTrainMode: " << fTrainMode << endl;
195 cout << " fMarkCells: " << fMarkCells << endl;
196 cout << " fMinL0Time: " << fMinL0Time << endl;
197 cout << " fMaxL0Time: " << fMaxL0Time << endl;
38727e64 198 cout << " fMcMode: " << fMcMode << endl;
199 cout << " fGeom: " << fGeom << endl;
200 cout << " fReco: " << fReco << endl;
201
202 if (!fGeom)
203 fGeom = new AliEMCALGeoUtils(fGeoName,"EMCAL");
204 else {
205 if (fGeom->GetMatrixForSuperModule(0))
206 fIsGeoMatsSet = kTRUE;
207 }
208 if (!fReco)
209 fReco = new AliEMCALRecoUtils();
b3ee6797 210 fTrClassNamesArr = fTrClassNames.Tokenize(" ");
717fe7de 211 fOutput = new TList();
212 fOutput->SetOwner();
296ea9b4 213 fSelTracks = new TObjArray;
b6c599fe 214 fSelPrimTracks = new TObjArray;
ea3fd2d5 215
f5d4ab70 216 if (fDoNtuple) {
217 TFile *f = OpenFile(1);
6bf90832 218 if (f) {
219 f->SetCompressionLevel(2);
788ca675 220 fNtuple = new TTree(Form("tree%.0fto%.0f",fCentFrom,fCentTo), "StandaloneTree");
6bf90832 221 fNtuple->SetDirectory(f);
3a952328 222 if (fTrainMode) {
223 fNtuple->SetAutoFlush(-2*1024*1024);
224 fNtuple->SetAutoSave(0);
225 } else {
226 fNtuple->SetAutoFlush(-32*1024*1024);
227 fNtuple->SetAutoSave(0);
228 }
229
788ca675 230 fHeader = new AliStaHeader;
231 fNtuple->Branch("header", &fHeader, 16*1024, 99);
232 fPrimVert = new AliStaVertex;
233 fNtuple->Branch("primv", &fPrimVert, 16*1024, 99);
234 fSpdVert = new AliStaVertex;
235 fNtuple->Branch("spdv", &fSpdVert, 16*1024, 99);
236 fTpcVert = new AliStaVertex;
237 fNtuple->Branch("tpcv", &fTpcVert, 16*1024, 99);
238 if (TClass::GetClass("AliStaCluster"))
239 TClass::GetClass("AliStaCluster")->IgnoreTObjectStreamer();
3a952328 240 fClusters = new TClonesArray("AliStaCluster");
788ca675 241 fNtuple->Branch("clusters", &fClusters, 8*16*1024, 99);
3a952328 242 if (TClass::GetClass("AliStaTrigger"))
243 TClass::GetClass("AliStaTrigger")->IgnoreTObjectStreamer();
244 fTriggers = new TClonesArray("AliStaTrigger");
245 fNtuple->Branch("l0prim", &fTriggers, 16*1024, 99);
6bf90832 246 }
f5d4ab70 247 }
248
002dcebe 249 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
b6c599fe 250 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
251 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
002dcebe 252
d595acbb 253 // histograms
a49742b5 254 TH1::SetDefaultSumw2(kTRUE);
255 TH2::SetDefaultSumw2(kTRUE);
b3ee6797 256 fHCuts = new TH1F("hEventCuts","",5,0.5,5.5);
257 fHCuts->GetXaxis()->SetBinLabel(1,"All");
258 fHCuts->GetXaxis()->SetBinLabel(2,"PS");
259 fHCuts->GetXaxis()->SetBinLabel(3,Form("%s: %.0f-%.0f",fCentVar.Data(),fCentFrom,fCentTo));
260 fHCuts->GetXaxis()->SetBinLabel(4,"QFlag");
261 fHCuts->GetXaxis()->SetBinLabel(5,Form("zvtx: %.0f-%.0f",fVtxZMin,fVtxZMax));
fa443410 262 fOutput->Add(fHCuts);
263 fHVertexZ = new TH1F("hVertexZBeforeCut","",100,-25,25);
264 fHVertexZ->SetXTitle("z [cm]");
265 fOutput->Add(fHVertexZ);
76332037 266 fHVertexZ2 = new TH1F("hVertexZAfterCut","",100,-25,25);
267 fHVertexZ2->SetXTitle("z [cm]");
268 fOutput->Add(fHVertexZ2);
f2b8fca6 269 fHCent = new TH1F("hCentBeforeCut","",102,-1,101);
fa443410 270 fHCent->SetXTitle(fCentVar.Data());
76332037 271 fOutput->Add(fHCent);
f2b8fca6 272 fHCentQual = new TH1F("hCentAfterCut","",102,-1,101);
fa443410 273 fHCentQual->SetXTitle(fCentVar.Data());
274 fOutput->Add(fHCentQual);
2e4d8148 275 fHTclsBeforeCuts = new TH1F("hTclsBeforeCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
276 fHTclsAfterCuts = new TH1F("hTclsAfterCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
b3ee6797 277 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
278 const char *name = fTrClassNamesArr->At(i)->GetName();
279 fHTclsBeforeCuts->GetXaxis()->SetBinLabel(1+i,name);
280 fHTclsAfterCuts->GetXaxis()->SetBinLabel(1+i,name);
281 }
282 fOutput->Add(fHTclsBeforeCuts);
283 fOutput->Add(fHTclsAfterCuts);
90d5b88b 284
d595acbb 285 // histograms for cells
296ea9b4 286 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
287 fHColuRow = new TH2*[nsm];
288 fHColuRowE = new TH2*[nsm];
289 fHCellMult = new TH1*[nsm];
290 for (Int_t i = 0; i<nsm; ++i) {
2e4d8148 291 fHColuRow[i] = new TH2F(Form("hColRow_Mod%d",i),"",48,0,48,24,0.,24);
d595acbb 292 fHColuRow[i]->SetTitle(Form("Module %d: Occupancy", i));
293 fHColuRow[i]->SetXTitle("col (i#eta)");
294 fHColuRow[i]->SetYTitle("row (i#phi)");
2e4d8148 295 fHColuRowE[i] = new TH2F(Form("hColRowE_Mod%d", i),"",48,0,48,24,0,24);
90d5b88b 296 fHColuRowE[i]->SetTitle(Form("Module %d: Cell energy",i));
d595acbb 297 fHColuRowE[i]->SetXTitle("col (i#eta)");
298 fHColuRowE[i]->SetYTitle("row (i#phi)");
299 fHCellMult[i] = new TH1F(Form("hCellMult_Mod%d",i),"",1000,0,1000);
300 fHCellMult[i]->SetTitle(Form("Module %d: Cell multiplicity",i));
90d5b88b 301 fHCellMult[i]->SetXTitle("# of cells");
d595acbb 302 fOutput->Add(fHColuRow[i]);
303 fOutput->Add(fHColuRowE[i]);
304 fOutput->Add(fHCellMult[i]);
305 }
2e4d8148 306 fHCellE = new TH1F("hCellE","",250,0.,25.);
d595acbb 307 fHCellE->SetXTitle("E_{cell} [GeV]");
308 fOutput->Add(fHCellE);
fb5a0b69 309 fHCellH = new TH1F ("hCellHighestE","",250,0.,25.);
6eb6260e 310 fHCellH->SetXTitle("E^{max}_{cell} [GeV]");
d595acbb 311 fOutput->Add(fHCellH);
fb5a0b69 312 fHCellM = new TH1F ("hCellMeanEperHitCell","",250,0.,2.5);
6eb6260e 313 fHCellM->SetXTitle("#LT E_{cell}#GT [GeV]");
314 fOutput->Add(fHCellM);
fb5a0b69 315 fHCellM2 = new TH1F ("hCellMeanEperAllCells","",250,0.,1);
f4ec884e 316 fHCellM2->SetXTitle("1/N_{cells} #Sigma E_{cell} [GeV]");
a49742b5 317 fOutput->Add(fHCellM2);
90d5b88b 318
2e4d8148 319 fHCellFreqNoCut = new TH1*[nsm];
320 fHCellFreqCut100M = new TH1*[nsm];
321 fHCellFreqCut300M = new TH1*[nsm];
322 fHCellFreqE = new TH1*[nsm];
296ea9b4 323 for (Int_t i = 0; i<nsm; ++i){
324 Double_t lbin = i*24*48-0.5;
325 Double_t hbin = lbin+24*48;
2e4d8148 326 fHCellFreqNoCut[i] = new TH1F(Form("hCellFreqNoCut_SM%d",i),
327 Form("Frequency SM%d (no cut);id;#",i), 1152, lbin, hbin);
328 fHCellFreqCut100M[i] = new TH1F(Form("hCellFreqCut100M_SM%d",i),
329 Form("Frequency SM%d (>0.1GeV);id;#",i), 1152, lbin, hbin);
330 fHCellFreqCut300M[i] = new TH1F(Form("hCellFreqCut300M_SM%d",i),
331 Form("Frequency SM%d (>0.3GeV);id;#",i), 1152, lbin, hbin);
332 fHCellFreqE[i] = new TH1F(Form("hCellFreqE_SM%d",i),
333 Form("Frequency SM%d (E weighted);id;#",i), 1152, lbin, hbin);
296ea9b4 334 fOutput->Add(fHCellFreqNoCut[i]);
2e4d8148 335 fOutput->Add(fHCellFreqCut100M[i]);
336 fOutput->Add(fHCellFreqCut300M[i]);
337 fOutput->Add(fHCellFreqE[i]);
296ea9b4 338 }
3a952328 339 if (!fMarkCells.IsNull()) {
296ea9b4 340 fHCellCheckE = new TH1*[24*48*nsm];
408dc04c 341 memset(fHCellCheckE,0,24*48*nsm*sizeof(TH1*));
3a952328 342 TObjArray *cells = fMarkCells.Tokenize(" ");
343 Int_t n = cells->GetEntries();
344 Int_t *tcs = new Int_t[n];
345 for (Int_t i=0;i<n;++i) {
346 TString name(cells->At(i)->GetName());
347 tcs[i]=name.Atoi();
348 }
349 for (Int_t i = 0; i<n; ++i) {
296ea9b4 350 Int_t c=tcs[i];
351 if (c<24*48*nsm) {
3a952328 352 fHCellCheckE[i] = new TH1F(Form("hCellE_id%d",c), Form("Cell %d;E [GeV/c];#",c), 1000, 0, 10);
296ea9b4 353 fOutput->Add(fHCellCheckE[i]);
354 }
355 }
3a952328 356 delete cells;
357 delete [] tcs;
296ea9b4 358 }
359
d595acbb 360 // histograms for clusters
3a952328 361 if (!fTrainMode) {
362 fHClustEccentricity = new TH1F("hClustEccentricity","",100,-0.1,1.1);
363 fHClustEccentricity->SetXTitle("#epsilon_{C}");
364 fOutput->Add(fHClustEccentricity);
365 fHClustEtaPhi = new TH2F("hClustEtaPhi","",500,-0.8,0.8,100*nsm,phimin,phimax);
366 fHClustEtaPhi->SetXTitle("#eta");
367 fHClustEtaPhi->SetYTitle("#varphi");
368 fOutput->Add(fHClustEtaPhi);
369 fHClustEnergyPt = new TH2F("hClustEnergyPt","",250,0,50,250,0,50);
370 fHClustEnergyPt->SetXTitle("E [GeV]");
371 fHClustEnergyPt->SetYTitle("p_{T} [GeV/c]");
372 fOutput->Add(fHClustEnergyPt);
373 fHClustEnergySigma = new TH2F("hClustEnergySigma","",100,0,100,500,0,50);
374 fHClustEnergySigma->SetXTitle("E_{C} * #sigma_{max} [GeV*cm]");
375 fHClustEnergySigma->SetYTitle("E_{C} [GeV]");
376 fOutput->Add(fHClustEnergySigma);
377 fHClustSigmaSigma = new TH2F("hClustSigmaSigma","",500,0,50,500,0,50);
378 fHClustSigmaSigma->SetXTitle("#lambda_{0} [cm]");
379 fHClustSigmaSigma->SetYTitle("#sigma_{max} [cm]");
380 fOutput->Add(fHClustSigmaSigma);
381 fHClustNCellEnergyRatio = new TH2F("hClustNCellEnergyRatio","",27,-0.5,26.5,101,-0.05,1.05);
382 fHClustNCellEnergyRatio->SetXTitle("N_{cells}");
383 fHClustNCellEnergyRatio->SetYTitle("E^{max}_{cell}/E_{clus}");
384 fOutput->Add(fHClustNCellEnergyRatio);
385 }
90d5b88b 386
b6c599fe 387 // histograms for track matching
388 fHMatchDr = new TH1F("hMatchDrDist",";dR [cm]",500,0,200);
389 fOutput->Add(fHMatchDr);
390 fHMatchDz = new TH1F("hMatchDzDist",";dZ [cm]",500,-100,100);
391 fOutput->Add(fHMatchDz);
392 fHMatchEp = new TH1F("hMatchEpDist",";E/p",100,0,10);
393 fOutput->Add(fHMatchEp);
3a952328 394
d595acbb 395 // histograms for pion candidates
3a952328 396 if (!fTrainMode) {
397 fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100*nsm,phimin,phimax);
398 fHPionEtaPhi->SetXTitle("#eta_{#gamma#gamma}");
399 fHPionEtaPhi->SetYTitle("#varphi_{#gamma#gamma}");
400 fOutput->Add(fHPionEtaPhi);
401 fHPionMggPt = new TH2F("hPionMggPt","",1000,0,2,100,0,20.0);
402 fHPionMggPt->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
403 fHPionMggPt->SetYTitle("p_{T}^{#gamma#gamma} [GeV/c]");
404 fOutput->Add(fHPionMggPt);
405 fHPionMggAsym = new TH2F("hPionMggAsym","",1000,0,2,100,0,1);
406 fHPionMggAsym->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
407 fHPionMggAsym->SetYTitle("Z_{#gamma#gamma} [GeV]");
408 fOutput->Add(fHPionMggAsym);
409 fHPionMggDgg = new TH2F("hPionMggDgg","",1000,0,2,100,0,10);
410 fHPionMggDgg->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
411 fHPionMggDgg->SetYTitle("opening angle [grad]");
412 fOutput->Add(fHPionMggDgg);
413 const Int_t nbins = 20;
414 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};
415 fPtRanges = new TAxis(nbins-1,xbins);
416 for (Int_t i = 0; i<=nbins; ++i) {
417 fHPionInvMasses[i] = new TH1F(Form("hPionInvMass%d",i),"",1000,0,2);
418 fHPionInvMasses[i]->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
419 if (i==0)
420 fHPionInvMasses[i]->SetTitle(Form("0 < p_{T}^{#gamma#gamma} <%.1f",xbins[0]));
421 else if (i==nbins)
422 fHPionInvMasses[i]->SetTitle(Form("p_{T}^{#gamma#gamma} > 50"));
423 else
424 fHPionInvMasses[i]->SetTitle(Form("%.1f < p_{T}^{#gamma#gamma} <%.1f",xbins[i-1],xbins[i]));
425 fOutput->Add(fHPionInvMasses[i]);
426 }
fa443410 427 }
296ea9b4 428
717fe7de 429 PostData(1, fOutput);
ea3fd2d5 430}
431
432//________________________________________________________________________
433void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *)
434{
717fe7de 435 // Called for each event.
436
437 if (!InputEvent())
438 return;
439
43cfaa06 440 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
441 fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent());
44310bce 442 UInt_t offtrigger = 0;
43cfaa06 443 if (fEsdEv) {
444 am->LoadBranch("AliESDRun.");
445 am->LoadBranch("AliESDHeader.");
9809ed8c 446 UInt_t mask1 = fEsdEv->GetESDRun()->GetDetectorsInDAQ();
447 UInt_t mask2 = fEsdEv->GetESDRun()->GetDetectorsInReco();
448 Bool_t desc1 = (mask1 >> 18) & 0x1;
449 Bool_t desc2 = (mask2 >> 18) & 0x1;
450 if (desc1==0 || desc2==0) { //AliDAQ::OfflineModuleName(180=="EMCAL"
451 AliError(Form("EMCAL not in DAQ/RECO: %u (%u)/%u (%u)",
452 mask1, fEsdEv->GetESDRun()->GetDetectorsInReco(),
453 mask2, fEsdEv->GetESDRun()->GetDetectorsInDAQ()));
44310bce 454 return;
455 }
456 offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
43cfaa06 457 } else {
458 fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
408dc04c 459 if (!fAodEv) {
460 AliFatal("Neither ESD nor AOD event found");
461 return;
462 }
43cfaa06 463 am->LoadBranch("header");
44310bce 464 offtrigger = fAodEv->GetHeader()->GetOfflineTrigger();
465 }
38727e64 466 if (!fMcMode && (offtrigger & AliVEvent::kFastOnly)) {
44310bce 467 AliWarning(Form("EMCAL not in fast only partition"));
468 return;
43cfaa06 469 }
b6c599fe 470 if (fDoTrMatGeom && !AliGeomManager::GetGeometry()) { // get geometry
27c2e3d9 471 AliWarning("Accessing geometry from OCDB, this is not very efficient!");
472 AliCDBManager *cdb = AliCDBManager::Instance();
473 if (!cdb->IsDefaultStorageSet())
474 cdb->SetDefaultStorage("raw://");
475 Int_t runno = InputEvent()->GetRunNumber();
476 if (runno != cdb->GetRun())
477 cdb->SetRun(runno);
478 AliGeomManager::LoadGeometry();
479 }
480
481 if (!AliGeomManager::GetGeometry()&&!fIsGeoMatsSet) { // set misalignment matrices (stored in first event)
482 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
9809ed8c 483 for (Int_t i=0; i<nsm; ++i) {
484 const TGeoHMatrix *geom = 0;
485 if (fEsdEv)
486 geom = fEsdEv->GetESDRun()->GetEMCALMatrix(i);
487 else
488 geom = fAodEv->GetHeader()->GetEMCALMatrix(i);
489 if (!geom)
490 continue;
491 geom->Print();
492 fGeom->SetMisalMatrix(geom,i);
27c2e3d9 493 }
494 fIsGeoMatsSet = kTRUE;
495 }
496
497 if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
498 if (fEsdEv) {
499 const AliESDRun *erun = fEsdEv->GetESDRun();
500 AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(),
501 erun->GetCurrentDip(),
502 AliMagF::kConvLHC,
503 kFALSE,
504 erun->GetBeamEnergy(),
505 erun->GetBeamType());
506 TGeoGlobalMagField::Instance()->SetField(field);
507 } else {
508 Double_t pol = -1; //polarity
509 Double_t be = -1; //beam energy
510 AliMagF::BeamType_t btype = AliMagF::kBeamTypepp;
511 Int_t runno = fAodEv->GetRunNumber();
512 if (runno>=136851 && runno<138275) {
513 pol = -1;
514 be = 2760;
515 btype = AliMagF::kBeamTypeAA;
516 } else if (runno>=138275 && runno<=139517) {
517 pol = +1;
518 be = 2760;
519 btype = AliMagF::kBeamTypeAA;
520 } else {
521 AliError(Form("Do not know the bfield parameters for run %d! Using defaults!!!", runno));
522 }
523 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", pol, pol, AliMagF::k5kG, btype, be));
524 }
525 }
526
b3ee6797 527 Int_t cut = 1;
528 fHCuts->Fill(cut++);
529
530 TString trgclasses;
531 AliESDHeader *h = dynamic_cast<AliESDHeader*>(InputEvent()->GetHeader());
532 if (h) {
533 trgclasses = fEsdEv->GetFiredTriggerClasses();
534 } else {
535 AliAODHeader *h2 = dynamic_cast<AliAODHeader*>(InputEvent()->GetHeader());
536 if (h2)
537 trgclasses = h2->GetFiredTriggerClasses();
538 }
539 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
540 const char *name = fTrClassNamesArr->At(i)->GetName();
541 if (trgclasses.Contains(name))
542 fHTclsBeforeCuts->Fill(1+i);
543 }
544
545 UInt_t res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
546 if (res==0)
547 return;
548
fa443410 549 fHCuts->Fill(cut++);
286b47a5 550
717fe7de 551 const AliCentrality *centP = InputEvent()->GetCentrality();
552 Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar);
fa443410 553 fHCent->Fill(cent);
717fe7de 554 if (cent<fCentFrom||cent>fCentTo)
286b47a5 555 return;
43cfaa06 556
fa443410 557 fHCuts->Fill(cut++);
286b47a5 558
43cfaa06 559 if (fUseQualFlag) {
560 if (centP->GetQuality()>0)
561 return;
562 }
286b47a5 563
fa443410 564 fHCentQual->Fill(cent);
565 fHCuts->Fill(cut++);
717fe7de 566
43cfaa06 567 if (fEsdEv) {
fa443410 568 am->LoadBranch("PrimaryVertex.");
569 am->LoadBranch("SPDVertex.");
570 am->LoadBranch("TPCVertex.");
43cfaa06 571 } else {
572 fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
573 am->LoadBranch("vertices");
3c661da5 574 if (!fAodEv) return;
43cfaa06 575 }
576
577 const AliVVertex *vertex = InputEvent()->GetPrimaryVertex();
717fe7de 578 if (!vertex)
579 return;
580
fa443410 581 fHVertexZ->Fill(vertex->GetZ());
286b47a5 582
717fe7de 583 if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax)
584 return;
286b47a5 585
fa443410 586 fHCuts->Fill(cut++);
76332037 587 fHVertexZ2->Fill(vertex->GetZ());
717fe7de 588
b3ee6797 589 // count number of accepted events
590 ++fNEvs;
591
592 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
593 const char *name = fTrClassNamesArr->At(i)->GetName();
594 if (trgclasses.Contains(name))
595 fHTclsAfterCuts->Fill(1+i);
596 }
597
43cfaa06 598 fRecPoints = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
b6c599fe 599 fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set or if clusters are attached
43cfaa06 600 fEsdCells = 0; // will be set if ESD input used
b6c599fe 601 fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set or if clusters are attached
43cfaa06 602 // or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode
603 fAodCells = 0; // will be set if AOD input used
717fe7de 604
43cfaa06 605 // deal with special output from AliAnalysisTaskEMCALClusterizeFast first
606 Bool_t clusattached = 0;
607 Bool_t recalibrated = 0;
608 if (1 && !fClusName.IsNull()) {
609 AliAnalysisTaskEMCALClusterizeFast *cltask = 0;
610 TObjArray *ts = am->GetTasks();
611 cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName));
612 if (cltask && cltask->GetClusters()) {
d595acbb 613 fRecPoints = const_cast<TObjArray*>(cltask->GetClusters());
43cfaa06 614 clusattached = cltask->GetAttachClusters();
615 if (cltask->GetCalibData()!=0)
616 recalibrated = kTRUE;
617 }
618 }
38727e64 619 if (1 && !fClusName.IsNull()) {
620 TList *l = 0;
621 if (AODEvent())
622 l = AODEvent()->GetList();
623 else if (fAodEv)
624 l = fAodEv->GetList();
717fe7de 625 if (l) {
38727e64 626 fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
717fe7de 627 }
ea3fd2d5 628 }
717fe7de 629
630 if (fEsdEv) { // ESD input mode
43cfaa06 631 if (1 && (!fRecPoints||clusattached)) {
632 if (!clusattached)
633 am->LoadBranch("CaloClusters");
717fe7de 634 TList *l = fEsdEv->GetList();
717fe7de 635 if (l) {
38727e64 636 fEsdClusters = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
ea3fd2d5 637 }
638 }
43cfaa06 639 if (1) {
640 if (!recalibrated)
641 am->LoadBranch("EMCALCells.");
717fe7de 642 fEsdCells = fEsdEv->GetEMCALCells();
643 }
644 } else if (fAodEv) { // AOD input mode
43cfaa06 645 if (1 && (!fAodClusters || clusattached)) {
646 if (!clusattached)
647 am->LoadBranch("caloClusters");
717fe7de 648 TList *l = fAodEv->GetList();
717fe7de 649 if (l) {
38727e64 650 fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
ea3fd2d5 651 }
652 }
43cfaa06 653 if (1) {
654 if (!recalibrated)
655 am->LoadBranch("emcalCells");
717fe7de 656 fAodCells = fAodEv->GetEMCALCells();
657 }
658 } else {
659 AliFatal("Impossible to not have either pointer to ESD or AOD event");
660 }
ea3fd2d5 661
43cfaa06 662 if (1) {
663 AliDebug(2,Form("fRecPoints set: %p", fRecPoints));
664 AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters));
665 AliDebug(2,Form("fEsdCells set: %p", fEsdCells));
666 AliDebug(2,Form("fAodClusters set: %p", fAodClusters));
667 AliDebug(2,Form("fAodCells set: %p", fAodCells));
668 }
669
a49742b5 670 if (fDoAfterburner)
671 ClusterAfterburner();
6eb6260e 672
38727e64 673 CalcMcInfo();
3a952328 674 CalcCaloTriggers();
b6c599fe 675 CalcPrimTracks();
3a952328 676 CalcTracks();
296ea9b4 677 CalcClusterProps();
678
286b47a5 679 FillCellHists();
3a952328 680 if (!fTrainMode) {
681 FillClusHists();
682 FillPionHists();
683 FillOtherHists();
684 }
38727e64 685 FillMcHists();
788ca675 686 FillNtuple();
ea3fd2d5 687
3a952328 688 if (fTrainMode) {
689 fSelTracks->Clear();
690 fSelPrimTracks->Clear();
691 fTriggers->Clear();
692 }
de4cee41 693
717fe7de 694 PostData(1, fOutput);
ea3fd2d5 695}
696
697//________________________________________________________________________
698void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *)
699{
717fe7de 700 // Terminate called at the end of analysis.
f5d4ab70 701
702 if (fNtuple) {
703 TFile *f = OpenFile(1);
704 if (f)
705 fNtuple->Write();
706 }
d9f26424 707
b3ee6797 708 AliInfo(Form("%s: Accepted %lld events", GetName(), fNEvs));
ea3fd2d5 709}
ea3fd2d5 710
296ea9b4 711//________________________________________________________________________
3a952328 712void AliAnalysisTaskEMCALPi0PbPb::CalcCaloTriggers()
296ea9b4 713{
3a952328 714 // Calculate triggers
296ea9b4 715
38727e64 716 if (fAodEv)
717 return; // information not available in AOD
718
3a952328 719 fTriggers->Clear();
296ea9b4 720
3a952328 721 AliVCaloCells *cells = fEsdCells;
722 if (!cells)
723 cells = fAodCells;
724 if (!cells)
725 return;
726
727 Int_t ncells = cells->GetNumberOfCells();
728 if (ncells<=0)
729 return;
730
731 if (ncells>fNAmpInTrigger) {
732 delete [] fAmpInTrigger;
733 fAmpInTrigger = new Float_t[ncells];
734 fNAmpInTrigger = ncells;
296ea9b4 735 }
3a952328 736 for (Int_t i=0;i<ncells;++i)
737 fAmpInTrigger[i] = 0;
296ea9b4 738
3a952328 739 std::map<Short_t,Short_t> map;
740 for (Short_t pos=0;pos<ncells;++pos) {
741 Short_t id = cells->GetCellNumber(pos);
742 map[id]=pos;
743 }
744
745 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
746 am->LoadBranch("EMCALTrigger.");
747
748 AliESDCaloTrigger *triggers = fEsdEv->GetCaloTrigger("EMCAL");
749 if (!triggers)
750 return;
751 if (triggers->GetEntries()<=0)
b6c599fe 752 return;
753
3a952328 754 triggers->Reset();
755 Int_t ntrigs=0;
756 while (triggers->Next()) {
757 Int_t gCol=0, gRow=0, ntimes=0;
758 triggers->GetPosition(gCol,gRow);
759 triggers->GetNL0Times(ntimes);
760 if (ntimes<1)
761 continue;
762 Float_t amp=0;
763 triggers->GetAmplitude(amp);
764 Int_t find = -1;
765 fGeom->GetAbsFastORIndexFromPositionInEMCAL(gCol,gRow,find);
766 if (find<0)
767 continue;
768 Int_t cidx[4] = {-1};
769 Bool_t ret = fGeom->GetCellIndexFromFastORIndex(find, cidx);
770 if (!ret)
771 continue;
772 Int_t trgtimes[25];
773 triggers->GetL0Times(trgtimes);
774 Int_t mintime = trgtimes[0];
775 Int_t maxtime = trgtimes[0];
776 Bool_t trigInTimeWindow = 0;
777 for (Int_t i=0;i<ntimes;++i) {
778 if (trgtimes[i]<mintime)
779 mintime = trgtimes[i];
780 if (maxtime<trgtimes[i])
781 maxtime = trgtimes[i];
f2b8fca6 782 if ((fMinL0Time<=trgtimes[i]) && (fMaxL0Time>=trgtimes[i]))
3a952328 783 trigInTimeWindow = 1;
784 }
785
786 Double_t tenergy = 0;
787 Double_t tphi=0;
788 Double_t teta=0;
789 for (Int_t i=0;i<3;++i) {
790 Short_t pos = -1;
791 std::map<Short_t,Short_t>::iterator it = map.find(cidx[i]);
792 if (it!=map.end())
793 pos = it->second;
794 if (pos<0)
b6c599fe 795 continue;
3a952328 796 if (trigInTimeWindow)
5fe1ca23 797 fAmpInTrigger[pos] = amp;
3a952328 798 Float_t eta=-1, phi=-1;
799 fGeom->EtaPhiFromIndex(cidx[i],eta,phi);
800 Double_t en= cells->GetAmplitude(pos);
801 tenergy+=en;
802 teta+=eta*en;
803 tphi+=phi*en;
804 }
805 teta/=tenergy;
806 tphi/=tenergy;
807 AliStaTrigger *trignew = static_cast<AliStaTrigger*>(fTriggers->New(ntrigs++));
808 trignew->fE = tenergy;
809 trignew->fEta = teta;
810 trignew->fPhi = tphi;
811 trignew->fAmp = amp;
812 trignew->fMinTime = mintime;
813 trignew->fMaxTime = maxtime;
814 }
815}
816
817//________________________________________________________________________
818void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
819{
820 // Calculate cluster properties
821
822 fClusters->Clear();
823
824 AliVCaloCells *cells = fEsdCells;
825 if (!cells)
826 cells = fAodCells;
827 if (!cells)
828 return;
829
830 TObjArray *clusters = fEsdClusters;
831 if (!clusters)
832 clusters = fAodClusters;
833 if (!clusters)
834 return;
835
836 Int_t ncells = cells->GetNumberOfCells();
837 Int_t nclus = clusters->GetEntries();
838 Int_t ntrks = fSelTracks->GetEntries();
839 Bool_t btracks[6][ntrks];
840 memset(btracks,0,sizeof(btracks));
841
842 std::map<Short_t,Short_t> map;
843 for (Short_t pos=0;pos<ncells;++pos) {
844 Short_t id = cells->GetCellNumber(pos);
845 map[id]=pos;
846 }
847
848 for(Int_t i=0, ncl=0; i<nclus; ++i) {
849 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
850 if (!clus)
851 continue;
852 if (!clus->IsEMCAL())
853 continue;
854 if (clus->E()<fMinE)
855 continue;
856
857 Float_t clsPos[3] = {0};
858 clus->GetPosition(clsPos);
859 TVector3 clsVec(clsPos);
860 Double_t vertex[3] = {0};
861 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
862 TLorentzVector clusterVec;
863 clus->GetMomentum(clusterVec,vertex);
864 Double_t clsEta = clusterVec.Eta();
865
866 AliStaCluster *cl = static_cast<AliStaCluster*>(fClusters->New(ncl++));
867 cl->fE = clus->E();
868 cl->fR = clsVec.Perp();
869 cl->fEta = clsVec.Eta();
870 cl->fPhi = clsVec.Phi();
871 cl->fN = clus->GetNCells();
872 cl->fN1 = GetNCells(clus,0.100);
873 cl->fN3 = GetNCells(clus,0.300);
874 Short_t id = -1;
875 Double_t emax = GetMaxCellEnergy(clus, id);
876 cl->fIdMax = id;
877 cl->fEmax = emax;
878 if (clus->GetDistanceToBadChannel()<10000)
879 cl->fDbc = clus->GetDistanceToBadChannel();
880 if (!TMath::IsNaN(clus->GetDispersion()))
881 cl->fDisp = clus->GetDispersion();
882 if (!TMath::IsNaN(clus->GetM20()))
0fbe8d4f 883 cl->fM20 = clus->GetM20();
3a952328 884 if (!TMath::IsNaN(clus->GetM02()))
0fbe8d4f 885 cl->fM02 = clus->GetM02();
3a952328 886 Double_t maxAxis = 0, minAxis = 0;
887 GetSigma(clus,maxAxis,minAxis);
888 clus->SetTOF(maxAxis); // store sigma in TOF
889 cl->fSig = maxAxis;
890 Double_t clusterEcc = 0;
891 if (maxAxis > 0)
892 clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
893 clus->SetChi2(clusterEcc); // store ecc in chi2
894 cl->fEcc = clusterEcc;
895 cl->fTrIso = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
896 cl->fTrIso1 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 1);
897 cl->fTrIso2 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 2);
898 cl->fCeCore = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.05);
899 cl->fCeIso = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
0fbe8d4f 900
38727e64 901 if (fAmpInTrigger) { // fill trigger info if present
902 Double_t trigpen = 0;
903 Double_t trignen = 0;
904 for(Int_t j=0; j<cl->fN; ++j) {
905 Short_t cid = TMath::Abs(clus->GetCellAbsId(j));
906 Short_t pos = -1;
907 std::map<Short_t,Short_t>::iterator it = map.find(cid);
908 if (it!=map.end())
909 pos = it->second;
910 if (pos<0)
911 continue;
912 if (fAmpInTrigger[pos]>0)
913 trigpen += cells->GetAmplitude(pos);
914 else if (fAmpInTrigger[pos]<0)
915 trignen += cells->GetAmplitude(pos);
916 }
917 if (trigpen>0) {
918 cl->fIsTrigM = 1;
919 cl->fTrigE = trigpen;
920 }
921 if (trignen>0) {
922 cl->fIsTrigM = 1;
923 cl->fTrigMaskE = trignen;
924 }
3a952328 925 }
0fbe8d4f 926 cl->fIsShared = IsShared(clus);
3a952328 927
928 // track matching
929 Double_t mind2 = 1e10;
930 for(Int_t j = 0; j<ntrks; ++j) {
931 AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
b6c599fe 932 if (!track)
933 continue;
3a952328 934
935 if (TMath::Abs(clsEta-track->Eta())>0.5)
b6c599fe 936 continue;
3a952328 937
938 TVector3 vec(clsPos);
939 Int_t index = (Int_t)(vec.Phi()*TMath::RadToDeg()/20);
940 if (btracks[index-4][j]) {
b6c599fe 941 continue;
3a952328 942 }
b6c599fe 943
3a952328 944 Float_t tmpR=-1, tmpZ=-1;
945 if (!fDoTrMatGeom) {
946 AliExternalTrackParam *tParam = 0;
947 if (fEsdEv) {
948 AliESDtrack *esdTrack = static_cast<AliESDtrack*>(track);
949 tParam = new AliExternalTrackParam(*esdTrack->GetTPCInnerParam());
950 } else
951 tParam = new AliExternalTrackParam(track);
952
953 Double_t bfield[3] = {0};
954 track->GetBxByBz(bfield);
955 Double_t alpha = (index+0.5)*20*TMath::DegToRad();
956 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
957 tParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
958 Bool_t ret = tParam->PropagateToBxByBz(vec.X(), bfield);
959 if (!ret) {
960 btracks[index-4][j]=1;
961 delete tParam;
962 continue;
b6c599fe 963 }
3a952328 964 Double_t trkPos[3] = {0};
965 tParam->GetXYZ(trkPos); //Get the extrapolated global position
38727e64 966 tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2) +
967 TMath::Power(clsPos[1]-trkPos[1],2) +
968 TMath::Power(clsPos[2]-trkPos[2],2) );
3a952328 969 tmpZ = clsPos[2]-trkPos[2];
970 delete tParam;
971 } else {
972 if (TMath::Abs(clsEta-track->Eta())>fIsoDist)
973 continue;
974 AliExternalTrackParam tParam(track);
975 if (!fReco->ExtrapolateTrackToCluster(&tParam, clus, tmpR, tmpZ))
976 continue;
b6c599fe 977 }
3a952328 978
979 Double_t d2 = tmpR;
980 if (mind2>d2) {
981 mind2=d2;
982 cl->fTrDz = tmpZ;
983 cl->fTrDr = TMath::Sqrt(tmpR*tmpR-tmpZ*tmpZ);
984 cl->fTrEp = clus->E()/track->P();
f3582e89 985 cl->fIsTrackM = 1;
3a952328 986 }
987 }
988
f3582e89 989 if (cl->fIsTrackM) {
3a952328 990 fHMatchDr->Fill(cl->fTrDr);
991 fHMatchDz->Fill(cl->fTrDz);
992 fHMatchEp->Fill(cl->fTrEp);
b6c599fe 993 }
994 }
995}
996
997//________________________________________________________________________
998void AliAnalysisTaskEMCALPi0PbPb::CalcPrimTracks()
999{
1000 // Calculate track properties.
1001
1002 fSelPrimTracks->Clear();
1003
1004 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1005 TClonesArray *tracks = 0;
1006 if (fEsdEv) {
1007 am->LoadBranch("Tracks");
1008 TList *l = fEsdEv->GetList();
1009 tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
1010 } else {
1011 am->LoadBranch("tracks");
1012 TList *l = fAodEv->GetList();
1013 tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1014 }
1015
296ea9b4 1016 if (!tracks)
1017 return;
1018
0ec74551 1019 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
b6c599fe 1020 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-fIsoDist*1.25;
1021 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+fIsoDist*1.25;
0ec74551 1022
296ea9b4 1023 if (fEsdEv) {
b6c599fe 1024 fSelPrimTracks->SetOwner(kTRUE);
0ec74551 1025 am->LoadBranch("PrimaryVertex.");
1026 const AliESDVertex *vtx = fEsdEv->GetPrimaryVertexSPD();
1027 am->LoadBranch("SPDVertex.");
1028 const AliESDVertex *vtxSPD = fEsdEv->GetPrimaryVertexSPD();
1029 am->LoadBranch("Tracks");
1030 const Int_t Ntracks = tracks->GetEntries();
1031 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1032 AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1033 if (!track) {
1034 AliWarning(Form("Could not receive track %d\n", iTracks));
1035 continue;
1036 }
1037 if (fTrCuts && !fTrCuts->IsSelected(track))
296ea9b4 1038 continue;
1039 Double_t eta = track->Eta();
1040 if (eta<-1||eta>1)
1041 continue;
0ec74551 1042 if (track->Phi()<phimin||track->Phi()>phimax)
1043 continue;
2e4d8148 1044
0ec74551 1045 AliESDtrack copyt(*track);
1046 Double_t bfield[3];
1047 copyt.GetBxByBz(bfield);
1048 AliExternalTrackParam tParam;
1049 Bool_t relate = copyt.RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&tParam);
1050 if (!relate)
1051 continue;
1052 copyt.Set(tParam.GetX(),tParam.GetAlpha(),tParam.GetParameter(),tParam.GetCovariance());
2e4d8148 1053
0ec74551 1054 Double_t p[3] = { 0. };
1055 copyt.GetPxPyPz(p);
1056 Double_t pos[3] = { 0. };
1057 copyt.GetXYZ(pos);
1058 Double_t covTr[21] = { 0. };
1059 copyt.GetCovarianceXYZPxPyPz(covTr);
1060 Double_t pid[10] = { 0. };
1061 copyt.GetESDpid(pid);
1062 AliAODTrack *aTrack = new AliAODTrack(copyt.GetID(),
1063 copyt.GetLabel(),
1064 p,
1065 kTRUE,
1066 pos,
1067 kFALSE,
1068 covTr,
1069 (Short_t)copyt.GetSign(),
1070 copyt.GetITSClusterMap(),
1071 pid,
1072 0,/*fPrimaryVertex,*/
1073 kTRUE, // check if this is right
1074 vtx->UsesTrack(copyt.GetID()));
1075 aTrack->SetTPCClusterMap(copyt.GetTPCClusterMap());
1076 aTrack->SetTPCSharedMap (copyt.GetTPCSharedMap());
1077 Float_t ndf = copyt.GetTPCNcls() + 1 - 5 ;
1078 if(ndf>0)
1079 aTrack->SetChi2perNDF(copyt.GetTPCchi2()/ndf);
1080 else
1081 aTrack->SetChi2perNDF(-1);
1082 aTrack->SetFlags(copyt.GetStatus());
1083 aTrack->SetTPCPointsF(copyt.GetTPCNclsF());
b6c599fe 1084 fSelPrimTracks->Add(aTrack);
296ea9b4 1085 }
1086 } else {
1087 Int_t ntracks = tracks->GetEntries();
1088 for (Int_t i=0; i<ntracks; ++i) {
1089 AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
1090 if (!track)
1091 continue;
1092 Double_t eta = track->Eta();
1093 if (eta<-1||eta>1)
1094 continue;
0ec74551 1095 if (track->Phi()<phimin||track->Phi()>phimax)
1096 continue;
b6c599fe 1097 if(track->GetTPCNcls()<fMinNClusPerTr)
296ea9b4 1098 continue;
b6c599fe 1099 //todo: Learn how to set/filter AODs for prim/sec tracks
1100 fSelPrimTracks->Add(track);
296ea9b4 1101 }
1102 }
1103}
1104
1105//________________________________________________________________________
3a952328 1106void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
296ea9b4 1107{
3a952328 1108 // Calculate track properties (including secondaries).
b6c599fe 1109
3a952328 1110 fSelTracks->Clear();
296ea9b4 1111
3a952328 1112 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1113 TClonesArray *tracks = 0;
1114 if (fEsdEv) {
1115 am->LoadBranch("Tracks");
1116 TList *l = fEsdEv->GetList();
1117 tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
1118 } else {
1119 am->LoadBranch("tracks");
1120 TList *l = fAodEv->GetList();
1121 tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1122 }
296ea9b4 1123
3a952328 1124 if (!tracks)
1125 return;
296ea9b4 1126
3a952328 1127 if (fEsdEv) {
1128 const Int_t Ntracks = tracks->GetEntries();
1129 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1130 AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1131 if (!track) {
1132 AliWarning(Form("Could not receive track %d\n", iTracks));
1133 continue;
1134 }
1135 if (fTrCuts && !fTrCuts->IsSelected(track))
1136 continue;
1137 Double_t eta = track->Eta();
1138 if (eta<-1||eta>1)
1139 continue;
1140 fSelTracks->Add(track);
1141 }
1142 } else {
1143 Int_t ntracks = tracks->GetEntries();
1144 for (Int_t i=0; i<ntracks; ++i) {
1145 AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
296ea9b4 1146 if (!track)
1147 continue;
3a952328 1148 Double_t eta = track->Eta();
1149 if (eta<-1||eta>1)
c4236a58 1150 continue;
3a952328 1151 if(track->GetTPCNcls()<fMinNClusPerTr)
b6c599fe 1152 continue;
2e4d8148 1153
3a952328 1154 if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL
2e4d8148 1155 AliExternalTrackParam tParam(track);
3a952328 1156 if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 438, 0.139, 1, kTRUE)) {
1157 Double_t trkPos[3];
1158 tParam.GetXYZ(trkPos);
1159 track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]);
1160 }
b6c599fe 1161 }
3a952328 1162 fSelTracks->Add(track);
c4236a58 1163 }
296ea9b4 1164 }
1165}
1166
323834f0 1167//________________________________________________________________________
3a952328 1168void AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner()
323834f0 1169{
296ea9b4 1170 // Run custer reconstruction afterburner.
323834f0 1171
1172 AliVCaloCells *cells = fEsdCells;
1173 if (!cells)
1174 cells = fAodCells;
1175
1176 if (!cells)
1177 return;
1178
1179 Int_t ncells = cells->GetNumberOfCells();
1180 if (ncells<=0)
1181 return;
1182
1183 Double_t cellMeanE = 0, cellSigE = 0;
1184 for (Int_t i = 0; i<ncells; ++i) {
1185 Double_t cellE = cells->GetAmplitude(i);
1186 cellMeanE += cellE;
1187 cellSigE += cellE*cellE;
1188 }
1189 cellMeanE /= ncells;
1190 cellSigE /= ncells;
1191 cellSigE -= cellMeanE*cellMeanE;
1192 if (cellSigE<0)
1193 cellSigE = 0;
1194 cellSigE = TMath::Sqrt(cellSigE / ncells);
1195
1196 Double_t subE = cellMeanE - 7*cellSigE;
1197 if (subE<0)
1198 return;
1199
1200 for (Int_t i = 0; i<ncells; ++i) {
1201 Short_t id=-1;
1202 Double_t amp=0,time=0;
1203 if (!cells->GetCell(i, id, amp, time))
1204 continue;
1205 amp -= cellMeanE;
1206 if (amp<0.001)
1207 amp = 0;
1208 cells->SetCell(i, id, amp, time);
1209 }
1210
1211 TObjArray *clusters = fEsdClusters;
1212 if (!clusters)
1213 clusters = fAodClusters;
323834f0 1214 if (!clusters)
1215 return;
1216
1217 Int_t nclus = clusters->GetEntries();
1218 for (Int_t i = 0; i<nclus; ++i) {
1219 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1220 if (!clus->IsEMCAL())
1221 continue;
1222 Int_t nc = clus->GetNCells();
1223 Double_t clusE = 0;
1224 UShort_t ids[100] = {0};
1225 Double_t fra[100] = {0};
1226 for (Int_t j = 0; j<nc; ++j) {
1227 Short_t id = TMath::Abs(clus->GetCellAbsId(j));
1228 Double_t cen = cells->GetCellAmplitude(id);
1229 clusE += cen;
1230 if (cen>0) {
1231 ids[nc] = id;
1232 ++nc;
1233 }
1234 }
1235 if (clusE<=0) {
1236 clusters->RemoveAt(i);
1237 continue;
1238 }
1239
1240 for (Int_t j = 0; j<nc; ++j) {
1241 Short_t id = ids[j];
1242 Double_t cen = cells->GetCellAmplitude(id);
1243 fra[j] = cen/clusE;
1244 }
1245 clus->SetE(clusE);
1246 AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus);
1247 if (aodclus) {
1248 aodclus->Clear("");
1249 aodclus->SetNCells(nc);
1250 aodclus->SetCellsAmplitudeFraction(fra);
1251 aodclus->SetCellsAbsId(ids);
1252 }
1253 }
1254 clusters->Compress();
1255}
1256
286b47a5 1257//________________________________________________________________________
1258void AliAnalysisTaskEMCALPi0PbPb::FillCellHists()
1259{
1260 // Fill histograms related to cell properties.
d595acbb 1261
90d5b88b 1262 AliVCaloCells *cells = fEsdCells;
1263 if (!cells)
1264 cells = fAodCells;
1265
296ea9b4 1266 if (!cells)
1267 return;
90d5b88b 1268
296ea9b4 1269 Int_t cellModCount[12] = {0};
1270 Double_t cellMaxE = 0;
1271 Double_t cellMeanE = 0;
1272 Int_t ncells = cells->GetNumberOfCells();
1273 if (ncells==0)
1274 return;
1275
1276 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
1277
1278 for (Int_t i = 0; i<ncells; ++i) {
1279 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
1280 Double_t cellE = cells->GetAmplitude(i);
1281 fHCellE->Fill(cellE);
1282 if (cellE>cellMaxE)
1283 cellMaxE = cellE;
1284 cellMeanE += cellE;
1285
1286 Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1;
1287 Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta);
1288 if (!ret) {
1289 AliError(Form("Could not get cell index for %d", absID));
1290 continue;
1291 }
1292 ++cellModCount[iSM];
1293 Int_t iPhi=-1, iEta=-1;
1294 fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta);
1295 fHColuRow[iSM]->Fill(iEta,iPhi,1);
1296 fHColuRowE[iSM]->Fill(iEta,iPhi,cellE);
1297 fHCellFreqNoCut[iSM]->Fill(absID);
2e4d8148 1298 if (cellE > 0.1) fHCellFreqCut100M[iSM]->Fill(absID);
1299 if (cellE > 0.3) fHCellFreqCut300M[iSM]->Fill(absID);
296ea9b4 1300 if (fHCellCheckE && fHCellCheckE[absID])
1301 fHCellCheckE[absID]->Fill(cellE);
2e4d8148 1302 fHCellFreqE[iSM]->Fill(absID, cellE);
296ea9b4 1303 }
1304 fHCellH->Fill(cellMaxE);
1305 cellMeanE /= ncells;
1306 fHCellM->Fill(cellMeanE);
1307 fHCellM2->Fill(cellMeanE*ncells/24/48/nsm);
1308 for (Int_t i=0; i<nsm; ++i)
1309 fHCellMult[i]->Fill(cellModCount[i]);
286b47a5 1310}
1311
1312//________________________________________________________________________
1313void AliAnalysisTaskEMCALPi0PbPb::FillClusHists()
1314{
90d5b88b 1315 // Fill histograms related to cluster properties.
fa443410 1316
90d5b88b 1317 TObjArray *clusters = fEsdClusters;
1318 if (!clusters)
1319 clusters = fAodClusters;
296ea9b4 1320 if (!clusters)
1321 return;
90d5b88b 1322
296ea9b4 1323 Double_t vertex[3] = {0};
1324 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1325
1326 Int_t nclus = clusters->GetEntries();
1327 for(Int_t i = 0; i<nclus; ++i) {
1328 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1329 if (!clus)
1330 continue;
1331 if (!clus->IsEMCAL())
1332 continue;
90d5b88b 1333 TLorentzVector clusterVec;
296ea9b4 1334 clus->GetMomentum(clusterVec,vertex);
3a952328 1335 Double_t maxAxis = clus->GetTOF(); //sigma
1336 Double_t clusterEcc = clus->Chi2(); //eccentricity
296ea9b4 1337 fHClustEccentricity->Fill(clusterEcc);
1338 fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi());
1339 fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt());
1340 fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E());
1341 fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis);
1342 fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E());
788ca675 1343 }
1344}
b3ee6797 1345
38727e64 1346//________________________________________________________________________
1347void AliAnalysisTaskEMCALPi0PbPb::CalcMcInfo()
1348{
1349 // Get Mc truth particle information.
1350
1351 if (!fMcMode)
1352 return;
1353
1354// todo: treat mc differently for AOD and ESD. Check with Amorsch when we get MC handler also for AOD case.
1355
1356 AliMCEvent *mcEvent = MCEvent();
1357 if (!mcEvent)
1358 return;
1359
1360 const AliVVertex *evtVtx = mcEvent->GetPrimaryVertex();
1361 if (!evtVtx)
1362 return;
1363
1364 mcEvent->PreReadAll();
1365 AliStack *stack = mcEvent->Stack();
1366
1367 Int_t nTracks = mcEvent->GetNumberOfPrimaries();
1368 for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
1369 AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
1370 if (!mcP)
1371 continue;
1372
1373 // pion or eta meson
1374 TParticle *tP = mcP->Particle();
1375 if(tP->GetPdgCode() == 111) {
1376 } else if(tP->GetPdgCode() == 221) {
1377 } else
1378 continue;
1379
1380 // primary particle
1381 Double_t dR = TMath::Sqrt((tP->Vx()-evtVtx->GetX())*(tP->Vx()-evtVtx->GetX())+(tP->Vy()-evtVtx->GetY())*(tP->Vy()-evtVtx->GetY()));
1382 if(dR > 1.0)
1383 continue;
1384
1385 // do not account dalitz decays
1386 if(tP->GetNDaughters()!=2)
1387 continue;
1388
1389 // kinematic cuts
1390 Double_t pt = tP->Pt() ;
1391 if (pt<0.5) {
1392 continue;
1393 }
1394 Double_t rap = tP->Y();
1395 if (TMath::Abs(rap)>1){
1396 continue;
1397 }
1398
1399 // get daughters
1400 TParticle *tD1 = 0;
1401 TClonesArray *tD1refs = 0;
1402 mcEvent->GetParticleAndTR(tP->GetFirstDaughter(),tD1,tD1refs);
1403
1404 TParticle *tD2 = 0;
1405 TClonesArray *tD2refs = 0;
1406 mcEvent->GetParticleAndTR(tP->GetLastDaughter(),tD2,tD2refs);
1407
1408 TParticle *tD11 = 0;
1409 TClonesArray *tD11refs = 0;
1410 TParticle *tD12 = 0;
1411 TClonesArray *tD12refs = 0;
1412 Bool_t d1Dec = 0;
1413 if (tD1->GetNDaughters()>2) {
1414 AliWarning(Form("Decay photon has more than two daughters: %d", tD1->GetNDaughters()));
1415 return;
1416 }
1417 if (tD1->GetNDaughters()>=1)
1418 mcEvent->GetParticleAndTR(tD1->GetFirstDaughter(),tD11,tD11refs);
1419 if (tD1->GetNDaughters()==2) {
1420 mcEvent->GetParticleAndTR(tD1->GetLastDaughter(),tD12,tD12refs);
1421 d1Dec = 1;
1422 }
1423
1424 TParticle *tD21 = 0;
1425 TClonesArray *tD21refs = 0;
1426 TParticle *tD22 = 0;
1427 TClonesArray *tD22refs = 0;
1428 Bool_t d2Dec = 0;
1429 if (tD2->GetNDaughters()>2) {
1430 AliWarning(Form("Decay photon has more than two daughters: %d", tD2->GetNDaughters()));
1431 return;
1432 }
1433 if (tD2->GetNDaughters()>=1)
1434 mcEvent->GetParticleAndTR(tD2->GetFirstDaughter(),tD21,tD21refs);
1435 if (tD2->GetNDaughters()==2) {
1436 mcEvent->GetParticleAndTR(tD2->GetLastDaughter(),tD22,tD22refs);
1437 d2Dec = 1;
1438 }
1439
1440 Bool_t pDec = d1Dec && d2Dec;
1441 if (!pDec)
1442 continue;
1443
1444 if (tD11) {
1445 cout << " Dau 11 with " << tD11->GetNDaughters() << endl;
1446 tD11->Print();
1447 }
1448
1449 if (0) {
1450 tP->Print();
1451 tD1->Print();
1452 tD2->Print();
1453 if (tD11) {
1454 cout << " Dau 11 with " << tD11->GetNDaughters() << endl;
1455 tD11->Print();
1456 if (tD11refs) {
1457 Int_t nrefs = tD11refs->GetEntries();
1458 for (Int_t j=0;j<nrefs;++j) {
1459 AliTrackReference *tr = (AliTrackReference*)tD11refs->At(j);
1460 tr->SetUserId(tr->DetectorId());
1461 tr->Print();
1462 }
1463 }
1464 }
1465 if (tD12){
1466 cout << " Dau 12 with " << tD12->GetNDaughters() << endl;
1467 tD12->Print();
1468 if (tD12refs) {
1469 Int_t nrefs = tD12refs->GetEntries();
1470 for (Int_t j=0;j<nrefs;++j) {
1471 AliTrackReference *tr = (AliTrackReference*)tD12refs->At(j);
1472 tr->SetUserId(tr->DetectorId());
1473 tr->Print();
1474 }
1475 }
1476 }
1477 if (tD21) {
1478 cout << " Dau 21 with " << tD21->GetNDaughters() << endl;
1479 tD21->Print();
1480 if (tD21refs) {
1481 Int_t nrefs = tD21refs->GetEntries();
1482 for (Int_t j=0;j<nrefs;++j) {
1483 AliTrackReference *tr = (AliTrackReference*)tD21refs->At(j);
1484 tr->SetUserId(tr->DetectorId());
1485 tr->Print();
1486 }
1487 }
1488 }
1489 if (tD22) {
1490 cout << " Dau 22 with " << tD22->GetNDaughters() << endl;
1491 tD22->Print();
1492 if (tD22refs) {
1493 Int_t nrefs = tD22refs->GetEntries();
1494 for (Int_t j=0;j<nrefs;++j) {
1495 AliTrackReference *tr = (AliTrackReference*)tD22refs->At(j);
1496 tr->SetUserId(tr->DetectorId());
1497 tr->Print();
1498 }
1499 }
1500 }
1501 }
1502#if 0
1503 TParticle * gamma1 = fStack->Particle(particle->GetFirstDaughter());
1504 TParticle * gamma2 = fStack->Particle(particle->GetLastDaughter());
1505 //Number of pi0s decayed into acceptance
1506 Bool_t inAcc1 = (TMath::Abs(gamma1->Eta())<0.9) ;
1507 Bool_t inAcc2 = (TMath::Abs(gamma2->Eta())<0.9) ;
1508 Int_t mod1,mod2 ;
1509 Double_t x=0.,z=0. ;
1510 Bool_t hitPHOS1 = fPHOSgeom->ImpactOnEmc(gamma1, mod1, z,x) ;
1511 Bool_t hitPHOS2 = fPHOSgeom->ImpactOnEmc(gamma2, mod2, z,x) ;
1512 Bool_t hitEMCAL1= fEMCALgeom->Impact(gamma1) ;
1513 Bool_t hitEMCAL2= fEMCALgeom->Impact(gamma2) ;
1514
1515
1516 TParticle *mcP2 = 0;
1517 TClonesArray *refs = 0;
1518 mcEvent->GetParticleAndTR(iTrack,mcP2,refs);
1519 if (!refs)
1520 continue;
1521
1522 Int_t nrefs = refs->GetEntries();
1523 if (nrefs == 0)
1524 continue;
1525
1526 AliTrackReference *trEmcal = 0;
1527 for (Int_t j=0;j<nrefs;++j) {
1528 AliTrackReference *tr = (AliTrackReference*)refs->At(j);
1529 Int_t did = tr->DetectorId();
1530 if (did!=AliTrackReference::kEMCAL)
1531 continue;
1532 trEmcal = tr;
1533 break;
1534 }
1535
1536 if (!trEmcal)
1537 continue;
1538
1539 Double_t distA = TMath::Sqrt((trEmcal->X()-tP->Vx())*(trEmcal->X()-tP->Vx()) +
1540 (trEmcal->Y()-tP->Vy())*(trEmcal->Y()-tP->Vy()) +
1541 (trEmcal->Z()-tP->Vz())*(trEmcal->Z()-tP->Vz()));
1542 cout << "dist label to particle " << distA << " --------- " << endl;
1543 trEmcal->Print();
1544 tP->Print();
1545
1546 Int_t moi = tP->GetFirstMother();
1547 if (moi>=0) {
1548 AliMCParticle *mcMo = static_cast<AliMCParticle*>(mcEvent->GetTrack(moi));
1549 TParticle *tMo = mcMo->Particle();
1550 Double_t dist3d = TMath::Sqrt((tMo->Vx()-tP->Vx())*(tMo->Vx()-tP->Vx()) +
1551 (tMo->Vy()-tP->Vy())*(tMo->Vy()-tP->Vy()) +
1552 (tMo->Vz()-tP->Vz())*(tMo->Vz()-tP->Vz()));
1553 cout << "dist 3d " << dist3d << endl;
1554 tMo->Print();
1555 while (1) {
1556 Int_t newm = tMo->GetFirstMother();
1557 if (newm>=0) {
1558 AliMCParticle *newmo = static_cast<AliMCParticle*>(mcEvent->GetTrack(newm));
1559 tMo = newmo->Particle();
1560 tMo->Print();
1561 } else {
1562 break;
1563 }
1564 }
1565 cout << " ------- " << endl;
1566 }
1567#endif
1568 }
1569}
1570
788ca675 1571//________________________________________________________________________
1572void AliAnalysisTaskEMCALPi0PbPb::FillNtuple()
1573{
1574 // Fill ntuple.
1575
1576 if (!fNtuple)
1577 return;
1578
1579 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1580 if (fAodEv) {
1581 fHeader->fRun = fAodEv->GetRunNumber();
1582 fHeader->fOrbit = fAodEv->GetHeader()->GetOrbitNumber();
1583 fHeader->fPeriod = fAodEv->GetHeader()->GetPeriodNumber();
1584 fHeader->fBx = fAodEv->GetHeader()->GetBunchCrossNumber();
1585 fHeader->fL0 = fAodEv->GetHeader()->GetL0TriggerInputs();
1586 fHeader->fL1 = fAodEv->GetHeader()->GetL1TriggerInputs();
1587 fHeader->fL2 = fAodEv->GetHeader()->GetL2TriggerInputs();
1588 fHeader->fTrClassMask = fAodEv->GetHeader()->GetTriggerMask();
1589 fHeader->fTrCluster = fAodEv->GetHeader()->GetTriggerCluster();
1590 fHeader->fOffTriggers = fAodEv->GetHeader()->GetOfflineTrigger();
1591 fHeader->fFiredTriggers = fAodEv->GetHeader()->GetFiredTriggerClasses();
1592 } else {
1593 fHeader->fRun = fEsdEv->GetRunNumber();
1594 fHeader->fOrbit = fEsdEv->GetHeader()->GetOrbitNumber();
1595 fHeader->fPeriod = fEsdEv->GetESDRun()->GetPeriodNumber();
1596 fHeader->fBx = fEsdEv->GetHeader()->GetBunchCrossNumber();
1597 fHeader->fL0 = fEsdEv->GetHeader()->GetL0TriggerInputs();
1598 fHeader->fL1 = fEsdEv->GetHeader()->GetL1TriggerInputs();
1599 fHeader->fL2 = fEsdEv->GetHeader()->GetL2TriggerInputs();
1600 fHeader->fTrClassMask = fEsdEv->GetHeader()->GetTriggerMask();
1601 fHeader->fTrCluster = fEsdEv->GetHeader()->GetTriggerCluster();
1602 fHeader->fOffTriggers = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
1603 fHeader->fFiredTriggers = fEsdEv->GetFiredTriggerClasses();
5fe1ca23 1604 Float_t v0CorrR = 0;
1605 fHeader->fV0 = AliESDUtils::GetCorrV0(fEsdEv,v0CorrR);
1606 const AliMultiplicity *mult = fEsdEv->GetMultiplicity();
1607 if (mult)
1608 fHeader->fCl1 = mult->GetNumberOfITSClusters(1);
1609 fHeader->fTr = AliESDtrackCuts::GetReferenceMultiplicity(fEsdEv,1);
788ca675 1610 }
1611 AliCentrality *cent = InputEvent()->GetCentrality();
1612 fHeader->fV0Cent = cent->GetCentralityPercentileUnchecked("V0M");
1613 fHeader->fCl1Cent = cent->GetCentralityPercentileUnchecked("CL1");
1614 fHeader->fTrCent = cent->GetCentralityPercentileUnchecked("TRK");
1615 fHeader->fCqual = cent->GetQuality();
b6c599fe 1616
1617 AliEventplane *ep = InputEvent()->GetEventplane();
1618 if (ep) {
1619 if (ep->GetQVector())
1620 fHeader->fPsi = ep->GetQVector()->Phi()/2. ;
f2b8fca6 1621 else
1622 fHeader->fPsi = -1;
b6c599fe 1623 if (ep->GetQsub1()&&ep->GetQsub2())
1624 fHeader->fPsiRes = ep->GetQsub1()->Phi()/2.-ep->GetQsub2()->Phi()/2.;
1625 else
1626 fHeader->fPsiRes = 0;
1627 }
1628
788ca675 1629 Double_t val = 0;
1630 TString trgclasses(fHeader->fFiredTriggers);
1631 for (Int_t j = 0; j<fTrClassNamesArr->GetEntries(); ++j) {
1632 const char *name = fTrClassNamesArr->At(j)->GetName();
1633 if (trgclasses.Contains(name))
1634 val += TMath::Power(2,j);
1635 }
1636 fHeader->fTcls = (UInt_t)val;
1637
5fe1ca23 1638 fHeader->fNSelTr = fSelTracks->GetEntries();
1639 fHeader->fNSelPrimTr = fSelPrimTracks->GetEntries();
1640
1641 fHeader->fNCells = 0;
1642 fHeader->fNCells1 = 0;
1643 fHeader->fNCells2 = 0;
1644 fHeader->fNCells5 = 0;
1645 fHeader->fMaxCellE = 0;
1646
1647 AliVCaloCells *cells = fEsdCells;
1648 if (!cells)
1649 cells = fAodCells;
1650
1651 if (cells) {
1652 Int_t ncells = cells->GetNumberOfCells();
1653 for(Int_t j=0; j<ncells; ++j) {
1654 Double_t cellen = cells->GetAmplitude(j);
1655 if (cellen>1)
1656 ++fHeader->fNCells1;
1657 if (cellen>2)
1658 ++fHeader->fNCells2;
1659 if (cellen>5)
1660 ++fHeader->fNCells5;
1661 if (cellen>fHeader->fMaxCellE)
1662 fHeader->fMaxCellE = cellen;
1663 }
1664 fHeader->fNCells = ncells;
1665 }
1666
1667 fHeader->fNClus = 0;
1668 fHeader->fNClus1 = 0;
1669 fHeader->fNClus2 = 0;
1670 fHeader->fNClus5 = 0;
1671 fHeader->fMaxClusE = 0;
1672
1673 TObjArray *clusters = fEsdClusters;
1674 if (!clusters)
1675 clusters = fAodClusters;
1676
1677 if (clusters) {
1678 Int_t nclus = clusters->GetEntries();
1679 for(Int_t j=0; j<nclus; ++j) {
1680 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(j));
1681 Double_t clusen = clus->E();
1682 if (clusen>1)
1683 ++fHeader->fNClus1;
1684 if (clusen>2)
1685 ++fHeader->fNClus2;
1686 if (clusen>5)
1687 ++fHeader->fNClus5;
1688 if (clusen>fHeader->fMaxClusE)
1689 fHeader->fMaxClusE = clusen;
1690 }
1691 fHeader->fNClus = nclus;
1692 }
1693
1694
788ca675 1695 if (fAodEv) {
1696 am->LoadBranch("vertices");
1697 AliAODVertex *pv = fAodEv->GetPrimaryVertex();
1698 FillVertex(fPrimVert, pv);
1699 AliAODVertex *sv = fAodEv->GetPrimaryVertexSPD();
1700 FillVertex(fSpdVert, sv);
1701 } else {
1702 am->LoadBranch("PrimaryVertex.");
1703 const AliESDVertex *pv = fEsdEv->GetPrimaryVertexTracks();
1704 FillVertex(fPrimVert, pv);
1705 am->LoadBranch("SPDVertex.");
1706 const AliESDVertex *sv = fEsdEv->GetPrimaryVertexSPD();
1707 FillVertex(fSpdVert, sv);
1708 am->LoadBranch("TPCVertex.");
1709 const AliESDVertex *tv = fEsdEv->GetPrimaryVertexTPC();
1710 FillVertex(fTpcVert, tv);
fa443410 1711 }
788ca675 1712
788ca675 1713 fNtuple->Fill();
286b47a5 1714}
1715
1716//________________________________________________________________________
1717void AliAnalysisTaskEMCALPi0PbPb::FillPionHists()
1718{
1719 // Fill histograms related to pions.
286b47a5 1720
296ea9b4 1721 Double_t vertex[3] = {0};
fa443410 1722 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1723
90d5b88b 1724 TObjArray *clusters = fEsdClusters;
1725 if (!clusters)
1726 clusters = fAodClusters;
fa443410 1727
90d5b88b 1728 if (clusters) {
1729 TLorentzVector clusterVec1;
1730 TLorentzVector clusterVec2;
1731 TLorentzVector pionVec;
1732
1733 Int_t nclus = clusters->GetEntries();
fa443410 1734 for (Int_t i = 0; i<nclus; ++i) {
90d5b88b 1735 AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i));
fa443410 1736 if (!clus1)
1737 continue;
1738 if (!clus1->IsEMCAL())
1739 continue;
3c661da5 1740 if (clus1->E()<fMinE)
6eb6260e 1741 continue;
a49742b5 1742 if (clus1->GetNCells()<fNminCells)
1743 continue;
f224d35b 1744 if (GetMaxCellEnergy(clus1)/clus1->E()<fMinErat)
1745 continue;
3c661da5 1746 if (clus1->Chi2()<fMinEcc) // eccentricity cut
f224d35b 1747 continue;
fa443410 1748 clus1->GetMomentum(clusterVec1,vertex);
1749 for (Int_t j = i+1; j<nclus; ++j) {
90d5b88b 1750 AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j));
fa443410 1751 if (!clus2)
1752 continue;
1753 if (!clus2->IsEMCAL())
1754 continue;
3c661da5 1755 if (clus2->E()<fMinE)
6eb6260e 1756 continue;
a49742b5 1757 if (clus2->GetNCells()<fNminCells)
1758 continue;
f224d35b 1759 if (GetMaxCellEnergy(clus2)/clus2->E()<fMinErat)
1760 continue;
3c661da5 1761 if (clus2->Chi2()<fMinEcc) // eccentricity cut
f224d35b 1762 continue;
fa443410 1763 clus2->GetMomentum(clusterVec2,vertex);
1764 pionVec = clusterVec1 + clusterVec2;
1765 Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E();
6eb6260e 1766 Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect());
d595acbb 1767 if (pionZgg < fAsymMax) {
1768 fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi());
1769 fHPionMggPt->Fill(pionVec.M(),pionVec.Pt());
1770 fHPionMggAsym->Fill(pionVec.M(),pionZgg);
6eb6260e 1771 fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle);
d595acbb 1772 Int_t bin = fPtRanges->FindBin(pionVec.Pt());
1773 fHPionInvMasses[bin]->Fill(pionVec.M());
1774 }
fa443410 1775 }
1776 }
90d5b88b 1777 }
fa443410 1778}
d595acbb 1779
38727e64 1780//________________________________________________________________________
1781void AliAnalysisTaskEMCALPi0PbPb::FillMcHists()
1782{
1783 // Fill additional MC information histograms.
1784
1785 if (!fMcMode)
1786 return;
1787
1788 // check if aod or esd mc mode and the fill histos
1789}
1790
6eb6260e 1791//________________________________________________________________________
323834f0 1792void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists()
6eb6260e 1793{
788ca675 1794 // Fill other histograms.
1795}
1796
1797//__________________________________________________________________________________________________
1798void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliESDVertex *esdv)
1799{
1800 // Fill vertex from ESD vertex info.
1801
1802 v->fVx = esdv->GetX();
1803 v->fVy = esdv->GetY();
1804 v->fVz = esdv->GetZ();
1805 v->fVc = esdv->GetNContributors();
1806 v->fDisp = esdv->GetDispersion();
1807 v->fZres = esdv->GetZRes();
1808 v->fChi2 = esdv->GetChi2();
1809 v->fSt = esdv->GetStatus();
1810 v->fIs3D = esdv->IsFromVertexer3D();
1811 v->fIsZ = esdv->IsFromVertexerZ();
1812}
1813
1814//__________________________________________________________________________________________________
1815void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliAODVertex *aodv)
1816{
1817 // Fill vertex from AOD vertex info.
1818
1819 v->fVx = aodv->GetX();
1820 v->fVy = aodv->GetY();
1821 v->fVz = aodv->GetZ();
1822 v->fVc = aodv->GetNContributors();
1823 v->fChi2 = aodv->GetChi2();
6eb6260e 1824}
1825
d595acbb 1826//________________________________________________________________________
296ea9b4 1827Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
1828{
1829 // Compute isolation based on cell content.
1830
1831 AliVCaloCells *cells = fEsdCells;
1832 if (!cells)
1833 cells = fAodCells;
1834 if (!cells)
1835 return 0;
1836
1837 Double_t cellIsolation = 0;
1838 Double_t rad2 = radius*radius;
1839 Int_t ncells = cells->GetNumberOfCells();
1840 for (Int_t i = 0; i<ncells; ++i) {
1841 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
0fbe8d4f 1842 Float_t eta=-1, phi=-1;
296ea9b4 1843 fGeom->EtaPhiFromIndex(absID,eta,phi);
0fbe8d4f 1844 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
296ea9b4 1845 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1846 if(dist>rad2)
1847 continue;
0fbe8d4f 1848 Double_t cellE = cells->GetAmplitude(i);
296ea9b4 1849 cellIsolation += cellE;
1850 }
1851 return cellIsolation;
1852}
1853
1854//________________________________________________________________________
0fbe8d4f 1855Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellEnergy(const AliVCluster *cluster) const
1856{
1857 // Get maximum energy of attached cell.
1858
1859 Double_t ret = 0;
1860 Int_t ncells = cluster->GetNCells();
1861 if (fEsdCells) {
1862 for (Int_t i=0; i<ncells; i++) {
1863 Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1864 ret += e;
1865 }
1866 } else {
1867 for (Int_t i=0; i<ncells; i++) {
1868 Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1869 ret += e;
1870 }
1871 }
1872 return ret;
1873}
1874
1875//________________________________________________________________________
1876Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
d595acbb 1877{
90d5b88b 1878 // Get maximum energy of attached cell.
1879
788ca675 1880 id = -1;
90d5b88b 1881 Double_t maxe = 0;
90d5b88b 1882 Int_t ncells = cluster->GetNCells();
1883 if (fEsdCells) {
1884 for (Int_t i=0; i<ncells; i++) {
27422847 1885 Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
f0897c18 1886 if (e>maxe) {
90d5b88b 1887 maxe = e;
788ca675 1888 id = cluster->GetCellAbsId(i);
f0897c18 1889 }
90d5b88b 1890 }
1891 } else {
1892 for (Int_t i=0; i<ncells; i++) {
27422847 1893 Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
90d5b88b 1894 if (e>maxe)
1895 maxe = e;
788ca675 1896 id = cluster->GetCellAbsId(i);
90d5b88b 1897 }
1898 }
6eb6260e 1899 return maxe;
d595acbb 1900}
1901
1902//________________________________________________________________________
0fbe8d4f 1903void AliAnalysisTaskEMCALPi0PbPb::GetSigma(const AliVCluster *c, Double_t& sigmaMax, Double_t &sigmaMin) const
d595acbb 1904{
90d5b88b 1905 // Calculate the (E) weighted variance along the longer (eigen) axis.
1906
6bf90832 1907 sigmaMax = 0; // cluster variance along its longer axis
1908 sigmaMin = 0; // cluster variance along its shorter axis
1909 Double_t Ec = c->E(); // cluster energy
296ea9b4 1910 if(Ec<=0)
1911 return;
6bf90832 1912 Double_t Xc = 0 ; // cluster first moment along X
1913 Double_t Yc = 0 ; // cluster first moment along Y
1914 Double_t Sxx = 0 ; // cluster second central moment along X (variance_X^2)
1915 Double_t Sxy = 0 ; // cluster second central moment along Y (variance_Y^2)
1916 Double_t Syy = 0 ; // cluster covariance^2
90d5b88b 1917
1918 AliVCaloCells *cells = fEsdCells;
1919 if (!cells)
1920 cells = fAodCells;
1921
6eb6260e 1922 if (!cells)
1923 return;
1924
6bf90832 1925 Int_t ncells = c->GetNCells();
6eb6260e 1926 if (ncells==1)
1927 return;
1928
6eb6260e 1929 for(Int_t j=0; j<ncells; ++j) {
6bf90832 1930 Int_t id = TMath::Abs(c->GetCellAbsId(j));
6eb6260e 1931 Double_t cellen = cells->GetCellAmplitude(id);
3a952328 1932 TVector3 pos;
1933 fGeom->GetGlobal(id,pos);
6eb6260e 1934 Xc += cellen*pos.X();
1935 Yc += cellen*pos.Y();
1936 Sxx += cellen*pos.X()*pos.X();
1937 Syy += cellen*pos.Y()*pos.Y();
1938 Sxy += cellen*pos.X()*pos.Y();
1939 }
1940 Xc /= Ec;
1941 Yc /= Ec;
1942 Sxx /= Ec;
1943 Syy /= Ec;
1944 Sxy /= Ec;
1945 Sxx -= Xc*Xc;
1946 Syy -= Yc*Yc;
1947 Sxy -= Xc*Yc;
f0897c18 1948 Sxx = TMath::Abs(Sxx);
1949 Syy = TMath::Abs(Syy);
296ea9b4 1950 sigmaMax = (Sxx + Syy + TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
1951 sigmaMax = TMath::Sqrt(TMath::Abs(sigmaMax));
1952 sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
1953 sigmaMin = TMath::Sqrt(TMath::Abs(sigmaMin));
d595acbb 1954}
90d5b88b 1955
6bf90832 1956//________________________________________________________________________
0fbe8d4f 1957Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(const AliVCluster *c, Double_t emin) const
6bf90832 1958{
1959 // Calculate number of attached cells above emin.
1960
6bf90832 1961 AliVCaloCells *cells = fEsdCells;
1962 if (!cells)
1963 cells = fAodCells;
6bf90832 1964 if (!cells)
296ea9b4 1965 return 0;
6bf90832 1966
296ea9b4 1967 Int_t n = 0;
6bf90832 1968 Int_t ncells = c->GetNCells();
1969 for(Int_t j=0; j<ncells; ++j) {
1970 Int_t id = TMath::Abs(c->GetCellAbsId(j));
1971 Double_t cellen = cells->GetCellAmplitude(id);
1972 if (cellen>=emin)
1973 ++n;
1974 }
1975 return n;
1976}
1977
296ea9b4 1978//________________________________________________________________________
b6c599fe 1979Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
296ea9b4 1980{
1981 // Compute isolation based on tracks.
1982
1983 Double_t trkIsolation = 0;
1984 Double_t rad2 = radius*radius;
b6c599fe 1985 Int_t ntrks = fSelPrimTracks->GetEntries();
296ea9b4 1986 for(Int_t j = 0; j<ntrks; ++j) {
1987 AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
1988 if (!track)
1989 continue;
b6c599fe 1990 if (track->Pt()<pt)
1991 continue;
296ea9b4 1992 Float_t eta = track->Eta();
1993 Float_t phi = track->Phi();
0fbe8d4f 1994 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
296ea9b4 1995 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1996 if(dist>rad2)
1997 continue;
1998 trkIsolation += track->Pt();
1999 }
2000 return trkIsolation;
2001}
3a952328 2002
0fbe8d4f 2003//________________________________________________________________________
2004Bool_t AliAnalysisTaskEMCALPi0PbPb::IsShared(const AliVCluster *c) const
2005{
2006 // Returns if cluster shared across super modules.
2007
2008 AliVCaloCells *cells = fEsdCells;
2009 if (!cells)
2010 cells = fAodCells;
2011 if (!cells)
2012 return 0;
2013
2014 Int_t n = -1;
2015 Int_t ncells = c->GetNCells();
2016 for(Int_t j=0; j<ncells; ++j) {
2017 Int_t id = TMath::Abs(c->GetCellAbsId(j));
2018 Int_t got = id / (24*48);
2019 if (n==-1) {
2020 n = got;
2021 continue;
2022 }
2023 if (got!=n)
2024 return 1;
2025 }
2026 return 0;
2027}
3a952328 2028
38727e64 2029//________________________________________________________________________
2030void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(AliMCParticle *p) const
2031{
2032 // Print daugher information. TODO
2033}
2034
2035//________________________________________________________________________
2036void AliAnalysisTaskEMCALPi0PbPb::PrintTrackRefs(AliMCParticle *p) const
2037{
2038 // Print track ref array. TODO
2039}
2040