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 "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"
  32 #include "AliLog.h"
  33 #include "AliMagF.h"
  34 #include "AliTrackerBase.h"
  35
  36 ClassImp(AliAnalysisTaskEMCALPi0PbPb)
  37
  38 //________________________________________________________________________
  39 AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name)
  40   : AliAnalysisTaskSE(name),
  41     fCentVar("V0M"),
  42     fCentFrom(0),
  43     fCentTo(100),
  44     fVtxZMin(-10),
  45     fVtxZMax(+10),
  46     fUseQualFlag(1),
  47     fClusName(),
  48     fDoNtuple(0),
  49     fDoAfterburner(0),
  50     fAsymMax(1),
  51     fNminCells(1),
  52     fMinE(0.100),
  53     fMinErat(0),
  54     fMinEcc(0),
  55     fGeoName("EMCAL_FIRSTYEARV1"),
  56     fMinNClustPerTrack(50),
  57     fMinPtPerTrack(1.0),
  58     fIsoDist(0.2),
  59     fTrClassNames(""),
  60     fTrCuts(0),
  61     fDoTrackMatWithGeom(0),
  62     fDoConstrain(0),
  63     fNEvs(0),
  64     fGeom(0),
  65     fReco(0),
  66     fOutput(0),
  67     fTrClassNamesArr(0),
  68     fEsdEv(0),
  69     fAodEv(0),
  70     fRecPoints(0),
  71     fEsdClusters(0),
  72     fEsdCells(0),
  73     fAodClusters(0),
  74     fAodCells(0),
  75     fPtRanges(0),
  76     fSelTracks(0),
  77     fNtuple(0),
  78     fHeader(0),
  79     fPrimVert(0),
  80     fSpdVert(0),
  81     fTpcVert(0),
  82     fClusters(0),
  83     fHCuts(0x0),
  84     fHVertexZ(0x0),
  85     fHVertexZ2(0x0),
  86     fHCent(0x0),
  87     fHCentQual(0x0),
  88     fHTclsBeforeCuts(0x0),
  89     fHTclsAfterCuts(0x0),
  90     fHColuRow(0x0),
  91     fHColuRowE(0x0),
  92     fHCellMult(0x0),
  93     fHCellE(0x0),
  94     fHCellH(0x0),
  95     fHCellM(0x0),
  96     fHCellM2(0x0),
  97     fHCellFreqNoCut(0x0),
  98     fHCellFreqCut100M(0x0),
  99     fHCellFreqCut300M(0x0),
 100     fHCellFreqE(0x0),
 101     fHCellCheckE(0x0),
 102     fHClustEccentricity(0),
 103     fHClustEtaPhi(0x0),
 104     fHClustEnergyPt(0x0),
 105     fHClustEnergySigma(0x0),
 106     fHClustSigmaSigma(0x0),
 107     fHClustNCellEnergyRatio(0x0),
 108     fHPionEtaPhi(0x0),
 109     fHPionMggPt(0x0),
 110     fHPionMggAsym(0x0),
 111     fHPionMggDgg(0x0)
 112 {
 113   // Constructor.
 114
 115   if (!name)
 116     return;
 117   SetName(name);
 118   DefineInput(0, TChain::Class());
 119   DefineOutput(1, TList::Class());
 120   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,EMCALCells.,Tracks "
 121                "AOD:header,vertices,emcalCells,tracks";
 122 }
 123
 124 //________________________________________________________________________
 125 AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb()
 126 {
 127   // Destructor.
 128
 129   if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 130     delete fOutput; fOutput = 0;
 131   }
 132   delete fPtRanges; fPtRanges = 0;
 133   delete fGeom; fGeom = 0;
 134   delete fReco; fReco = 0;
 135   delete fTrClassNamesArr;
 136   delete fSelTracks;
 137   delete [] fHColuRow;
 138   delete [] fHColuRowE;
 139   delete [] fHCellMult;
 140   delete [] fHCellFreqNoCut;
 141   delete [] fHCellFreqCut100M;
 142   delete [] fHCellFreqCut300M;
 143 }
 144
 145 //________________________________________________________________________
 146 void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects()
 147 {
 148   // Create user objects here.
 149
 150   fGeom = new AliEMCALGeoUtils(fGeoName,"EMCAL");
 151   fReco = new AliEMCALRecoUtils();
 152   fTrClassNamesArr = fTrClassNames.Tokenize(" ");
 153   fOutput = new TList();
 154   fOutput->SetOwner();
 155   fSelTracks = new TObjArray;
 156
 157   if (fDoNtuple) {
 158     TFile *f = OpenFile(1);
 159     if (f) {
 160       f->SetCompressionLevel(2);
 161       fNtuple = new TTree(Form("tree%.0fto%.0f",fCentFrom,fCentTo), "StandaloneTree");
 162       fNtuple->SetDirectory(f);
 163       fNtuple->SetAutoFlush(-1024*1024*1024);
 164       fNtuple->SetAutoSave(-1024*1024*1024);
 165       fHeader = new AliStaHeader;
 166       fNtuple->Branch("header", &fHeader, 16*1024, 99);
 167       fPrimVert = new AliStaVertex;
 168       fNtuple->Branch("primv", &fPrimVert, 16*1024, 99);
 169       fSpdVert = new AliStaVertex;
 170       fNtuple->Branch("spdv", &fSpdVert, 16*1024, 99);
 171       fTpcVert = new AliStaVertex;
 172       fNtuple->Branch("tpcv", &fTpcVert, 16*1024, 99);
 173       if (TClass::GetClass("AliStaCluster"))
 174         TClass::GetClass("AliStaCluster")->IgnoreTObjectStreamer();
 175       fClusters = new TClonesArray("AliStaCluster",1000);
 176       fNtuple->Branch("clusters", &fClusters, 8*16*1024, 99);
 177     }
 178   }
 179
 180   AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
 181   Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-0.25;
 182   Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+0.25;
 183
 184   // histograms
 185   TH1::SetDefaultSumw2(kTRUE);
 186   TH2::SetDefaultSumw2(kTRUE);
 187   fHCuts = new TH1F("hEventCuts","",5,0.5,5.5);
 188   fHCuts->GetXaxis()->SetBinLabel(1,"All");
 189   fHCuts->GetXaxis()->SetBinLabel(2,"PS");
 190   fHCuts->GetXaxis()->SetBinLabel(3,Form("%s: %.0f-%.0f",fCentVar.Data(),fCentFrom,fCentTo));
 191   fHCuts->GetXaxis()->SetBinLabel(4,"QFlag");
 192   fHCuts->GetXaxis()->SetBinLabel(5,Form("zvtx: %.0f-%.0f",fVtxZMin,fVtxZMax));
 193   fOutput->Add(fHCuts);
 194   fHVertexZ = new TH1F("hVertexZBeforeCut","",100,-25,25);
 195   fHVertexZ->SetXTitle("z [cm]");
 196   fOutput->Add(fHVertexZ);
 197   fHVertexZ2 = new TH1F("hVertexZAfterCut","",100,-25,25);
 198   fHVertexZ2->SetXTitle("z [cm]");
 199   fOutput->Add(fHVertexZ2);
 200   fHCent = new TH1F("hCentBeforeCut","",101,-1,100);
 201   fHCent->SetXTitle(fCentVar.Data());
 202   fOutput->Add(fHCent);
 203   fHCentQual = new TH1F("hCentAfterCut","",101,-1,100);
 204   fHCentQual->SetXTitle(fCentVar.Data());
 205   fOutput->Add(fHCentQual);
 206   fHTclsBeforeCuts = new TH1F("hTclsBeforeCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
 207   fHTclsAfterCuts = new TH1F("hTclsAfterCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
 208   for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
 209     const char *name = fTrClassNamesArr->At(i)->GetName();
 210     fHTclsBeforeCuts->GetXaxis()->SetBinLabel(1+i,name);
 211     fHTclsAfterCuts->GetXaxis()->SetBinLabel(1+i,name);
 212   }
 213   fOutput->Add(fHTclsBeforeCuts);
 214   fOutput->Add(fHTclsAfterCuts);
 215
 216   // histograms for cells
 217   Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
 218   fHColuRow   = new TH2*[nsm];
 219   fHColuRowE  = new TH2*[nsm];
 220   fHCellMult  = new TH1*[nsm];
 221   for (Int_t i = 0; i<nsm; ++i) {
 222     fHColuRow[i] = new TH2F(Form("hColRow_Mod%d",i),"",48,0,48,24,0.,24);
 223     fHColuRow[i]->SetTitle(Form("Module %d: Occupancy", i));
 224     fHColuRow[i]->SetXTitle("col (i#eta)");
 225     fHColuRow[i]->SetYTitle("row (i#phi)");
 226     fHColuRowE[i] = new TH2F(Form("hColRowE_Mod%d", i),"",48,0,48,24,0,24);
 227     fHColuRowE[i]->SetTitle(Form("Module %d: Cell energy",i));
 228     fHColuRowE[i]->SetXTitle("col (i#eta)");
 229     fHColuRowE[i]->SetYTitle("row (i#phi)");
 230     fHCellMult[i] = new TH1F(Form("hCellMult_Mod%d",i),"",1000,0,1000);
 231     fHCellMult[i]->SetTitle(Form("Module %d: Cell multiplicity",i));
 232     fHCellMult[i]->SetXTitle("# of cells");
 233     fOutput->Add(fHColuRow[i]);
 234     fOutput->Add(fHColuRowE[i]);
 235     fOutput->Add(fHCellMult[i]);
 236   }
 237   fHCellE = new TH1F("hCellE","",250,0.,25.);
 238   fHCellE->SetXTitle("E_{cell} [GeV]");
 239   fOutput->Add(fHCellE);
 240   fHCellH = new TH1F ("hCellHighestE","",250,0.,25.);
 241   fHCellH->SetXTitle("E^{max}_{cell} [GeV]");
 242   fOutput->Add(fHCellH);
 243   fHCellM = new TH1F ("hCellMeanEperHitCell","",250,0.,2.5);
 244   fHCellM->SetXTitle("#LT E_{cell}#GT [GeV]");
 245   fOutput->Add(fHCellM);
 246   fHCellM2 = new TH1F ("hCellMeanEperAllCells","",250,0.,1);
 247   fHCellM2->SetXTitle("1/N_{cells} #Sigma E_{cell} [GeV]");
 248   fOutput->Add(fHCellM2);
 249
 250   fHCellFreqNoCut   = new TH1*[nsm];
 251   fHCellFreqCut100M = new TH1*[nsm];
 252   fHCellFreqCut300M = new TH1*[nsm];
 253   fHCellFreqE       = new TH1*[nsm];
 254   for (Int_t i = 0; i<nsm; ++i){
 255     Double_t lbin = i*24*48-0.5;
 256     Double_t hbin = lbin+24*48;
 257     fHCellFreqNoCut[i]   = new TH1F(Form("hCellFreqNoCut_SM%d",i),
 258                                     Form("Frequency SM%d (no cut);id;#",i),  1152, lbin, hbin);
 259     fHCellFreqCut100M[i] = new TH1F(Form("hCellFreqCut100M_SM%d",i),
 260                                     Form("Frequency SM%d (>0.1GeV);id;#",i), 1152, lbin, hbin);
 261     fHCellFreqCut300M[i] = new TH1F(Form("hCellFreqCut300M_SM%d",i),
 262                                     Form("Frequency SM%d (>0.3GeV);id;#",i), 1152, lbin, hbin);
 263     fHCellFreqE[i]       = new TH1F(Form("hCellFreqE_SM%d",i),
 264                                     Form("Frequency SM%d (E weighted);id;#",i), 1152, lbin, hbin);
 265     fOutput->Add(fHCellFreqNoCut[i]);
 266     fOutput->Add(fHCellFreqCut100M[i]);
 267     fOutput->Add(fHCellFreqCut300M[i]);
 268     fOutput->Add(fHCellFreqE[i]);
 269   }
 270   if (1) {
 271     fHCellCheckE = new TH1*[24*48*nsm];
 272     memset(fHCellCheckE,0,24*48*nsm*sizeof(TH1*));
 273     Int_t tcs[1] = {4102};
 274     for (UInt_t i = 0; i<sizeof(tcs)/sizeof(Int_t); ++i){
 275       Int_t c=tcs[i];
 276       if (c<24*48*nsm) {
 277         fHCellCheckE[i] = new TH1F(Form("hCellE_id%d",c), Form("Cell %d;E [GeV/c];#",c), 500, 0, 8);
 278         fOutput->Add(fHCellCheckE[i]);
 279       }
 280     }
 281   }
 282
 283   // histograms for clusters
 284   fHClustEccentricity = new TH1F("hClustEccentricity","",100,-0.1,1.1);
 285   fHClustEccentricity->SetXTitle("#epsilon_{C}");
 286   fOutput->Add(fHClustEccentricity);
 287   fHClustEtaPhi = new TH2F("hClustEtaPhi","",500,-0.8,0.8,100*nsm,phimin,phimax);
 288   fHClustEtaPhi->SetXTitle("#eta");
 289   fHClustEtaPhi->SetYTitle("#varphi");
 290   fOutput->Add(fHClustEtaPhi);
 291   fHClustEnergyPt = new TH2F("hClustEnergyPt","",250,0,50,250,0,50);
 292   fHClustEnergyPt->SetXTitle("E [GeV]");
 293   fHClustEnergyPt->SetYTitle("p_{T} [GeV/c]");
 294   fOutput->Add(fHClustEnergyPt);
 295   fHClustEnergySigma = new TH2F("hClustEnergySigma","",100,0,100,500,0,50);
 296   fHClustEnergySigma->SetXTitle("E_{C} * #sigma_{max} [GeV*cm]");
 297   fHClustEnergySigma->SetYTitle("E_{C} [GeV]");
 298   fOutput->Add(fHClustEnergySigma);
 299   fHClustSigmaSigma = new TH2F("hClustSigmaSigma","",500,0,50,500,0,50);
 300   fHClustSigmaSigma->SetXTitle("#lambda_{0} [cm]");
 301   fHClustSigmaSigma->SetYTitle("#sigma_{max} [cm]");
 302   fOutput->Add(fHClustSigmaSigma);
 303   fHClustNCellEnergyRatio = new TH2F("hClustNCellEnergyRatio","",27,-0.5,26.5,101,-0.05,1.05);
 304   fHClustNCellEnergyRatio->SetXTitle("N_{cells}");
 305   fHClustNCellEnergyRatio->SetYTitle("E^{max}_{cell}/E_{clus}");
 306   fOutput->Add(fHClustNCellEnergyRatio);
 307
 308   // histograms for pion candidates
 309   fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100*nsm,phimin,phimax);
 310   fHPionEtaPhi->SetXTitle("#eta_{#gamma#gamma}");
 311   fHPionEtaPhi->SetYTitle("#varphi_{#gamma#gamma}");
 312   fOutput->Add(fHPionEtaPhi);
 313   fHPionMggPt = new TH2F("hPionMggPt","",1000,0,2,100,0,20.0);
 314   fHPionMggPt->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 315   fHPionMggPt->SetYTitle("p_{T}^{#gamma#gamma} [GeV/c]");
 316   fOutput->Add(fHPionMggPt);
 317   fHPionMggAsym = new TH2F("hPionMggAsym","",1000,0,2,100,0,1);
 318   fHPionMggAsym->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 319   fHPionMggAsym->SetYTitle("Z_{#gamma#gamma} [GeV]");
 320   fOutput->Add(fHPionMggAsym);
 321   fHPionMggDgg = new TH2F("hPionMggDgg","",1000,0,2,100,0,10);
 322   fHPionMggDgg->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 323   fHPionMggDgg->SetYTitle("opening angle [grad]");
 324   fOutput->Add(fHPionMggDgg);
 325   const Int_t nbins = 20;
 326   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};
 327   fPtRanges = new TAxis(nbins-1,xbins);
 328   for (Int_t i = 0; i<=nbins; ++i) {
 329     fHPionInvMasses[i] = new TH1F(Form("hPionInvMass%d",i),"",1000,0,2);
 330     fHPionInvMasses[i]->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 331     if (i==0)
 332       fHPionInvMasses[i]->SetTitle(Form("0 < p_{T}^{#gamma#gamma} <%.1f",xbins[0]));
 333     else if (i==nbins)
 334       fHPionInvMasses[i]->SetTitle(Form("p_{T}^{#gamma#gamma} > 50"));
 335     else
 336       fHPionInvMasses[i]->SetTitle(Form("%.1f < p_{T}^{#gamma#gamma} <%.1f",xbins[i-1],xbins[i]));
 337     fOutput->Add(fHPionInvMasses[i]);
 338   }
 339
 340   PostData(1, fOutput);
 341 }
 342
 343 //________________________________________________________________________
 344 void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *)
 345 {
 346   // Called for each event.
 347
 348   if (!InputEvent())
 349     return;
 350
 351   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 352   fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent());
 353   if (fEsdEv) {
 354     am->LoadBranch("AliESDRun.");
 355     am->LoadBranch("AliESDHeader.");
 356   } else {
 357     fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
 358     if (!fAodEv) {
 359       AliFatal("Neither ESD nor AOD event found");
 360       return;
 361     }
 362     am->LoadBranch("header");
 363   }
 364
 365   Int_t cut = 1;
 366   fHCuts->Fill(cut++);
 367
 368   TString trgclasses;
 369   AliESDHeader *h = dynamic_cast<AliESDHeader*>(InputEvent()->GetHeader());
 370   if (h) {
 371     trgclasses = fEsdEv->GetFiredTriggerClasses();
 372   } else {
 373     AliAODHeader *h2 = dynamic_cast<AliAODHeader*>(InputEvent()->GetHeader());
 374     if (h2)
 375       trgclasses = h2->GetFiredTriggerClasses();
 376   }
 377   for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
 378     const char *name = fTrClassNamesArr->At(i)->GetName();
 379     if (trgclasses.Contains(name))
 380       fHTclsBeforeCuts->Fill(1+i);
 381   }
 382
 383   UInt_t res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 384   if (res==0)
 385     return;
 386
 387   if (fHCuts->GetBinContent(2)==0) {
 388     if (fDoTrackMatWithGeom && !AliGeomManager::GetGeometry()) { // get geometry
 389       AliWarning("Accessing geometry from OCDB, this is not very efficient!");
 390       AliCDBManager *cdb = AliCDBManager::Instance();
 391       if (!cdb->IsDefaultStorageSet())
 392         cdb->SetDefaultStorage("raw://");
 393       Int_t runno = InputEvent()->GetRunNumber();
 394       if (runno != cdb->GetRun())
 395         cdb->SetRun(runno);
 396       AliGeomManager::LoadGeometry();
 397     }
 398
 399     if (!fGeom->GetMatrixForSuperModule(0)) { // set misalignment matrices (stored in first event)
 400       if (fEsdEv) {
 401         for (Int_t i=0; i<fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++i)
 402           fGeom->SetMisalMatrix(fEsdEv->GetESDRun()->GetEMCALMatrix(i),i);
 403       } else {
 404         for (Int_t i=0; i<fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++i)
 405           fGeom->SetMisalMatrix(fAodEv->GetHeader()->GetEMCALMatrix(i),i);
 406       }
 407     }
 408
 409     if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
 410       if (fEsdEv) {
 411         const AliESDRun *erun = fEsdEv->GetESDRun();
 412         AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(),
 413                                                  erun->GetCurrentDip(),
 414                                                  AliMagF::kConvLHC,
 415                                                  kFALSE,
 416                                                  erun->GetBeamEnergy(),
 417                                                  erun->GetBeamType());
 418         TGeoGlobalMagField::Instance()->SetField(field);
 419       } else {
 420         Double_t pol = -1; //polarity
 421         Double_t be = -1;  //beam energy
 422         AliMagF::BeamType_t btype = AliMagF::kBeamTypepp;
 423         Int_t runno = fAodEv->GetRunNumber();
 424         if (runno>=136851 && runno<138275) {
 425           pol = -1;
 426           be = 2760;
 427           btype = AliMagF::kBeamTypeAA;
 428         } else if (runno>=138275 && runno<=139517) {
 429           pol = +1;
 430           be = 2760;
 431           btype = AliMagF::kBeamTypeAA;
 432         } else {
 433           AliError(Form("Do not know the bfield parameters for run %d! Using defaults!!!", runno));
 434         }
 435         TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", pol, pol, AliMagF::k5kG, btype, be));
 436       }
 437     }
 438   }
 439
 440   fHCuts->Fill(cut++);
 441
 442   const AliCentrality *centP = InputEvent()->GetCentrality();
 443   Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar);
 444   fHCent->Fill(cent);
 445   if (cent<fCentFrom||cent>fCentTo)
 446     return;
 447
 448   fHCuts->Fill(cut++);
 449
 450   if (fUseQualFlag) {
 451     if (centP->GetQuality()>0)
 452       return;
 453   }
 454
 455   fHCentQual->Fill(cent);
 456   fHCuts->Fill(cut++);
 457
 458   if (fEsdEv) {
 459     am->LoadBranch("PrimaryVertex.");
 460     am->LoadBranch("SPDVertex.");
 461     am->LoadBranch("TPCVertex.");
 462   } else {
 463     fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
 464     am->LoadBranch("vertices");
 465     if (!fAodEv) return;
 466   }
 467
 468   const AliVVertex *vertex = InputEvent()->GetPrimaryVertex();
 469   if (!vertex)
 470     return;
 471
 472   fHVertexZ->Fill(vertex->GetZ());
 473
 474   if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax)
 475     return;
 476
 477   fHCuts->Fill(cut++);
 478   fHVertexZ2->Fill(vertex->GetZ());
 479
 480   // count number of accepted events
 481   ++fNEvs;
 482
 483   for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
 484     const char *name = fTrClassNamesArr->At(i)->GetName();
 485     if (trgclasses.Contains(name))
 486       fHTclsAfterCuts->Fill(1+i);
 487   }
 488
 489   fRecPoints   = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
 490   fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set of if clusters are attached
 491   fEsdCells    = 0; // will be set if ESD input used
 492   fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set of if clusters are attached
 493                     //             or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode
 494   fAodCells    = 0; // will be set if AOD input used
 495
 496   // deal with special output from AliAnalysisTaskEMCALClusterizeFast first
 497   Bool_t clusattached = 0;
 498   Bool_t recalibrated = 0;
 499   if (1 && !fClusName.IsNull()) {
 500     AliAnalysisTaskEMCALClusterizeFast *cltask = 0;
 501     TObjArray *ts = am->GetTasks();
 502     cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName));
 503     if (cltask && cltask->GetClusters()) {
 504       fRecPoints = const_cast<TObjArray*>(cltask->GetClusters());
 505       clusattached = cltask->GetAttachClusters();
 506       if (cltask->GetCalibData()!=0)
 507         recalibrated = kTRUE;
 508     }
 509   }
 510   if (1 && AODEvent() && !fClusName.IsNull()) {
 511     TList *l = AODEvent()->GetList();
 512     TClonesArray *clus = 0;
 513     if (l) {
 514       clus = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
 515       fAodClusters = clus;
 516     }
 517   }
 518
 519   if (fEsdEv) { // ESD input mode
 520     if (1 && (!fRecPoints||clusattached)) {
 521       if (!clusattached)
 522         am->LoadBranch("CaloClusters");
 523       TList *l = fEsdEv->GetList();
 524       TClonesArray *clus = 0;
 525       if (l) {
 526         clus = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
 527         fEsdClusters = clus;
 528       }
 529     }
 530     if (1) {
 531       if (!recalibrated)
 532         am->LoadBranch("EMCALCells.");
 533       fEsdCells = fEsdEv->GetEMCALCells();
 534     }
 535   } else if (fAodEv) { // AOD input mode
 536     if (1 && (!fAodClusters || clusattached)) {
 537       if (!clusattached)
 538         am->LoadBranch("caloClusters");
 539       TList *l = fAodEv->GetList();
 540       TClonesArray *clus = 0;
 541       if (l) {
 542         clus = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
 543         fAodClusters = clus;
 544       }
 545     }
 546     if (1) {
 547       if (!recalibrated)
 548         am->LoadBranch("emcalCells");
 549       fAodCells = fAodEv->GetEMCALCells();
 550     }
 551   } else {
 552     AliFatal("Impossible to not have either pointer to ESD or AOD event");
 553   }
 554
 555   if (1) {
 556     AliDebug(2,Form("fRecPoints   set: %p", fRecPoints));
 557     AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters));
 558     AliDebug(2,Form("fEsdCells    set: %p", fEsdCells));
 559     AliDebug(2,Form("fAodClusters set: %p", fAodClusters));
 560     AliDebug(2,Form("fAodCells    set: %p", fAodCells));
 561   }
 562
 563   if (fDoAfterburner)
 564     ClusterAfterburner();
 565
 566   CalcTracks();
 567   CalcClusterProps();
 568
 569   FillCellHists();
 570   FillClusHists();
 571   FillPionHists();
 572   FillOtherHists();
 573   FillNtuple();
 574
 575   PostData(1, fOutput);
 576 }
 577
 578 //________________________________________________________________________
 579 void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *)
 580 {
 581   // Terminate called at the end of analysis.
 582
 583   if (fNtuple) {
 584     TFile *f = OpenFile(1);
 585     if (f)
 586       fNtuple->Write();
 587   }
 588
 589   AliInfo(Form("%s: Accepted %lld events", GetName(), fNEvs));
 590 }
 591
 592 //________________________________________________________________________
 593 void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
 594 {
 595   // Calculate track properties.
 596
 597   fSelTracks->Clear();
 598
 599   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 600   TClonesArray *tracks = 0;
 601   if (fEsdEv) {
 602     am->LoadBranch("Tracks");
 603     TList *l = fEsdEv->GetList();
 604     tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
 605   } else {
 606     am->LoadBranch("tracks");
 607     TList *l = fAodEv->GetList();
 608     tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
 609   }
 610
 611   if (!tracks)
 612     return;
 613
 614   AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
 615   Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-0.25;
 616   Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+0.25;
 617
 618   if (fEsdEv) {
 619     if (fDoConstrain)
 620       fSelTracks->SetOwner(kTRUE);
 621     am->LoadBranch("PrimaryVertex.");
 622     const AliESDVertex *vtx = fEsdEv->GetPrimaryVertexSPD();
 623     am->LoadBranch("SPDVertex.");
 624     const AliESDVertex *vtxSPD = fEsdEv->GetPrimaryVertexSPD();
 625     am->LoadBranch("Tracks");
 626     const Int_t Ntracks = tracks->GetEntries();
 627     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
 628       AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
 629       if (!track) {
 630         AliWarning(Form("Could not receive track %d\n", iTracks));
 631         continue;
 632       }
 633       if (fTrCuts && !fTrCuts->IsSelected(track))
 634         continue;
 635       Double_t eta = track->Eta();
 636       if (eta<-1||eta>1)
 637         continue;
 638       if (track->Phi()<phimin||track->Phi()>phimax)
 639         continue;
 640       if(track->GetTPCNcls()<fMinNClustPerTrack)
 641         continue;
 642
 643       if (!fDoConstrain) {
 644         fSelTracks->Add(track);
 645         continue;
 646       }
 647
 648       AliESDtrack copyt(*track);
 649       Double_t bfield[3];
 650       copyt.GetBxByBz(bfield);
 651       AliExternalTrackParam tParam;
 652       Bool_t relate = copyt.RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&tParam);
 653       if (!relate)
 654         continue;
 655       copyt.Set(tParam.GetX(),tParam.GetAlpha(),tParam.GetParameter(),tParam.GetCovariance());
 656
 657       Double_t p[3]      = { 0. };
 658       copyt.GetPxPyPz(p);
 659       Double_t pos[3]    = { 0. };
 660       copyt.GetXYZ(pos);
 661       Double_t covTr[21] = { 0. };
 662       copyt.GetCovarianceXYZPxPyPz(covTr);
 663       Double_t pid[10]   = { 0. };
 664       copyt.GetESDpid(pid);
 665       AliAODTrack *aTrack = new AliAODTrack(copyt.GetID(),
 666                                             copyt.GetLabel(),
 667                                             p,
 668                                             kTRUE,
 669                                             pos,
 670                                             kFALSE,
 671                                             covTr,
 672                                             (Short_t)copyt.GetSign(),
 673                                             copyt.GetITSClusterMap(),
 674                                             pid,
 675                                             0,/*fPrimaryVertex,*/
 676                                             kTRUE, // check if this is right
 677                                             vtx->UsesTrack(copyt.GetID()));
 678       aTrack->SetTPCClusterMap(copyt.GetTPCClusterMap());
 679       aTrack->SetTPCSharedMap (copyt.GetTPCSharedMap());
 680       Float_t ndf = copyt.GetTPCNcls() + 1 - 5 ;
 681       if(ndf>0)
 682         aTrack->SetChi2perNDF(copyt.GetTPCchi2()/ndf);
 683       else
 684         aTrack->SetChi2perNDF(-1);
 685       aTrack->SetFlags(copyt.GetStatus());
 686       aTrack->SetTPCPointsF(copyt.GetTPCNclsF());
 687       fSelTracks->Add(aTrack);
 688     }
 689   } else {
 690     Int_t ntracks = tracks->GetEntries();
 691     for (Int_t i=0; i<ntracks; ++i) {
 692       AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
 693       if (!track)
 694         continue;
 695       Double_t eta = track->Eta();
 696       if (eta<-1||eta>1)
 697         continue;
 698       if (track->Phi()<phimin||track->Phi()>phimax)
 699         continue;
 700       if(track->GetTPCNcls()<fMinNClustPerTrack)
 701         continue;
 702
 703       if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL
 704         AliExternalTrackParam tParam(track);
 705         if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 438, 0.139, 1, kTRUE)) {
 706           Double_t trkPos[3];
 707           tParam.GetXYZ(trkPos);
 708           track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]);
 709         }
 710       }
 711       fSelTracks->Add(track);
 712     }
 713   }
 714 }
 715
 716 //________________________________________________________________________
 717 void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
 718 {
 719   // Calculate cluster properties
 720
 721   TObjArray *clusters = fEsdClusters;
 722   if (!clusters)
 723     clusters = fAodClusters;
 724   if (!clusters)
 725     return;
 726
 727   Int_t nclus = clusters->GetEntries();
 728   Int_t ntrks = fSelTracks->GetEntries();
 729
 730   for(Int_t i = 0; i<nclus; ++i) {
 731     fClusProps[i].Reset();
 732
 733     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
 734     if (!clus)
 735       continue;
 736     if (!clus->IsEMCAL())
 737       continue;
 738     if (clus->E()<fMinE)
 739       continue;
 740
 741     Float_t  clsPos[3] = {0};
 742     clus->GetPosition(clsPos);
 743     TVector3 clsVec(clsPos);
 744     Double_t vertex[3] = {0};
 745     InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 746     TLorentzVector clusterVec;
 747     clus->GetMomentum(clusterVec,vertex);
 748     Double_t clsEta = clusterVec.Eta();
 749
 750     Double_t mind2 = 1e10;
 751     for(Int_t j = 0; j<ntrks; ++j) {
 752       AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
 753       if (!track)
 754         continue;
 755       Double_t pt  = track->Pt();
 756       if (pt<fMinPtPerTrack)
 757         continue;
 758       if (TMath::Abs(clsEta-track->Eta())>0.5)
 759         continue;
 760
 761       Float_t tmpR=-1, tmpZ=-1;
 762       if (!fDoTrackMatWithGeom) {
 763
 764         AliExternalTrackParam *tParam = 0;
 765         if (fEsdEv) {
 766           AliESDtrack *esdTrack = static_cast<AliESDtrack*>(track);
 767           tParam = new AliExternalTrackParam(*esdTrack->GetTPCInnerParam());
 768         } else
 769           tParam = new AliExternalTrackParam(track);
 770
 771         Double_t bfield[3];
 772         track->GetBxByBz(bfield);
 773         TVector3 vec(clsPos);
 774         Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
 775         vec.RotateZ(-alpha);  //Rotate the cluster to the local extrapolation coordinate system
 776         tParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
 777         Bool_t ret = tParam->PropagateToBxByBz(vec.X(), bfield);
 778         if (!ret) {
 779           delete tParam;
 780           continue;
 781         }
 782         Double_t trkPos[3];
 783         tParam->GetXYZ(trkPos); //Get the extrapolated global position
 784         tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
 785         tmpZ = clsPos[2]-trkPos[2];
 786         delete tParam;
 787       } else {
 788         if (TMath::Abs(clsEta-track->Eta())>fIsoDist)
 789           continue;
 790         AliExternalTrackParam tParam(track);
 791         if (!fReco->ExtrapolateTrackToCluster(&tParam, clus, tmpR, tmpZ))
 792           continue;
 793       }
 794
 795       Double_t d2 = tmpR;
 796       if (mind2>d2) {
 797         mind2=d2;
 798         fClusProps[i].fTrIndex = j;
 799         fClusProps[i].fTrDz    = tmpZ;
 800         fClusProps[i].fTrDr    = TMath::Sqrt(tmpR*tmpR-tmpZ*tmpZ);
 801         fClusProps[i].fTrDist  = d2;
 802         fClusProps[i].fTrEp    = clus->E()/track->P();
 803       }
 804     }
 805
 806     if (0 && (fClusProps[i].fTrIndex>=0)) {
 807       cout << i << " " << fClusProps[i].fTrIndex << ": Dr " << fClusProps[i].fTrDr << " " << " Dz " << fClusProps[i].fTrDz << endl;
 808     }
 809
 810     fClusProps[i].fTrIso      = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
 811     fClusProps[i].fTrLowPtIso = 0;
 812     fClusProps[i].fCellIso    = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
 813   }
 814 }
 815
 816 //________________________________________________________________________
 817 void  AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner()
 818 {
 819   // Run custer reconstruction afterburner.
 820
 821   AliVCaloCells *cells = fEsdCells;
 822   if (!cells)
 823     cells = fAodCells;
 824
 825   if (!cells)
 826     return;
 827
 828   Int_t ncells = cells->GetNumberOfCells();
 829   if (ncells<=0)
 830     return;
 831
 832   Double_t cellMeanE = 0, cellSigE = 0;
 833   for (Int_t i = 0; i<ncells; ++i) {
 834     Double_t cellE = cells->GetAmplitude(i);
 835     cellMeanE += cellE;
 836     cellSigE += cellE*cellE;
 837   }
 838   cellMeanE /= ncells;
 839   cellSigE /= ncells;
 840   cellSigE -= cellMeanE*cellMeanE;
 841   if (cellSigE<0)
 842     cellSigE = 0;
 843   cellSigE = TMath::Sqrt(cellSigE / ncells);
 844
 845   Double_t subE = cellMeanE - 7*cellSigE;
 846   if (subE<0)
 847     return;
 848
 849   for (Int_t i = 0; i<ncells; ++i) {
 850     Short_t id=-1;
 851     Double_t amp=0,time=0;
 852     if (!cells->GetCell(i, id, amp, time))
 853       continue;
 854     amp -= cellMeanE;
 855     if (amp<0.001)
 856       amp = 0;
 857     cells->SetCell(i, id, amp, time);
 858   }
 859
 860   TObjArray *clusters = fEsdClusters;
 861   if (!clusters)
 862     clusters = fAodClusters;
 863   if (!clusters)
 864     return;
 865
 866   Int_t nclus = clusters->GetEntries();
 867   for (Int_t i = 0; i<nclus; ++i) {
 868     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
 869     if (!clus->IsEMCAL())
 870       continue;
 871     Int_t nc = clus->GetNCells();
 872     Double_t clusE = 0;
 873     UShort_t ids[100] = {0};
 874     Double_t fra[100] = {0};
 875     for (Int_t j = 0; j<nc; ++j) {
 876       Short_t id = TMath::Abs(clus->GetCellAbsId(j));
 877       Double_t cen = cells->GetCellAmplitude(id);
 878       clusE += cen;
 879       if (cen>0) {
 880         ids[nc] = id;
 881         ++nc;
 882       }
 883     }
 884     if (clusE<=0) {
 885       clusters->RemoveAt(i);
 886       continue;
 887     }
 888
 889     for (Int_t j = 0; j<nc; ++j) {
 890       Short_t id = ids[j];
 891       Double_t cen = cells->GetCellAmplitude(id);
 892       fra[j] = cen/clusE;
 893     }
 894     clus->SetE(clusE);
 895     AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus);
 896     if (aodclus) {
 897       aodclus->Clear("");
 898       aodclus->SetNCells(nc);
 899       aodclus->SetCellsAmplitudeFraction(fra);
 900       aodclus->SetCellsAbsId(ids);
 901     }
 902   }
 903   clusters->Compress();
 904 }
 905
 906 //________________________________________________________________________
 907 void AliAnalysisTaskEMCALPi0PbPb::FillCellHists()
 908 {
 909   // Fill histograms related to cell properties.
 910
 911   AliVCaloCells *cells = fEsdCells;
 912   if (!cells)
 913     cells = fAodCells;
 914
 915   if (!cells)
 916     return;
 917
 918   Int_t cellModCount[12] = {0};
 919   Double_t cellMaxE = 0;
 920   Double_t cellMeanE = 0;
 921   Int_t ncells = cells->GetNumberOfCells();
 922   if (ncells==0)
 923     return;
 924
 925   Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
 926
 927   for (Int_t i = 0; i<ncells; ++i) {
 928     Int_t absID    = TMath::Abs(cells->GetCellNumber(i));
 929     Double_t cellE = cells->GetAmplitude(i);
 930     fHCellE->Fill(cellE);
 931     if (cellE>cellMaxE)
 932       cellMaxE = cellE;
 933     cellMeanE += cellE;
 934
 935     Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1;
 936     Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta);
 937     if (!ret) {
 938       AliError(Form("Could not get cell index for %d", absID));
 939       continue;
 940     }
 941     ++cellModCount[iSM];
 942     Int_t iPhi=-1, iEta=-1;
 943     fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta);
 944     fHColuRow[iSM]->Fill(iEta,iPhi,1);
 945     fHColuRowE[iSM]->Fill(iEta,iPhi,cellE);
 946     fHCellFreqNoCut[iSM]->Fill(absID);
 947     if (cellE > 0.1) fHCellFreqCut100M[iSM]->Fill(absID);
 948     if (cellE > 0.3) fHCellFreqCut300M[iSM]->Fill(absID);
 949     if (fHCellCheckE && fHCellCheckE[absID])
 950       fHCellCheckE[absID]->Fill(cellE);
 951     fHCellFreqE[iSM]->Fill(absID, cellE);
 952   }
 953   fHCellH->Fill(cellMaxE);
 954   cellMeanE /= ncells;
 955   fHCellM->Fill(cellMeanE);
 956   fHCellM2->Fill(cellMeanE*ncells/24/48/nsm);
 957   for (Int_t i=0; i<nsm; ++i)
 958     fHCellMult[i]->Fill(cellModCount[i]);
 959 }
 960
 961 //________________________________________________________________________
 962 void AliAnalysisTaskEMCALPi0PbPb::FillClusHists()
 963 {
 964   // Fill histograms related to cluster properties.
 965
 966   TObjArray *clusters = fEsdClusters;
 967   if (!clusters)
 968     clusters = fAodClusters;
 969   if (!clusters)
 970     return;
 971
 972   Double_t vertex[3] = {0};
 973   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 974
 975   Int_t nclus = clusters->GetEntries();
 976   for(Int_t i = 0; i<nclus; ++i) {
 977     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
 978     if (!clus)
 979       continue;
 980     if (!clus->IsEMCAL())
 981       continue;
 982     TLorentzVector clusterVec;
 983     clus->GetMomentum(clusterVec,vertex);
 984     Double_t maxAxis = 0, minAxis = 0;
 985     GetSigma(clus,maxAxis,minAxis);
 986     clus->SetTOF(maxAxis);     // store sigma in TOF
 987     Double_t clusterEcc = 0;
 988     if (maxAxis > 0)
 989       clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
 990     clus->SetChi2(clusterEcc); // store ecc in chi2
 991     fHClustEccentricity->Fill(clusterEcc);
 992     fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi());
 993     fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt());
 994     fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E());
 995     fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis);
 996     fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E());
 997   }
 998 }
 999
