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