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