1000 //________________________________________________________________________
1001 void AliAnalysisTaskEMCALPi0PbPb::FillNtuple()
1002 {
1003   // Fill ntuple.
1004
1005   if (!fNtuple)
1006     return;
1007
1008   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1009   if (fAodEv) {
1010     fHeader->fRun            = fAodEv->GetRunNumber();
1011     fHeader->fOrbit          = fAodEv->GetHeader()->GetOrbitNumber();
1012     fHeader->fPeriod         = fAodEv->GetHeader()->GetPeriodNumber();
1013     fHeader->fBx             = fAodEv->GetHeader()->GetBunchCrossNumber();
1014     fHeader->fL0             = fAodEv->GetHeader()->GetL0TriggerInputs();
1015     fHeader->fL1             = fAodEv->GetHeader()->GetL1TriggerInputs();
1016     fHeader->fL2             = fAodEv->GetHeader()->GetL2TriggerInputs();
1017     fHeader->fTrClassMask    = fAodEv->GetHeader()->GetTriggerMask();
1018     fHeader->fTrCluster      = fAodEv->GetHeader()->GetTriggerCluster();
1019     fHeader->fOffTriggers    = fAodEv->GetHeader()->GetOfflineTrigger();
1020     fHeader->fFiredTriggers  = fAodEv->GetHeader()->GetFiredTriggerClasses();
1021   } else {
1022     fHeader->fRun            = fEsdEv->GetRunNumber();
1023     fHeader->fOrbit          = fEsdEv->GetHeader()->GetOrbitNumber();
1024     fHeader->fPeriod         = fEsdEv->GetESDRun()->GetPeriodNumber();
1025     fHeader->fBx             = fEsdEv->GetHeader()->GetBunchCrossNumber();
1026     fHeader->fL0             = fEsdEv->GetHeader()->GetL0TriggerInputs();
1027     fHeader->fL1             = fEsdEv->GetHeader()->GetL1TriggerInputs();
1028     fHeader->fL2             = fEsdEv->GetHeader()->GetL2TriggerInputs();
1029     fHeader->fTrClassMask    = fEsdEv->GetHeader()->GetTriggerMask();
1030     fHeader->fTrCluster      = fEsdEv->GetHeader()->GetTriggerCluster();
1031     fHeader->fOffTriggers    = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
1032     fHeader->fFiredTriggers  = fEsdEv->GetFiredTriggerClasses();
1033   }
1034   AliCentrality *cent = InputEvent()->GetCentrality();
1035   fHeader->fV0Cent    = cent->GetCentralityPercentileUnchecked("V0M");
1036   fHeader->fCl1Cent   = cent->GetCentralityPercentileUnchecked("CL1");
1037   fHeader->fTrCent    = cent->GetCentralityPercentileUnchecked("TRK");
1038   fHeader->fCqual     = cent->GetQuality();
1039
1040   Double_t val = 0;
1041   TString trgclasses(fHeader->fFiredTriggers);
1042   for (Int_t j = 0; j<fTrClassNamesArr->GetEntries(); ++j) {
1043     const char *name = fTrClassNamesArr->At(j)->GetName();
1044     if (trgclasses.Contains(name))
1045       val += TMath::Power(2,j);
1046   }
1047   fHeader->fTcls = (UInt_t)val;
1048
1049   if (fAodEv) {
1050     am->LoadBranch("vertices");
1051     AliAODVertex *pv = fAodEv->GetPrimaryVertex();
1052     FillVertex(fPrimVert, pv);
1053     AliAODVertex *sv = fAodEv->GetPrimaryVertexSPD();
1054     FillVertex(fSpdVert, sv);
1055   } else {
1056     am->LoadBranch("PrimaryVertex.");
1057     const AliESDVertex *pv = fEsdEv->GetPrimaryVertexTracks();
1058     FillVertex(fPrimVert, pv);
1059     am->LoadBranch("SPDVertex.");
1060     const AliESDVertex *sv = fEsdEv->GetPrimaryVertexSPD();
1061     FillVertex(fSpdVert, sv);
1062     am->LoadBranch("TPCVertex.");
1063     const AliESDVertex *tv = fEsdEv->GetPrimaryVertexTPC();
1064     FillVertex(fTpcVert, tv);
1065   }
1066
1067   TObjArray *clusters = fEsdClusters;
1068   if (!clusters)
1069     clusters = fAodClusters;
1070   if (!clusters)
1071     return;
1072
1073   fClusters->Clear();
1074   Int_t nclus = clusters->GetEntries();
1075   for(Int_t i=0,ncl=0; i<nclus; ++i) {
1076     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1077     if (!clus)
1078       continue;
1079     if (!clus->IsEMCAL())
1080       continue;
1081     if (clus->E()<fMinE)
1082       continue;
1083
1084     AliStaCluster *cl = static_cast<AliStaCluster*>(fClusters->New(ncl++));
1085     cl->fE        = clus->E();
1086     Float_t pos[3];
1087     clus->GetPosition(pos);
1088     TVector3 vpos(pos);
1089     cl->fR        = vpos.Perp();
1090     cl->fEta      = vpos.Eta();
1091     cl->fPhi      = vpos.Phi();
1092     cl->fN        =  clus->GetNCells();
1093     cl->fN1       = GetNCells(clus,0.100);
1094     cl->fN3       = GetNCells(clus,0.300);
1095     Short_t id    = -1;
1096     Double_t emax = GetMaxCellEnergy(clus, id);
1097     cl->fIdMax    = id;
1098     cl->fEmax     = emax;
1099     cl->fDbc      = clus->GetDistanceToBadChannel();;
1100     cl->fDisp     = clus->GetDispersion();
1101     cl->fM20      = clus->GetM20();
1102     cl->fM02      = clus->GetM02();
1103     cl->fEcc      = clus->Chi2();   //eccentricity
1104     cl->fSig      = clus->GetTOF(); //sigma
1105     cl->fTrDz     = fClusProps[i].fTrDz;
1106     cl->fTrDr     = fClusProps[i].fTrDr;;
1107     cl->fTrEp     = fClusProps[i].fTrEp;;
1108     cl->fTrIso    = fClusProps[i].fTrIso;
1109     cl->fCeIso    = fClusProps[i].fCellIso;
1110   }
1111   fNtuple->Fill();
1112 }
1113
1114 //________________________________________________________________________
1115 void AliAnalysisTaskEMCALPi0PbPb::FillPionHists()
1116 {
1117   // Fill histograms related to pions.
1118
1119   Double_t vertex[3] = {0};
1120   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1121
1122   TObjArray *clusters = fEsdClusters;
1123   if (!clusters)
1124     clusters = fAodClusters;
1125
1126   if (clusters) {
1127     TLorentzVector clusterVec1;
1128     TLorentzVector clusterVec2;
1129     TLorentzVector pionVec;
1130
1131     Int_t nclus = clusters->GetEntries();
1132     for (Int_t i = 0; i<nclus; ++i) {
1133       AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i));
1134       if (!clus1)
1135         continue;
1136       if (!clus1->IsEMCAL())
1137         continue;
1138       if (clus1->E()<fMinE)
1139         continue;
1140       if (clus1->GetNCells()<fNminCells)
1141         continue;
1142       if (GetMaxCellEnergy(clus1)/clus1->E()<fMinErat)
1143         continue;
1144       if (clus1->Chi2()<fMinEcc) // eccentricity cut
1145         continue;
1146       clus1->GetMomentum(clusterVec1,vertex);
1147       for (Int_t j = i+1; j<nclus; ++j) {
1148         AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j));
1149         if (!clus2)
1150           continue;
1151         if (!clus2->IsEMCAL())
1152           continue;
1153         if (clus2->E()<fMinE)
1154           continue;
1155         if (clus2->GetNCells()<fNminCells)
1156           continue;
1157         if (GetMaxCellEnergy(clus2)/clus2->E()<fMinErat)
1158           continue;
1159         if (clus2->Chi2()<fMinEcc) // eccentricity cut
1160           continue;
1161         clus2->GetMomentum(clusterVec2,vertex);
1162         pionVec = clusterVec1 + clusterVec2;
1163         Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E();
1164         Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect());
1165         if (pionZgg < fAsymMax) {
1166           fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi());
1167           fHPionMggPt->Fill(pionVec.M(),pionVec.Pt());
1168           fHPionMggAsym->Fill(pionVec.M(),pionZgg);
1169           fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle);
1170           Int_t bin = fPtRanges->FindBin(pionVec.Pt());
1171           fHPionInvMasses[bin]->Fill(pionVec.M());
1172         }
1173       }
1174     }
1175   }
1176 }
1177
1178 //________________________________________________________________________
1179 void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists()
1180 {
1181   // Fill other histograms.
1182 }
1183
1184 //__________________________________________________________________________________________________
1185 void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliESDVertex *esdv)
1186 {
1187   // Fill vertex from ESD vertex info.
1188
1189   v->fVx   = esdv->GetX();
1190   v->fVy   = esdv->GetY();
1191   v->fVz   = esdv->GetZ();
1192   v->fVc   = esdv->GetNContributors();
1193   v->fDisp = esdv->GetDispersion();
1194   v->fZres = esdv->GetZRes();
1195   v->fChi2 = esdv->GetChi2();
1196   v->fSt   = esdv->GetStatus();
1197   v->fIs3D = esdv->IsFromVertexer3D();
1198   v->fIsZ  = esdv->IsFromVertexerZ();
1199 }
1200
1201 //__________________________________________________________________________________________________
1202 void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliAODVertex *aodv)
1203 {
1204   // Fill vertex from AOD vertex info.
1205
1206   v->fVx   = aodv->GetX();
1207   v->fVy   = aodv->GetY();
1208   v->fVz   = aodv->GetZ();
1209   v->fVc   = aodv->GetNContributors();
1210   v->fChi2 = aodv->GetChi2();
1211 }
1212
1213 //________________________________________________________________________
1214 Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
1215 {
1216   // Compute isolation based on cell content.
1217
1218   AliVCaloCells *cells = fEsdCells;
1219   if (!cells)
1220     cells = fAodCells;
1221   if (!cells)
1222     return 0;
1223
1224   Double_t cellIsolation = 0;
1225   Double_t rad2 = radius*radius;
1226   Int_t ncells = cells->GetNumberOfCells();
1227   for (Int_t i = 0; i<ncells; ++i) {
1228     Int_t absID    = TMath::Abs(cells->GetCellNumber(i));
1229     Double_t cellE = cells->GetAmplitude(i);
1230     Float_t eta, phi;
1231     fGeom->EtaPhiFromIndex(absID,eta,phi);
1232     Double_t phidiff = TVector2::Phi_0_2pi(phi-cPhi);
1233     Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1234     if(dist>rad2)
1235       continue;
1236     cellIsolation += cellE;
1237   }
1238   return cellIsolation;
1239 }
1240
1241 //________________________________________________________________________
1242 Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(AliVCluster *cluster, Short_t &id) const
1243 {
1244   // Get maximum energy of attached cell.
1245
1246   id = -1;
1247   Double_t maxe = 0;
1248   Int_t ncells = cluster->GetNCells();
1249   if (fEsdCells) {
1250     for (Int_t i=0; i<ncells; i++) {
1251       Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1252       if (e>maxe) {
1253         maxe = e;
1254         id   = cluster->GetCellAbsId(i);
1255       }
1256     }
1257   } else {
1258     for (Int_t i=0; i<ncells; i++) {
1259       Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1260       if (e>maxe)
1261         maxe = e;
1262         id   = cluster->GetCellAbsId(i);
1263     }
1264   }
1265   return maxe;
1266 }
1267
1268 //________________________________________________________________________
1269 void AliAnalysisTaskEMCALPi0PbPb::GetSigma(AliVCluster *c, Double_t& sigmaMax, Double_t &sigmaMin) const
1270 {
1271   // Calculate the (E) weighted variance along the longer (eigen) axis.
1272
1273   sigmaMax = 0;          // cluster variance along its longer axis
1274   sigmaMin = 0;          // cluster variance along its shorter axis
1275   Double_t Ec  = c->E(); // cluster energy
1276   if(Ec<=0)
1277     return;
1278   Double_t Xc  = 0 ;     // cluster first moment along X
1279   Double_t Yc  = 0 ;     // cluster first moment along Y
1280   Double_t Sxx = 0 ;     // cluster second central moment along X (variance_X^2)
1281   Double_t Sxy = 0 ;     // cluster second central moment along Y (variance_Y^2)
1282   Double_t Syy = 0 ;     // cluster covariance^2
1283
1284   AliVCaloCells *cells = fEsdCells;
1285   if (!cells)
1286     cells = fAodCells;
1287
1288   if (!cells)
1289     return;
1290
1291   Int_t ncells = c->GetNCells();
1292   if (ncells==1)
1293     return;
1294
1295   TVector3 pos;
1296   for(Int_t j=0; j<ncells; ++j) {
1297     Int_t id = TMath::Abs(c->GetCellAbsId(j));
1298     fGeom->GetGlobal(id,pos);
1299     Double_t cellen = cells->GetCellAmplitude(id);
1300     Xc  += cellen*pos.X();
1301     Yc  += cellen*pos.Y();
1302     Sxx += cellen*pos.X()*pos.X();
1303     Syy += cellen*pos.Y()*pos.Y();
1304     Sxy += cellen*pos.X()*pos.Y();
1305   }
1306   Xc  /= Ec;
1307   Yc  /= Ec;
1308   Sxx /= Ec;
1309   Syy /= Ec;
1310   Sxy /= Ec;
1311   Sxx -= Xc*Xc;
1312   Syy -= Yc*Yc;
1313   Sxy -= Xc*Yc;
1314   Sxx = TMath::Abs(Sxx);
1315   Syy = TMath::Abs(Syy);
1316   sigmaMax = (Sxx + Syy + TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
1317   sigmaMax = TMath::Sqrt(TMath::Abs(sigmaMax));
1318   sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
1319   sigmaMin = TMath::Sqrt(TMath::Abs(sigmaMin));
1320 }
1321
1322 //________________________________________________________________________
1323 Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(AliVCluster *c, Double_t emin) const
1324 {
1325   // Calculate number of attached cells above emin.
1326
1327   AliVCaloCells *cells = fEsdCells;
1328   if (!cells)
1329     cells = fAodCells;
1330   if (!cells)
1331     return 0;
1332
1333   Int_t n = 0;
1334   Int_t ncells = c->GetNCells();
1335   for(Int_t j=0; j<ncells; ++j) {
1336     Int_t id = TMath::Abs(c->GetCellAbsId(j));
1337     Double_t cellen = cells->GetCellAmplitude(id);
1338     if (cellen>=emin)
1339       ++n;
1340   }
1341   return n;
1342 }
1343
1344 //________________________________________________________________________
1345 Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
1346 {
1347   // Compute isolation based on tracks.
1348
1349   Double_t trkIsolation = 0;
1350   Double_t rad2 = radius*radius;
1351   Int_t ntrks = fSelTracks->GetEntries();
1352   for(Int_t j = 0; j<ntrks; ++j) {
1353     AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
1354     if (!track)
1355       continue;
1356     Float_t eta = track->Eta();
1357     Float_t phi = track->Phi();
1358     Double_t phidiff = TVector2::Phi_0_2pi(phi-cPhi);
1359     Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1360     if(dist>rad2)
1361       continue;
1362     trkIsolation += track->Pt();
1363   }
1364   return trkIsolation;
1365 }