3 #include "AliAnalysisTaskEMCALPi0PbPb.h"
6 #include <TClonesArray.h>
8 #include <TGeoGlobalMagField.h>
12 #include <TLorentzVector.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 "AliEMCALRecoUtils.h"
26 #include "AliESDEvent.h"
27 #include "AliESDVertex.h"
28 #include "AliESDtrack.h"
29 #include "AliESDtrackCuts.h"
30 #include "AliGeomManager.h"
31 #include "AliInputEventHandler.h"
34 #include "AliTrackerBase.h"
36 ClassImp(AliAnalysisTaskEMCALPi0PbPb)
38 //________________________________________________________________________
39 AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name)
40 : AliAnalysisTaskSE(name),
55 fGeoName("EMCAL_FIRSTYEARV1"),
56 fMinNClustPerTrack(50),
61 fDoTrackMatWithGeom(0),
89 fHTclsBeforeCuts(0x0),
99 fHCellFreqCut100M(0x0),
100 fHCellFreqCut300M(0x0),
103 fHClustEccentricity(0),
105 fHClustEnergyPt(0x0),
106 fHClustEnergySigma(0x0),
107 fHClustSigmaSigma(0x0),
108 fHClustNCellEnergyRatio(0x0),
119 DefineInput(0, TChain::Class());
120 DefineOutput(1, TList::Class());
121 fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,EMCALCells.,Tracks "
122 "AOD:header,vertices,emcalCells,tracks";
125 //________________________________________________________________________
126 AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb()
130 if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
131 delete fOutput; fOutput = 0;
133 delete fPtRanges; fPtRanges = 0;
134 delete fGeom; fGeom = 0;
135 delete fReco; fReco = 0;
136 delete fTrClassNamesArr;
139 delete [] fHColuRowE;
140 delete [] fHCellMult;
141 delete [] fHCellFreqNoCut;
142 delete [] fHCellFreqCut100M;
143 delete [] fHCellFreqCut300M;
146 //________________________________________________________________________
147 void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects()
149 // Create user objects here.
151 fGeom = new AliEMCALGeoUtils(fGeoName,"EMCAL");
152 fReco = new AliEMCALRecoUtils();
153 fTrClassNamesArr = fTrClassNames.Tokenize(" ");
154 fOutput = new TList();
156 fSelTracks = new TObjArray;
159 TFile *f = OpenFile(1);
161 f->SetCompressionLevel(2);
162 fNtuple = new TTree(Form("tree%.0fto%.0f",fCentFrom,fCentTo), "StandaloneTree");
163 fNtuple->SetDirectory(f);
164 fNtuple->SetAutoFlush(-1024*1024*1024);
165 fNtuple->SetAutoSave(-1024*1024*1024);
166 fHeader = new AliStaHeader;
167 fNtuple->Branch("header", &fHeader, 16*1024, 99);
168 fPrimVert = new AliStaVertex;
169 fNtuple->Branch("primv", &fPrimVert, 16*1024, 99);
170 fSpdVert = new AliStaVertex;
171 fNtuple->Branch("spdv", &fSpdVert, 16*1024, 99);
172 fTpcVert = new AliStaVertex;
173 fNtuple->Branch("tpcv", &fTpcVert, 16*1024, 99);
174 if (TClass::GetClass("AliStaCluster"))
175 TClass::GetClass("AliStaCluster")->IgnoreTObjectStreamer();
176 fClusters = new TClonesArray("AliStaCluster",1000);
177 fNtuple->Branch("clusters", &fClusters, 8*16*1024, 99);
181 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
182 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-0.25;
183 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+0.25;
186 TH1::SetDefaultSumw2(kTRUE);
187 TH2::SetDefaultSumw2(kTRUE);
188 fHCuts = new TH1F("hEventCuts","",5,0.5,5.5);
189 fHCuts->GetXaxis()->SetBinLabel(1,"All");
190 fHCuts->GetXaxis()->SetBinLabel(2,"PS");
191 fHCuts->GetXaxis()->SetBinLabel(3,Form("%s: %.0f-%.0f",fCentVar.Data(),fCentFrom,fCentTo));
192 fHCuts->GetXaxis()->SetBinLabel(4,"QFlag");
193 fHCuts->GetXaxis()->SetBinLabel(5,Form("zvtx: %.0f-%.0f",fVtxZMin,fVtxZMax));
194 fOutput->Add(fHCuts);
195 fHVertexZ = new TH1F("hVertexZBeforeCut","",100,-25,25);
196 fHVertexZ->SetXTitle("z [cm]");
197 fOutput->Add(fHVertexZ);
198 fHVertexZ2 = new TH1F("hVertexZAfterCut","",100,-25,25);
199 fHVertexZ2->SetXTitle("z [cm]");
200 fOutput->Add(fHVertexZ2);
201 fHCent = new TH1F("hCentBeforeCut","",101,-1,100);
202 fHCent->SetXTitle(fCentVar.Data());
203 fOutput->Add(fHCent);
204 fHCentQual = new TH1F("hCentAfterCut","",101,-1,100);
205 fHCentQual->SetXTitle(fCentVar.Data());
206 fOutput->Add(fHCentQual);
207 fHTclsBeforeCuts = new TH1F("hTclsBeforeCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
208 fHTclsAfterCuts = new TH1F("hTclsAfterCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
209 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
210 const char *name = fTrClassNamesArr->At(i)->GetName();
211 fHTclsBeforeCuts->GetXaxis()->SetBinLabel(1+i,name);
212 fHTclsAfterCuts->GetXaxis()->SetBinLabel(1+i,name);
214 fOutput->Add(fHTclsBeforeCuts);
215 fOutput->Add(fHTclsAfterCuts);
217 // histograms for cells
218 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
219 fHColuRow = new TH2*[nsm];
220 fHColuRowE = new TH2*[nsm];
221 fHCellMult = new TH1*[nsm];
222 for (Int_t i = 0; i<nsm; ++i) {
223 fHColuRow[i] = new TH2F(Form("hColRow_Mod%d",i),"",48,0,48,24,0.,24);
224 fHColuRow[i]->SetTitle(Form("Module %d: Occupancy", i));
225 fHColuRow[i]->SetXTitle("col (i#eta)");
226 fHColuRow[i]->SetYTitle("row (i#phi)");
227 fHColuRowE[i] = new TH2F(Form("hColRowE_Mod%d", i),"",48,0,48,24,0,24);
228 fHColuRowE[i]->SetTitle(Form("Module %d: Cell energy",i));
229 fHColuRowE[i]->SetXTitle("col (i#eta)");
230 fHColuRowE[i]->SetYTitle("row (i#phi)");
231 fHCellMult[i] = new TH1F(Form("hCellMult_Mod%d",i),"",1000,0,1000);
232 fHCellMult[i]->SetTitle(Form("Module %d: Cell multiplicity",i));
233 fHCellMult[i]->SetXTitle("# of cells");
234 fOutput->Add(fHColuRow[i]);
235 fOutput->Add(fHColuRowE[i]);
236 fOutput->Add(fHCellMult[i]);
238 fHCellE = new TH1F("hCellE","",250,0.,25.);
239 fHCellE->SetXTitle("E_{cell} [GeV]");
240 fOutput->Add(fHCellE);
241 fHCellH = new TH1F ("hCellHighestE","",250,0.,25.);
242 fHCellH->SetXTitle("E^{max}_{cell} [GeV]");
243 fOutput->Add(fHCellH);
244 fHCellM = new TH1F ("hCellMeanEperHitCell","",250,0.,2.5);
245 fHCellM->SetXTitle("#LT E_{cell}#GT [GeV]");
246 fOutput->Add(fHCellM);
247 fHCellM2 = new TH1F ("hCellMeanEperAllCells","",250,0.,1);
248 fHCellM2->SetXTitle("1/N_{cells} #Sigma E_{cell} [GeV]");
249 fOutput->Add(fHCellM2);
251 fHCellFreqNoCut = new TH1*[nsm];
252 fHCellFreqCut100M = new TH1*[nsm];
253 fHCellFreqCut300M = new TH1*[nsm];
254 fHCellFreqE = new TH1*[nsm];
255 for (Int_t i = 0; i<nsm; ++i){
256 Double_t lbin = i*24*48-0.5;
257 Double_t hbin = lbin+24*48;
258 fHCellFreqNoCut[i] = new TH1F(Form("hCellFreqNoCut_SM%d",i),
259 Form("Frequency SM%d (no cut);id;#",i), 1152, lbin, hbin);
260 fHCellFreqCut100M[i] = new TH1F(Form("hCellFreqCut100M_SM%d",i),
261 Form("Frequency SM%d (>0.1GeV);id;#",i), 1152, lbin, hbin);
262 fHCellFreqCut300M[i] = new TH1F(Form("hCellFreqCut300M_SM%d",i),
263 Form("Frequency SM%d (>0.3GeV);id;#",i), 1152, lbin, hbin);
264 fHCellFreqE[i] = new TH1F(Form("hCellFreqE_SM%d",i),
265 Form("Frequency SM%d (E weighted);id;#",i), 1152, lbin, hbin);
266 fOutput->Add(fHCellFreqNoCut[i]);
267 fOutput->Add(fHCellFreqCut100M[i]);
268 fOutput->Add(fHCellFreqCut300M[i]);
269 fOutput->Add(fHCellFreqE[i]);
272 fHCellCheckE = new TH1*[24*48*nsm];
273 memset(fHCellCheckE,0,24*48*nsm*sizeof(TH1*));
274 Int_t tcs[1] = {4102};
275 for (UInt_t i = 0; i<sizeof(tcs)/sizeof(Int_t); ++i){
278 fHCellCheckE[i] = new TH1F(Form("hCellE_id%d",c), Form("Cell %d;E [GeV/c];#",c), 500, 0, 8);
279 fOutput->Add(fHCellCheckE[i]);
284 // histograms for clusters
285 fHClustEccentricity = new TH1F("hClustEccentricity","",100,-0.1,1.1);
286 fHClustEccentricity->SetXTitle("#epsilon_{C}");
287 fOutput->Add(fHClustEccentricity);
288 fHClustEtaPhi = new TH2F("hClustEtaPhi","",500,-0.8,0.8,100*nsm,phimin,phimax);
289 fHClustEtaPhi->SetXTitle("#eta");
290 fHClustEtaPhi->SetYTitle("#varphi");
291 fOutput->Add(fHClustEtaPhi);
292 fHClustEnergyPt = new TH2F("hClustEnergyPt","",250,0,50,250,0,50);
293 fHClustEnergyPt->SetXTitle("E [GeV]");
294 fHClustEnergyPt->SetYTitle("p_{T} [GeV/c]");
295 fOutput->Add(fHClustEnergyPt);
296 fHClustEnergySigma = new TH2F("hClustEnergySigma","",100,0,100,500,0,50);
297 fHClustEnergySigma->SetXTitle("E_{C} * #sigma_{max} [GeV*cm]");
298 fHClustEnergySigma->SetYTitle("E_{C} [GeV]");
299 fOutput->Add(fHClustEnergySigma);
300 fHClustSigmaSigma = new TH2F("hClustSigmaSigma","",500,0,50,500,0,50);
301 fHClustSigmaSigma->SetXTitle("#lambda_{0} [cm]");
302 fHClustSigmaSigma->SetYTitle("#sigma_{max} [cm]");
303 fOutput->Add(fHClustSigmaSigma);
304 fHClustNCellEnergyRatio = new TH2F("hClustNCellEnergyRatio","",27,-0.5,26.5,101,-0.05,1.05);
305 fHClustNCellEnergyRatio->SetXTitle("N_{cells}");
306 fHClustNCellEnergyRatio->SetYTitle("E^{max}_{cell}/E_{clus}");
307 fOutput->Add(fHClustNCellEnergyRatio);
309 // histograms for pion candidates
310 fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100*nsm,phimin,phimax);
311 fHPionEtaPhi->SetXTitle("#eta_{#gamma#gamma}");
312 fHPionEtaPhi->SetYTitle("#varphi_{#gamma#gamma}");
313 fOutput->Add(fHPionEtaPhi);
314 fHPionMggPt = new TH2F("hPionMggPt","",1000,0,2,100,0,20.0);
315 fHPionMggPt->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
316 fHPionMggPt->SetYTitle("p_{T}^{#gamma#gamma} [GeV/c]");
317 fOutput->Add(fHPionMggPt);
318 fHPionMggAsym = new TH2F("hPionMggAsym","",1000,0,2,100,0,1);
319 fHPionMggAsym->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
320 fHPionMggAsym->SetYTitle("Z_{#gamma#gamma} [GeV]");
321 fOutput->Add(fHPionMggAsym);
322 fHPionMggDgg = new TH2F("hPionMggDgg","",1000,0,2,100,0,10);
323 fHPionMggDgg->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
324 fHPionMggDgg->SetYTitle("opening angle [grad]");
325 fOutput->Add(fHPionMggDgg);
326 const Int_t nbins = 20;
327 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};
328 fPtRanges = new TAxis(nbins-1,xbins);
329 for (Int_t i = 0; i<=nbins; ++i) {
330 fHPionInvMasses[i] = new TH1F(Form("hPionInvMass%d",i),"",1000,0,2);
331 fHPionInvMasses[i]->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
333 fHPionInvMasses[i]->SetTitle(Form("0 < p_{T}^{#gamma#gamma} <%.1f",xbins[0]));
335 fHPionInvMasses[i]->SetTitle(Form("p_{T}^{#gamma#gamma} > 50"));
337 fHPionInvMasses[i]->SetTitle(Form("%.1f < p_{T}^{#gamma#gamma} <%.1f",xbins[i-1],xbins[i]));
338 fOutput->Add(fHPionInvMasses[i]);
341 PostData(1, fOutput);
344 //________________________________________________________________________
345 void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *)
347 // Called for each event.
352 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
353 fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent());
355 am->LoadBranch("AliESDRun.");
356 am->LoadBranch("AliESDHeader.");
358 fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
360 AliFatal("Neither ESD nor AOD event found");
363 am->LoadBranch("header");
366 if (fDoTrackMatWithGeom && !AliGeomManager::GetGeometry()) { // get geometry
367 AliWarning("Accessing geometry from OCDB, this is not very efficient!");
368 AliCDBManager *cdb = AliCDBManager::Instance();
369 if (!cdb->IsDefaultStorageSet())
370 cdb->SetDefaultStorage("raw://");
371 Int_t runno = InputEvent()->GetRunNumber();
372 if (runno != cdb->GetRun())
374 AliGeomManager::LoadGeometry();
377 if (!AliGeomManager::GetGeometry()&&!fIsGeoMatsSet) { // set misalignment matrices (stored in first event)
378 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
380 for (Int_t i=0; i<nsm; ++i)
381 fGeom->SetMisalMatrix(fEsdEv->GetESDRun()->GetEMCALMatrix(i),i);
383 for (Int_t i=0; i<nsm; ++i)
384 fGeom->SetMisalMatrix(fAodEv->GetHeader()->GetEMCALMatrix(i),i);
386 fIsGeoMatsSet = kTRUE;
389 if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
391 const AliESDRun *erun = fEsdEv->GetESDRun();
392 AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(),
393 erun->GetCurrentDip(),
396 erun->GetBeamEnergy(),
397 erun->GetBeamType());
398 TGeoGlobalMagField::Instance()->SetField(field);
400 Double_t pol = -1; //polarity
401 Double_t be = -1; //beam energy
402 AliMagF::BeamType_t btype = AliMagF::kBeamTypepp;
403 Int_t runno = fAodEv->GetRunNumber();
404 if (runno>=136851 && runno<138275) {
407 btype = AliMagF::kBeamTypeAA;
408 } else if (runno>=138275 && runno<=139517) {
411 btype = AliMagF::kBeamTypeAA;
413 AliError(Form("Do not know the bfield parameters for run %d! Using defaults!!!", runno));
415 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", pol, pol, AliMagF::k5kG, btype, be));
423 AliESDHeader *h = dynamic_cast<AliESDHeader*>(InputEvent()->GetHeader());
425 trgclasses = fEsdEv->GetFiredTriggerClasses();
427 AliAODHeader *h2 = dynamic_cast<AliAODHeader*>(InputEvent()->GetHeader());
429 trgclasses = h2->GetFiredTriggerClasses();
431 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
432 const char *name = fTrClassNamesArr->At(i)->GetName();
433 if (trgclasses.Contains(name))
434 fHTclsBeforeCuts->Fill(1+i);
437 UInt_t res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
443 const AliCentrality *centP = InputEvent()->GetCentrality();
444 Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar);
446 if (cent<fCentFrom||cent>fCentTo)
452 if (centP->GetQuality()>0)
456 fHCentQual->Fill(cent);
460 am->LoadBranch("PrimaryVertex.");
461 am->LoadBranch("SPDVertex.");
462 am->LoadBranch("TPCVertex.");
464 fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
465 am->LoadBranch("vertices");
469 const AliVVertex *vertex = InputEvent()->GetPrimaryVertex();
473 fHVertexZ->Fill(vertex->GetZ());
475 if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax)
479 fHVertexZ2->Fill(vertex->GetZ());
481 // count number of accepted events
484 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
485 const char *name = fTrClassNamesArr->At(i)->GetName();
486 if (trgclasses.Contains(name))
487 fHTclsAfterCuts->Fill(1+i);
490 fRecPoints = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
491 fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set of if clusters are attached
492 fEsdCells = 0; // will be set if ESD input used
493 fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set of if clusters are attached
494 // or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode
495 fAodCells = 0; // will be set if AOD input used
497 // deal with special output from AliAnalysisTaskEMCALClusterizeFast first
498 Bool_t clusattached = 0;
499 Bool_t recalibrated = 0;
500 if (1 && !fClusName.IsNull()) {
501 AliAnalysisTaskEMCALClusterizeFast *cltask = 0;
502 TObjArray *ts = am->GetTasks();
503 cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName));
504 if (cltask && cltask->GetClusters()) {
505 fRecPoints = const_cast<TObjArray*>(cltask->GetClusters());
506 clusattached = cltask->GetAttachClusters();
507 if (cltask->GetCalibData()!=0)
508 recalibrated = kTRUE;
511 if (1 && AODEvent() && !fClusName.IsNull()) {
512 TList *l = AODEvent()->GetList();
513 TClonesArray *clus = 0;
515 clus = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
520 if (fEsdEv) { // ESD input mode
521 if (1 && (!fRecPoints||clusattached)) {
523 am->LoadBranch("CaloClusters");
524 TList *l = fEsdEv->GetList();
525 TClonesArray *clus = 0;
527 clus = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
533 am->LoadBranch("EMCALCells.");
534 fEsdCells = fEsdEv->GetEMCALCells();
536 } else if (fAodEv) { // AOD input mode
537 if (1 && (!fAodClusters || clusattached)) {
539 am->LoadBranch("caloClusters");
540 TList *l = fAodEv->GetList();
541 TClonesArray *clus = 0;
543 clus = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
549 am->LoadBranch("emcalCells");
550 fAodCells = fAodEv->GetEMCALCells();
553 AliFatal("Impossible to not have either pointer to ESD or AOD event");
557 AliDebug(2,Form("fRecPoints set: %p", fRecPoints));
558 AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters));
559 AliDebug(2,Form("fEsdCells set: %p", fEsdCells));
560 AliDebug(2,Form("fAodClusters set: %p", fAodClusters));
561 AliDebug(2,Form("fAodCells set: %p", fAodCells));
565 ClusterAfterburner();
576 PostData(1, fOutput);
579 //________________________________________________________________________
580 void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *)
582 // Terminate called at the end of analysis.
585 TFile *f = OpenFile(1);
590 AliInfo(Form("%s: Accepted %lld events", GetName(), fNEvs));
593 //________________________________________________________________________
594 void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
596 // Calculate track properties.
600 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
601 TClonesArray *tracks = 0;
603 am->LoadBranch("Tracks");
604 TList *l = fEsdEv->GetList();
605 tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
607 am->LoadBranch("tracks");
608 TList *l = fAodEv->GetList();
609 tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
615 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
616 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-0.25;
617 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+0.25;
621 fSelTracks->SetOwner(kTRUE);
622 am->LoadBranch("PrimaryVertex.");
623 const AliESDVertex *vtx = fEsdEv->GetPrimaryVertexSPD();
624 am->LoadBranch("SPDVertex.");
625 const AliESDVertex *vtxSPD = fEsdEv->GetPrimaryVertexSPD();
626 am->LoadBranch("Tracks");
627 const Int_t Ntracks = tracks->GetEntries();
628 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
629 AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
631 AliWarning(Form("Could not receive track %d\n", iTracks));
634 if (fTrCuts && !fTrCuts->IsSelected(track))
636 Double_t eta = track->Eta();
639 if (track->Phi()<phimin||track->Phi()>phimax)
641 if(track->GetTPCNcls()<fMinNClustPerTrack)
645 fSelTracks->Add(track);
649 AliESDtrack copyt(*track);
651 copyt.GetBxByBz(bfield);
652 AliExternalTrackParam tParam;
653 Bool_t relate = copyt.RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&tParam);
656 copyt.Set(tParam.GetX(),tParam.GetAlpha(),tParam.GetParameter(),tParam.GetCovariance());
658 Double_t p[3] = { 0. };
660 Double_t pos[3] = { 0. };
662 Double_t covTr[21] = { 0. };
663 copyt.GetCovarianceXYZPxPyPz(covTr);
664 Double_t pid[10] = { 0. };
665 copyt.GetESDpid(pid);
666 AliAODTrack *aTrack = new AliAODTrack(copyt.GetID(),
673 (Short_t)copyt.GetSign(),
674 copyt.GetITSClusterMap(),
676 0,/*fPrimaryVertex,*/
677 kTRUE, // check if this is right
678 vtx->UsesTrack(copyt.GetID()));
679 aTrack->SetTPCClusterMap(copyt.GetTPCClusterMap());
680 aTrack->SetTPCSharedMap (copyt.GetTPCSharedMap());
681 Float_t ndf = copyt.GetTPCNcls() + 1 - 5 ;
683 aTrack->SetChi2perNDF(copyt.GetTPCchi2()/ndf);
685 aTrack->SetChi2perNDF(-1);
686 aTrack->SetFlags(copyt.GetStatus());
687 aTrack->SetTPCPointsF(copyt.GetTPCNclsF());
688 fSelTracks->Add(aTrack);
691 Int_t ntracks = tracks->GetEntries();
692 for (Int_t i=0; i<ntracks; ++i) {
693 AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
696 Double_t eta = track->Eta();
699 if (track->Phi()<phimin||track->Phi()>phimax)
701 if(track->GetTPCNcls()<fMinNClustPerTrack)
704 if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL
705 AliExternalTrackParam tParam(track);
706 if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 438, 0.139, 1, kTRUE)) {
708 tParam.GetXYZ(trkPos);
709 track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]);
712 fSelTracks->Add(track);
717 //________________________________________________________________________
718 void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
720 // Calculate cluster properties
722 TObjArray *clusters = fEsdClusters;
724 clusters = fAodClusters;
728 Int_t nclus = clusters->GetEntries();
729 Int_t ntrks = fSelTracks->GetEntries();
731 for(Int_t i = 0; i<nclus; ++i) {
732 fClusProps[i].Reset();
734 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
737 if (!clus->IsEMCAL())
742 Float_t clsPos[3] = {0};
743 clus->GetPosition(clsPos);
744 TVector3 clsVec(clsPos);
745 Double_t vertex[3] = {0};
746 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
747 TLorentzVector clusterVec;
748 clus->GetMomentum(clusterVec,vertex);
749 Double_t clsEta = clusterVec.Eta();
751 Double_t mind2 = 1e10;
752 for(Int_t j = 0; j<ntrks; ++j) {
753 AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
756 Double_t pt = track->Pt();
757 if (pt<fMinPtPerTrack)
759 if (TMath::Abs(clsEta-track->Eta())>0.5)
762 Float_t tmpR=-1, tmpZ=-1;
763 if (!fDoTrackMatWithGeom) {
765 AliExternalTrackParam *tParam = 0;
767 AliESDtrack *esdTrack = static_cast<AliESDtrack*>(track);
768 tParam = new AliExternalTrackParam(*esdTrack->GetTPCInnerParam());
770 tParam = new AliExternalTrackParam(track);
773 track->GetBxByBz(bfield);
774 TVector3 vec(clsPos);
775 Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
776 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
777 tParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
778 Bool_t ret = tParam->PropagateToBxByBz(vec.X(), bfield);
784 tParam->GetXYZ(trkPos); //Get the extrapolated global position
785 tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
786 tmpZ = clsPos[2]-trkPos[2];
789 if (TMath::Abs(clsEta-track->Eta())>fIsoDist)
791 AliExternalTrackParam tParam(track);
792 if (!fReco->ExtrapolateTrackToCluster(&tParam, clus, tmpR, tmpZ))
799 fClusProps[i].fTrIndex = j;
800 fClusProps[i].fTrDz = tmpZ;
801 fClusProps[i].fTrDr = TMath::Sqrt(tmpR*tmpR-tmpZ*tmpZ);
802 fClusProps[i].fTrDist = d2;
803 fClusProps[i].fTrEp = clus->E()/track->P();
807 if (0 && (fClusProps[i].fTrIndex>=0)) {
808 cout << i << " " << fClusProps[i].fTrIndex << ": Dr " << fClusProps[i].fTrDr << " " << " Dz " << fClusProps[i].fTrDz << endl;
811 fClusProps[i].fTrIso = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
812 fClusProps[i].fTrLowPtIso = 0;
813 fClusProps[i].fCellIso = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
817 //________________________________________________________________________
818 void AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner()
820 // Run custer reconstruction afterburner.
822 AliVCaloCells *cells = fEsdCells;
829 Int_t ncells = cells->GetNumberOfCells();
833 Double_t cellMeanE = 0, cellSigE = 0;
834 for (Int_t i = 0; i<ncells; ++i) {
835 Double_t cellE = cells->GetAmplitude(i);
837 cellSigE += cellE*cellE;
841 cellSigE -= cellMeanE*cellMeanE;
844 cellSigE = TMath::Sqrt(cellSigE / ncells);
846 Double_t subE = cellMeanE - 7*cellSigE;
850 for (Int_t i = 0; i<ncells; ++i) {
852 Double_t amp=0,time=0;
853 if (!cells->GetCell(i, id, amp, time))
858 cells->SetCell(i, id, amp, time);
861 TObjArray *clusters = fEsdClusters;
863 clusters = fAodClusters;
867 Int_t nclus = clusters->GetEntries();
868 for (Int_t i = 0; i<nclus; ++i) {
869 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
870 if (!clus->IsEMCAL())
872 Int_t nc = clus->GetNCells();
874 UShort_t ids[100] = {0};
875 Double_t fra[100] = {0};
876 for (Int_t j = 0; j<nc; ++j) {
877 Short_t id = TMath::Abs(clus->GetCellAbsId(j));
878 Double_t cen = cells->GetCellAmplitude(id);
886 clusters->RemoveAt(i);
890 for (Int_t j = 0; j<nc; ++j) {
892 Double_t cen = cells->GetCellAmplitude(id);
896 AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus);
899 aodclus->SetNCells(nc);
900 aodclus->SetCellsAmplitudeFraction(fra);
901 aodclus->SetCellsAbsId(ids);
904 clusters->Compress();
907 //________________________________________________________________________
908 void AliAnalysisTaskEMCALPi0PbPb::FillCellHists()
910 // Fill histograms related to cell properties.
912 AliVCaloCells *cells = fEsdCells;
919 Int_t cellModCount[12] = {0};
920 Double_t cellMaxE = 0;
921 Double_t cellMeanE = 0;
922 Int_t ncells = cells->GetNumberOfCells();
926 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
928 for (Int_t i = 0; i<ncells; ++i) {
929 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
930 Double_t cellE = cells->GetAmplitude(i);
931 fHCellE->Fill(cellE);
936 Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1;
937 Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta);
939 AliError(Form("Could not get cell index for %d", absID));
943 Int_t iPhi=-1, iEta=-1;
944 fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta);
945 fHColuRow[iSM]->Fill(iEta,iPhi,1);
946 fHColuRowE[iSM]->Fill(iEta,iPhi,cellE);
947 fHCellFreqNoCut[iSM]->Fill(absID);
948 if (cellE > 0.1) fHCellFreqCut100M[iSM]->Fill(absID);
949 if (cellE > 0.3) fHCellFreqCut300M[iSM]->Fill(absID);
950 if (fHCellCheckE && fHCellCheckE[absID])
951 fHCellCheckE[absID]->Fill(cellE);
952 fHCellFreqE[iSM]->Fill(absID, cellE);
954 fHCellH->Fill(cellMaxE);
956 fHCellM->Fill(cellMeanE);
957 fHCellM2->Fill(cellMeanE*ncells/24/48/nsm);
958 for (Int_t i=0; i<nsm; ++i)
959 fHCellMult[i]->Fill(cellModCount[i]);
962 //________________________________________________________________________
963 void AliAnalysisTaskEMCALPi0PbPb::FillClusHists()
965 // Fill histograms related to cluster properties.
967 TObjArray *clusters = fEsdClusters;
969 clusters = fAodClusters;
973 Double_t vertex[3] = {0};
974 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
976 Int_t nclus = clusters->GetEntries();
977 for(Int_t i = 0; i<nclus; ++i) {
978 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
981 if (!clus->IsEMCAL())
983 TLorentzVector clusterVec;
984 clus->GetMomentum(clusterVec,vertex);
985 Double_t maxAxis = 0, minAxis = 0;
986 GetSigma(clus,maxAxis,minAxis);
987 clus->SetTOF(maxAxis); // store sigma in TOF
988 Double_t clusterEcc = 0;
990 clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
991 clus->SetChi2(clusterEcc); // store ecc in chi2
992 fHClustEccentricity->Fill(clusterEcc);
993 fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi());
994 fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt());
995 fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E());
996 fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis);
997 fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E());
1001 //________________________________________________________________________
1002 void AliAnalysisTaskEMCALPi0PbPb::FillNtuple()
1009 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1011 fHeader->fRun = fAodEv->GetRunNumber();
1012 fHeader->fOrbit = fAodEv->GetHeader()->GetOrbitNumber();
1013 fHeader->fPeriod = fAodEv->GetHeader()->GetPeriodNumber();
1014 fHeader->fBx = fAodEv->GetHeader()->GetBunchCrossNumber();
1015 fHeader->fL0 = fAodEv->GetHeader()->GetL0TriggerInputs();
1016 fHeader->fL1 = fAodEv->GetHeader()->GetL1TriggerInputs();
1017 fHeader->fL2 = fAodEv->GetHeader()->GetL2TriggerInputs();
1018 fHeader->fTrClassMask = fAodEv->GetHeader()->GetTriggerMask();
1019 fHeader->fTrCluster = fAodEv->GetHeader()->GetTriggerCluster();
1020 fHeader->fOffTriggers = fAodEv->GetHeader()->GetOfflineTrigger();
1021 fHeader->fFiredTriggers = fAodEv->GetHeader()->GetFiredTriggerClasses();
1023 fHeader->fRun = fEsdEv->GetRunNumber();
1024 fHeader->fOrbit = fEsdEv->GetHeader()->GetOrbitNumber();
1025 fHeader->fPeriod = fEsdEv->GetESDRun()->GetPeriodNumber();
1026 fHeader->fBx = fEsdEv->GetHeader()->GetBunchCrossNumber();
1027 fHeader->fL0 = fEsdEv->GetHeader()->GetL0TriggerInputs();
1028 fHeader->fL1 = fEsdEv->GetHeader()->GetL1TriggerInputs();
1029 fHeader->fL2 = fEsdEv->GetHeader()->GetL2TriggerInputs();
1030 fHeader->fTrClassMask = fEsdEv->GetHeader()->GetTriggerMask();
1031 fHeader->fTrCluster = fEsdEv->GetHeader()->GetTriggerCluster();
1032 fHeader->fOffTriggers = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
1033 fHeader->fFiredTriggers = fEsdEv->GetFiredTriggerClasses();
1035 AliCentrality *cent = InputEvent()->GetCentrality();
1036 fHeader->fV0Cent = cent->GetCentralityPercentileUnchecked("V0M");
1037 fHeader->fCl1Cent = cent->GetCentralityPercentileUnchecked("CL1");
1038 fHeader->fTrCent = cent->GetCentralityPercentileUnchecked("TRK");
1039 fHeader->fCqual = cent->GetQuality();
1042 TString trgclasses(fHeader->fFiredTriggers);
1043 for (Int_t j = 0; j<fTrClassNamesArr->GetEntries(); ++j) {
1044 const char *name = fTrClassNamesArr->At(j)->GetName();
1045 if (trgclasses.Contains(name))
1046 val += TMath::Power(2,j);
1048 fHeader->fTcls = (UInt_t)val;
1051 am->LoadBranch("vertices");
1052 AliAODVertex *pv = fAodEv->GetPrimaryVertex();
1053 FillVertex(fPrimVert, pv);
1054 AliAODVertex *sv = fAodEv->GetPrimaryVertexSPD();
1055 FillVertex(fSpdVert, sv);
1057 am->LoadBranch("PrimaryVertex.");
1058 const AliESDVertex *pv = fEsdEv->GetPrimaryVertexTracks();
1059 FillVertex(fPrimVert, pv);
1060 am->LoadBranch("SPDVertex.");
1061 const AliESDVertex *sv = fEsdEv->GetPrimaryVertexSPD();
1062 FillVertex(fSpdVert, sv);
1063 am->LoadBranch("TPCVertex.");
1064 const AliESDVertex *tv = fEsdEv->GetPrimaryVertexTPC();
1065 FillVertex(fTpcVert, tv);
1068 TObjArray *clusters = fEsdClusters;
1070 clusters = fAodClusters;
1075 Int_t nclus = clusters->GetEntries();
1076 for(Int_t i=0,ncl=0; i<nclus; ++i) {
1077 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1080 if (!clus->IsEMCAL())
1082 if (clus->E()<fMinE)
1085 AliStaCluster *cl = static_cast<AliStaCluster*>(fClusters->New(ncl++));
1088 clus->GetPosition(pos);
1090 cl->fR = vpos.Perp();
1091 cl->fEta = vpos.Eta();
1092 cl->fPhi = vpos.Phi();
1093 cl->fN = clus->GetNCells();
1094 cl->fN1 = GetNCells(clus,0.100);
1095 cl->fN3 = GetNCells(clus,0.300);
1097 Double_t emax = GetMaxCellEnergy(clus, id);
1100 cl->fDbc = clus->GetDistanceToBadChannel();;
1101 cl->fDisp = clus->GetDispersion();
1102 cl->fM20 = clus->GetM20();
1103 cl->fM02 = clus->GetM02();
1104 cl->fEcc = clus->Chi2(); //eccentricity
1105 cl->fSig = clus->GetTOF(); //sigma
1106 cl->fTrDz = fClusProps[i].fTrDz;
1107 cl->fTrDr = fClusProps[i].fTrDr;;
1108 cl->fTrEp = fClusProps[i].fTrEp;;
1109 cl->fTrIso = fClusProps[i].fTrIso;
1110 cl->fCeIso = fClusProps[i].fCellIso;
1115 //________________________________________________________________________
1116 void AliAnalysisTaskEMCALPi0PbPb::FillPionHists()
1118 // Fill histograms related to pions.
1120 Double_t vertex[3] = {0};
1121 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1123 TObjArray *clusters = fEsdClusters;
1125 clusters = fAodClusters;
1128 TLorentzVector clusterVec1;
1129 TLorentzVector clusterVec2;
1130 TLorentzVector pionVec;
1132 Int_t nclus = clusters->GetEntries();
1133 for (Int_t i = 0; i<nclus; ++i) {
1134 AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i));
1137 if (!clus1->IsEMCAL())
1139 if (clus1->E()<fMinE)
1141 if (clus1->GetNCells()<fNminCells)
1143 if (GetMaxCellEnergy(clus1)/clus1->E()<fMinErat)
1145 if (clus1->Chi2()<fMinEcc) // eccentricity cut
1147 clus1->GetMomentum(clusterVec1,vertex);
1148 for (Int_t j = i+1; j<nclus; ++j) {
1149 AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j));
1152 if (!clus2->IsEMCAL())
1154 if (clus2->E()<fMinE)
1156 if (clus2->GetNCells()<fNminCells)
1158 if (GetMaxCellEnergy(clus2)/clus2->E()<fMinErat)
1160 if (clus2->Chi2()<fMinEcc) // eccentricity cut
1162 clus2->GetMomentum(clusterVec2,vertex);
1163 pionVec = clusterVec1 + clusterVec2;
1164 Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E();
1165 Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect());
1166 if (pionZgg < fAsymMax) {
1167 fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi());
1168 fHPionMggPt->Fill(pionVec.M(),pionVec.Pt());
1169 fHPionMggAsym->Fill(pionVec.M(),pionZgg);
1170 fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle);
1171 Int_t bin = fPtRanges->FindBin(pionVec.Pt());
1172 fHPionInvMasses[bin]->Fill(pionVec.M());
1179 //________________________________________________________________________
1180 void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists()
1182 // Fill other histograms.
1185 //__________________________________________________________________________________________________
1186 void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliESDVertex *esdv)
1188 // Fill vertex from ESD vertex info.
1190 v->fVx = esdv->GetX();
1191 v->fVy = esdv->GetY();
1192 v->fVz = esdv->GetZ();
1193 v->fVc = esdv->GetNContributors();
1194 v->fDisp = esdv->GetDispersion();
1195 v->fZres = esdv->GetZRes();
1196 v->fChi2 = esdv->GetChi2();
1197 v->fSt = esdv->GetStatus();
1198 v->fIs3D = esdv->IsFromVertexer3D();
1199 v->fIsZ = esdv->IsFromVertexerZ();
1202 //__________________________________________________________________________________________________
1203 void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliAODVertex *aodv)
1205 // Fill vertex from AOD vertex info.
1207 v->fVx = aodv->GetX();
1208 v->fVy = aodv->GetY();
1209 v->fVz = aodv->GetZ();
1210 v->fVc = aodv->GetNContributors();
1211 v->fChi2 = aodv->GetChi2();
1214 //________________________________________________________________________
1215 Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
1217 // Compute isolation based on cell content.
1219 AliVCaloCells *cells = fEsdCells;
1225 Double_t cellIsolation = 0;
1226 Double_t rad2 = radius*radius;
1227 Int_t ncells = cells->GetNumberOfCells();
1228 for (Int_t i = 0; i<ncells; ++i) {
1229 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
1230 Double_t cellE = cells->GetAmplitude(i);
1232 fGeom->EtaPhiFromIndex(absID,eta,phi);
1233 Double_t phidiff = TVector2::Phi_0_2pi(phi-cPhi);
1234 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1237 cellIsolation += cellE;
1239 return cellIsolation;
1242 //________________________________________________________________________
1243 Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(AliVCluster *cluster, Short_t &id) const
1245 // Get maximum energy of attached cell.
1249 Int_t ncells = cluster->GetNCells();
1251 for (Int_t i=0; i<ncells; i++) {
1252 Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1255 id = cluster->GetCellAbsId(i);
1259 for (Int_t i=0; i<ncells; i++) {
1260 Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1263 id = cluster->GetCellAbsId(i);
1269 //________________________________________________________________________
1270 void AliAnalysisTaskEMCALPi0PbPb::GetSigma(AliVCluster *c, Double_t& sigmaMax, Double_t &sigmaMin) const
1272 // Calculate the (E) weighted variance along the longer (eigen) axis.
1274 sigmaMax = 0; // cluster variance along its longer axis
1275 sigmaMin = 0; // cluster variance along its shorter axis
1276 Double_t Ec = c->E(); // cluster energy
1279 Double_t Xc = 0 ; // cluster first moment along X
1280 Double_t Yc = 0 ; // cluster first moment along Y
1281 Double_t Sxx = 0 ; // cluster second central moment along X (variance_X^2)
1282 Double_t Sxy = 0 ; // cluster second central moment along Y (variance_Y^2)
1283 Double_t Syy = 0 ; // cluster covariance^2
1285 AliVCaloCells *cells = fEsdCells;
1292 Int_t ncells = c->GetNCells();
1297 for(Int_t j=0; j<ncells; ++j) {
1298 Int_t id = TMath::Abs(c->GetCellAbsId(j));
1299 fGeom->GetGlobal(id,pos);
1300 Double_t cellen = cells->GetCellAmplitude(id);
1301 Xc += cellen*pos.X();
1302 Yc += cellen*pos.Y();
1303 Sxx += cellen*pos.X()*pos.X();
1304 Syy += cellen*pos.Y()*pos.Y();
1305 Sxy += cellen*pos.X()*pos.Y();
1315 Sxx = TMath::Abs(Sxx);
1316 Syy = TMath::Abs(Syy);
1317 sigmaMax = (Sxx + Syy + TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
1318 sigmaMax = TMath::Sqrt(TMath::Abs(sigmaMax));
1319 sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
1320 sigmaMin = TMath::Sqrt(TMath::Abs(sigmaMin));
1323 //________________________________________________________________________
1324 Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(AliVCluster *c, Double_t emin) const
1326 // Calculate number of attached cells above emin.
1328 AliVCaloCells *cells = fEsdCells;
1335 Int_t ncells = c->GetNCells();
1336 for(Int_t j=0; j<ncells; ++j) {
1337 Int_t id = TMath::Abs(c->GetCellAbsId(j));
1338 Double_t cellen = cells->GetCellAmplitude(id);
1345 //________________________________________________________________________
1346 Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
1348 // Compute isolation based on tracks.
1350 Double_t trkIsolation = 0;
1351 Double_t rad2 = radius*radius;
1352 Int_t ntrks = fSelTracks->GetEntries();
1353 for(Int_t j = 0; j<ntrks; ++j) {
1354 AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
1357 Float_t eta = track->Eta();
1358 Float_t phi = track->Phi();
1359 Double_t phidiff = TVector2::Phi_0_2pi(phi-cPhi);
1360 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1363 trkIsolation += track->Pt();
1365 return trkIsolation;