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