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 <TH1F.h>
   9 #include <TH2F.h>
  10 #include <TList.h>
  11 #include <TLorentzVector.h>
  12 #include <TNtuple.h>
  13 #include "AliAODEvent.h"
  14 #include "AliAODVertex.h"
  15 #include "AliAnalysisManager.h"
  16 #include "AliAnalysisTaskEMCALClusterizeFast.h"
  17 #include "AliCentrality.h"
  18 #include "AliEMCALGeoUtils.h"
  19 #include "AliESDEvent.h"
  20 #include "AliESDVertex.h"
  21 #include "AliLog.h"
  22
  23 ClassImp(AliAnalysisTaskEMCALPi0PbPb)
  24
  25 //________________________________________________________________________
  26 AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name)
  27   : AliAnalysisTaskSE(name),
  28     fAsymMax(1),
  29     fCentVar("V0M"),
  30     fCentFrom(0),
  31     fCentTo(100),
  32     fVtxZMin(-10),
  33     fVtxZMax(+10),
  34     fUseQualFlag(1),
  35     fClusName(),
  36     fDoNtuple(0),
  37     fDoAfterburner(0),
  38     fNminCells(1),
  39     fNEvs(0),
  40     fGeom(0),
  41     fOutput(0),
  42     fEsdEv(0),
  43     fAodEv(0),
  44     fRecPoints(0),
  45     fEsdClusters(0),
  46     fEsdCells(0),
  47     fAodClusters(0),
  48     fAodCells(0),
  49     fPtRanges(0),
  50     fNtuple(0),
  51     fHCuts(0x0),
  52     fHVertexZ(0x0),
  53     fHVertexZ2(0x0),
  54     fHCent(0x0),
  55     fHCentQual(0x0),
  56     fHColuRow(0x0),
  57     fHColuRowE(0x0),
  58     fHCellMult(0x0),
  59     fHCellE(0x0),
  60     fHCellH(0x0),
  61     fHCellM(0x0),
  62     fHCellM2(0x0),
  63     fHClustEccentricity(0),
  64     fHClustEtaPhi(0x0),
  65     fHClustEnergyPt(0x0),
  66     fHClustEnergySigma(0x0),
  67     fHClustSigmaSigma(0x0),
  68     fHClustNCellEnergyRatio(0x0),
  69     fHPionEtaPhi(0x0),
  70     fHPionMggPt(0x0),
  71     fHPionMggAsym(0x0),
  72     fHPionMggDgg(0x0)
  73 {
  74   // Constructor.
  75
  76   if (!name)
  77     return;
  78   SetName(name);
  79   DefineInput(0, TChain::Class());
  80   DefineOutput(1, TList::Class());
  81   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,EMCALCells. "
  82                "AOD:header,vertices,emcalCells";
  83 }
  84
  85 //________________________________________________________________________
  86 AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb()
  87 {
  88   // Destructor.
  89
  90   delete fOutput; fOutput = 0;
  91   delete fPtRanges; fPtRanges = 0;
  92   delete fGeom; fGeom = 0;
  93   delete [] fHColuRow;
  94   delete [] fHColuRowE;
  95   delete [] fHCellMult;
  96 }
  97
  98 //________________________________________________________________________
  99 void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects()
 100 {
 101   // Create user objects here.
 102
 103   fOutput = new TList();
 104   fOutput->SetOwner();
 105
 106   if (fDoNtuple) {
 107     TFile *f = OpenFile(1);
 108     if (f)
 109       fNtuple = new TNtuple(Form("nt%.0fto%.0f",fCentFrom,fCentTo),"nt",
 110                             "run:evt:cent:pt:e:emax:n:db:disp:mn:ms:chi:cpv:ecc:sig:eta:phi");
 111   }
 112
 113   // histograms
 114   TH1::SetDefaultSumw2(kTRUE);
 115   TH2::SetDefaultSumw2(kTRUE);
 116   fHCuts = new TH1F("hEventCuts","",4,0.5,4.5);
 117   fHCuts->GetXaxis()->SetBinLabel(1,"All (PS)");
 118   fHCuts->GetXaxis()->SetBinLabel(2,Form("%s: %.0f-%.0f",fCentVar.Data(),fCentFrom,fCentTo));
 119   fHCuts->GetXaxis()->SetBinLabel(3,"QFlag");
 120   fHCuts->GetXaxis()->SetBinLabel(4,Form("zvtx: %.0f-%.0f",fVtxZMin,fVtxZMax));
 121   fOutput->Add(fHCuts);
 122   fHVertexZ = new TH1F("hVertexZBeforeCut","",100,-25,25);
 123   fHVertexZ->SetXTitle("z [cm]");
 124   fOutput->Add(fHVertexZ);
 125   fHVertexZ2 = new TH1F("hVertexZAfterCut","",100,-25,25);
 126   fHVertexZ2->SetXTitle("z [cm]");
 127   fOutput->Add(fHVertexZ2);
 128   fHCent = new TH1F("hCentBeforeCut","",101,-1,100);
 129   fHCent->SetXTitle(fCentVar.Data());
 130   fOutput->Add(fHCent);
 131   fHCentQual = new TH1F("hCentAfterCut","",101,-1,100);
 132   fHCentQual->SetXTitle(fCentVar.Data());
 133   fOutput->Add(fHCentQual);
 134
 135   // histograms for cells
 136   fHColuRow   = new TH2F*[4];
 137   fHColuRowE  = new TH2F*[4];
 138   fHCellMult  = new TH1F*[4];
 139   for (Int_t i = 0; i<4; ++i) {
 140     fHColuRow[i] = new TH2F(Form("hColRow_Mod%d",i),"",49,0,49,25,0,25);
 141     fHColuRow[i]->SetTitle(Form("Module %d: Occupancy", i));
 142     fHColuRow[i]->SetXTitle("col (i#eta)");
 143     fHColuRow[i]->SetYTitle("row (i#phi)");
 144     fHColuRowE[i] = new TH2F(Form("hColRowE_Mod%d", i),"",49,0,49,25,0,25);
 145     fHColuRowE[i]->SetTitle(Form("Module %d: Cell energy",i));
 146     fHColuRowE[i]->SetXTitle("col (i#eta)");
 147     fHColuRowE[i]->SetYTitle("row (i#phi)");
 148     fHCellMult[i] = new TH1F(Form("hCellMult_Mod%d",i),"",1000,0,1000);
 149     fHCellMult[i]->SetTitle(Form("Module %d: Cell multiplicity",i));
 150     fHCellMult[i]->SetXTitle("# of cells");
 151     fOutput->Add(fHColuRow[i]);
 152     fOutput->Add(fHColuRowE[i]);
 153     fOutput->Add(fHCellMult[i]);
 154   }
 155   fHCellE = new TH1F("hCellE","",150,0.,15.);
 156   fHCellE->SetXTitle("E_{cell} [GeV]");
 157   fOutput->Add(fHCellE);
 158   fHCellH = new TH1F ("fHCellHighestE","",150,0.,15.);
 159   fHCellH->SetXTitle("E^{max}_{cell} [GeV]");
 160   fOutput->Add(fHCellH);
 161   fHCellM = new TH1F ("fHCellMeanEperHitCell","",250,0.,2.5);
 162   fHCellM->SetXTitle("#LT E_{cell}#GT [GeV]");
 163   fOutput->Add(fHCellM);
 164   fHCellM2 = new TH1F ("fHCellMeanEperAllCells","",250,0.,1);
 165   fHCellM2->SetXTitle("1/N_{cells} #Sigma E_{cell} [GeV]");
 166   fOutput->Add(fHCellM2);
 167
 168   // histograms for clusters
 169   fHClustEccentricity = new TH1F("hClustEccentricity","",100,-0.1,1.1);
 170   fHClustEccentricity->SetXTitle("#epsilon_{C}");
 171   fOutput->Add(fHClustEccentricity);
 172   fHClustEtaPhi = new TH2F("hClustEtaPhi","",500,-0.8,0.8,500,1.2,2.2);
 173   fHClustEtaPhi->SetXTitle("#eta");
 174   fHClustEtaPhi->SetYTitle("#varphi");
 175   fOutput->Add(fHClustEtaPhi);
 176   fHClustEnergyPt = new TH2F("hClustEnergyPt","",250,0,50,250,0,50);
 177   fHClustEnergyPt->SetXTitle("E [GeV]");
 178   fHClustEnergyPt->SetYTitle("p_{T} [GeV/c]");
 179   fOutput->Add(fHClustEnergyPt);
 180   fHClustEnergySigma = new TH2F("hClustEnergySigma","",100,0,100,500,0,50);
 181   fHClustEnergySigma->SetXTitle("E_{C} * #sigma_{max} [GeV*cm]");
 182   fHClustEnergySigma->SetYTitle("E_{C} [GeV]");
 183   fOutput->Add(fHClustEnergySigma);
 184   fHClustSigmaSigma = new TH2F("hClustSigmaSigma","",500,0,50,500,0,50);
 185   fHClustSigmaSigma->SetXTitle("#lambda_{0} [cm]");
 186   fHClustSigmaSigma->SetYTitle("#sigma_{max} [cm]");
 187   fOutput->Add(fHClustSigmaSigma);
 188   fHClustNCellEnergyRatio = new TH2F("hClustNCellEnergyRatio","",27,-0.5,26.5,101,-0.05,1.05);
 189   fHClustNCellEnergyRatio->SetXTitle("N_{cells}");
 190   fHClustNCellEnergyRatio->SetYTitle("E^{max}_{cell}/E_{clus}");
 191   fOutput->Add(fHClustNCellEnergyRatio);
 192
 193   // histograms for pion candidates
 194   fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100,1.2,2.2);
 195   fHPionEtaPhi->SetXTitle("#eta_{#gamma#gamma}");
 196   fHPionEtaPhi->SetYTitle("#varphi_{#gamma#gamma}");
 197   fOutput->Add(fHPionEtaPhi);
 198   fHPionMggPt = new TH2F("hPionMggPt","",1000,0,2,100,0,20.0);
 199   fHPionMggPt->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 200   fHPionMggPt->SetYTitle("p_{T}^{#gamma#gamma} [GeV/c]");
 201   fOutput->Add(fHPionMggPt);
 202   fHPionMggAsym = new TH2F("hPionMggAsym","",1000,0,2,100,0,1);
 203   fHPionMggAsym->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 204   fHPionMggAsym->SetYTitle("Z_{#gamma#gamma} [GeV]");
 205   fOutput->Add(fHPionMggAsym);
 206   fHPionMggDgg = new TH2F("hPionMggDgg","",1000,0,2,100,0,10);
 207   fHPionMggDgg->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 208   fHPionMggDgg->SetYTitle("opening angle [grad]");
 209   fOutput->Add(fHPionMggDgg);
 210   const Int_t nbins = 20;
 211   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};
 212   fPtRanges = new TAxis(nbins-1,xbins);
 213   for (Int_t i = 0; i<=nbins; ++i) {
 214     fHPionInvMasses[i] = new TH1F(Form("hPionInvMass%d",i),"",1000,0,2);
 215     fHPionInvMasses[i]->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 216     if (i==0)
 217       fHPionInvMasses[i]->SetTitle(Form("0 < p_{T}^{#gamma#gamma} <%.1f",xbins[0]));
 218     else if (i==nbins)
 219       fHPionInvMasses[i]->SetTitle(Form("p_{T}^{#gamma#gamma} > 50"));
 220     else
 221       fHPionInvMasses[i]->SetTitle(Form("%.1f < p_{T}^{#gamma#gamma} <%.1f",xbins[i-1],xbins[i]));
 222     fOutput->Add(fHPionInvMasses[i]);
 223   }
 224
 225   PostData(1, fOutput);
 226 }
 227
 228 //________________________________________________________________________
 229 void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *)
 230 {
 231   // Called for each event.
 232
 233   if (!InputEvent())
 234     return;
 235
 236   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 237   fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent());
 238   if (fEsdEv) {
 239     am->LoadBranch("AliESDRun.");
 240     am->LoadBranch("AliESDHeader.");
 241   } else {
 242     fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
 243     am->LoadBranch("header");
 244   }
 245
 246   if (!fGeom) { // set misalignment matrices (stored in first event)
 247     fGeom = new AliEMCALGeoUtils("EMCAL_FIRSTYEARV1","EMCAL");
 248     if (fEsdEv) {
 249       for (Int_t i=0; i<4; ++i)
 250         fGeom->SetMisalMatrix(fEsdEv->GetESDRun()->GetEMCALMatrix(i),i);
 251     } else {
 252       for (Int_t i=0; i<4; ++i)
 253         fGeom->SetMisalMatrix(fAodEv->GetHeader()->GetEMCALMatrix(i),i);
 254     }
 255     fGeom->GetEMCGeometry();
 256   }
 257
 258   Int_t cut = 1;
 259   fHCuts->Fill(cut++);
 260
 261   const AliCentrality *centP = InputEvent()->GetCentrality();
 262   Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar);
 263   fHCent->Fill(cent);
 264   if (cent<fCentFrom||cent>fCentTo)
 265     return;
 266
 267   fHCuts->Fill(cut++);
 268
 269   if (fUseQualFlag) {
 270     if (centP->GetQuality()>0)
 271       return;
 272   }
 273
 274   fHCentQual->Fill(cent);
 275   fHCuts->Fill(cut++);
 276
 277   // count number of accepted events
 278   ++fNEvs;
 279
 280   if (fEsdEv) {
 281     am->LoadBranch("PrimaryVertex.");
 282     am->LoadBranch("SPDVertex.");
 283     am->LoadBranch("TPCVertex.");
 284   } else {
 285     fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
 286     am->LoadBranch("vertices");
 287   }
 288
 289   const AliVVertex *vertex = InputEvent()->GetPrimaryVertex();
 290   if (!vertex)
 291     return;
 292
 293   fHVertexZ->Fill(vertex->GetZ());
 294
 295   if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax)
 296     return;
 297
 298   fHCuts->Fill(cut++);
 299   fHVertexZ2->Fill(vertex->GetZ());
 300
 301   fRecPoints   = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
 302   fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set of if clusters are attached
 303   fEsdCells    = 0; // will be set if ESD input used
 304   fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set of if clusters are attached
 305                     //             or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode
 306   fAodCells    = 0; // will be set if AOD input used
 307
 308   // deal with special output from AliAnalysisTaskEMCALClusterizeFast first
 309   Bool_t clusattached = 0;
 310   Bool_t recalibrated = 0;
 311   if (1 && !fClusName.IsNull()) {
 312     AliAnalysisTaskEMCALClusterizeFast *cltask = 0;
 313     TObjArray *ts = am->GetTasks();
 314     cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName));
 315     if (cltask && cltask->GetClusters()) {
 316       fRecPoints = const_cast<TObjArray*>(cltask->GetClusters());
 317       clusattached = cltask->GetAttachClusters();
 318       if (cltask->GetCalibData()!=0)
 319         recalibrated = kTRUE;
 320     }
 321   }
 322   if (1 && AODEvent() && !fClusName.IsNull()) {
 323     TList *l = AODEvent()->GetList();
 324     TClonesArray *clus = 0;
 325     if (l) {
 326       clus = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
 327       fAodClusters = clus;
 328     }
 329   }
 330
 331   if (fEsdEv) { // ESD input mode
 332     if (1 && (!fRecPoints||clusattached)) {
 333       if (!clusattached)
 334         am->LoadBranch("CaloClusters");
 335       TList *l = fEsdEv->GetList();
 336       TClonesArray *clus = 0;
 337       if (l) {
 338         clus = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
 339         fEsdClusters = clus;
 340       }
 341     }
 342     if (1) {
 343       if (!recalibrated)
 344         am->LoadBranch("EMCALCells.");
 345       fEsdCells = fEsdEv->GetEMCALCells();
 346     }
 347   } else if (fAodEv) { // AOD input mode
 348     if (1 && (!fAodClusters || clusattached)) {
 349       if (!clusattached)
 350         am->LoadBranch("caloClusters");
 351       TList *l = fAodEv->GetList();
 352       TClonesArray *clus = 0;
 353       if (l) {
 354         clus = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
 355         fAodClusters = clus;
 356       }
 357     }
 358     if (1) {
 359       if (!recalibrated)
 360         am->LoadBranch("emcalCells");
 361       fAodCells = fAodEv->GetEMCALCells();
 362     }
 363   } else {
 364     AliFatal("Impossible to not have either pointer to ESD or AOD event");
 365   }
 366
 367   if (1) {
 368     AliDebug(2,Form("fRecPoints   set: %p", fRecPoints));
 369     AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters));
 370     AliDebug(2,Form("fEsdCells    set: %p", fEsdCells));
 371     AliDebug(2,Form("fAodClusters set: %p", fAodClusters));
 372     AliDebug(2,Form("fAodCells    set: %p", fAodCells));
 373   }
 374
 375   if (fDoAfterburner)
 376     ClusterAfterburner();
 377
 378   FillCellHists();
 379   FillClusHists();
 380   FillPionHists();
 381
 382   PostData(1, fOutput);
 383 }
 384
 385 //________________________________________________________________________
 386 void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *)
 387 {
 388   // Terminate called at the end of analysis.
 389
 390   if (fNtuple) {
 391     TFile *f = OpenFile(1);
 392     if (f)
 393       fNtuple->Write();
 394   }
 395
 396   AliInfo(Form("\nAccepted %lld events", fNEvs));
 397 }
 398
 399 //________________________________________________________________________
 400 void AliAnalysisTaskEMCALPi0PbPb::FillCellHists()
 401 {
 402   // Fill histograms related to cell properties.
 403
 404   AliVCaloCells *cells = fEsdCells;
 405   if (!cells)
 406     cells = fAodCells;
 407
 408   if (cells) {
 409     Int_t cellModCount[4] = {0,0,0,0};
 410     Double_t cellMaxE = 0;
 411     Double_t cellMeanE = 0;
 412     Int_t ncells = cells->GetNumberOfCells();
 413     if (ncells==0)
 414       return;
 415
 416     for (Int_t i = 0; i<ncells; ++i ) {
 417       Int_t absID    = TMath::Abs(cells->GetCellNumber(i));
 418       Double_t cellE = cells->GetAmplitude(i);
 419       fHCellE->Fill(cellE);
 420       if (cellE>cellMaxE)
 421         cellMaxE = cellE;
 422       cellMeanE += cellE;
 423
 424       Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1;
 425       Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta);
 426       if (!ret) {
 427         AliError(Form("Could not get cell index for %d", absID));
 428         continue;
 429       }
 430       ++cellModCount[iSM];
 431       Int_t iPhi=-1, iEta=-1;
 432       fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta);
 433       fHColuRow[iSM]->Fill(iEta,iPhi,1);
 434       fHColuRowE[iSM]->Fill(iEta,iPhi,cellE);
 435     }
 436     fHCellH->Fill(cellMaxE);
 437     cellMeanE /= ncells;
 438     fHCellM->Fill(cellMeanE);
 439     fHCellM2->Fill(cellMeanE*ncells/24/48/4); //hard-coded but there is a way to figure out from geometry
 440     for (Int_t i=0; i<4; ++i)
 441       fHCellMult[i]->Fill(cellModCount[i]);
 442   }
 443 }
 444
 445 //________________________________________________________________________
 446 void AliAnalysisTaskEMCALPi0PbPb::FillClusHists()
 447 {
 448   // Fill histograms related to cluster properties.
 449   Double_t clusterEcc = 0;
 450
 451   Double_t vertex[3] = {0,0,0};
 452   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 453
 454   TObjArray *clusters = fEsdClusters;
 455   if (!clusters)
 456     clusters = fAodClusters;
 457
 458   if (clusters) {
 459     TLorentzVector clusterVec;
 460     Int_t nclus = clusters->GetEntries();
 461     for(Int_t i = 0; i<nclus; ++i) {
 462       AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
 463       if (!clus)
 464         continue;
 465       if (!clus->IsEMCAL())
 466         continue;
 467       clus->GetMomentum(clusterVec,vertex);
 468       Double_t maxAxis = 0, minAxis = 0;
 469       GetSigma(clus,maxAxis,minAxis);
 470       if (maxAxis > 0)
 471         clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
 472
 473       fHClustEccentricity->Fill(clusterEcc);
 474       fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi());
 475       fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt());
 476       fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E());
 477       fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis);
 478       fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E());
 479       if (fNtuple) {
 480         Float_t vals[17];
 481         vals[0]  = InputEvent()->GetRunNumber();
 482         vals[1]  = (((UInt_t)InputEvent()->GetOrbitNumber()  << 12) | (UInt_t)InputEvent()->GetBunchCrossNumber());
 483         if (vals[1]<1)
 484           vals[1] = fNEvs;
 485         vals[2]  = InputEvent()->GetCentrality()->GetCentralityPercentileUnchecked(fCentVar);
 486         vals[3]  = clusterVec.Pt();
 487         vals[4]  = clusterVec.E();
 488         vals[5]  = GetMaxCellEnergy(clus);
 489         vals[6]  = clus->GetNCells();
 490         vals[7]  = clus->GetDistanceToBadChannel();
 491         vals[8]  = clus->GetDispersion();
 492         vals[9]  = clus->GetM20();
 493         vals[10] = clus->GetM02();
 494         vals[11] = clus->Chi2();
 495         vals[12] = clus->GetEmcCpvDistance();
 496         vals[13] = clusterEcc;
 497         vals[14] = maxAxis;
 498         vals[15] = clusterVec.Eta();
 499         vals[16] = clusterVec.Phi();
 500         fNtuple->Fill(vals);
 501       }
 502     }
 503   }
 504 }
 505
 506 //________________________________________________________________________
 507 void AliAnalysisTaskEMCALPi0PbPb::FillPionHists()
 508 {
 509   // Fill histograms related to pions.
 510
 511   Double_t vertex[3] = {0,0,0};
 512   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 513
 514   TObjArray *clusters = fEsdClusters;
 515   if (!clusters)
 516     clusters = fAodClusters;
 517
 518   if (clusters) {
 519     TLorentzVector clusterVec1;
 520     TLorentzVector clusterVec2;
 521     TLorentzVector pionVec;
 522
 523     Int_t nclus = clusters->GetEntries();
 524     for (Int_t i = 0; i<nclus; ++i) {
 525       AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i));
 526       if (!clus1)
 527         continue;
 528       if (!clus1->IsEMCAL())
 529         continue;
 530       if (clus1->E()<0.010)
 531         continue;
 532       if (clus1->GetNCells()<fNminCells)
 533         continue;
 534       clus1->GetMomentum(clusterVec1,vertex);
 535       for (Int_t j = i+1; j<nclus; ++j) {
 536         AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j));
 537         if (!clus2)
 538           continue;
 539         if (!clus2->IsEMCAL())
 540           continue;
 541         if (clus2->E()<0.010)
 542           continue;
 543         if (clus2->GetNCells()<fNminCells)
 544           continue;
 545         clus2->GetMomentum(clusterVec2,vertex);
 546         pionVec = clusterVec1 + clusterVec2;
 547         Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E();
 548         Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect());
 549         if (pionZgg < fAsymMax) {
 550           fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi());
 551           fHPionMggPt->Fill(pionVec.M(),pionVec.Pt());
 552           fHPionMggAsym->Fill(pionVec.M(),pionZgg);
 553           fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle);
 554           Int_t bin = fPtRanges->FindBin(pionVec.Pt());
 555           fHPionInvMasses[bin]->Fill(pionVec.M());
 556         }
 557       }
 558     }
 559   }
 560 }
 561
 562 //________________________________________________________________________
 563 void  AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner()
 564 {
 565   //
 566
 567   AliVCaloCells *cells = fEsdCells;
 568   if (!cells)
 569     cells = fAodCells;
 570
 571   if (!cells)
 572     return;
 573
 574   Int_t ncells = cells->GetNumberOfCells();
 575   if (ncells<=0)
 576     return;
 577
 578   Double_t cellMeanE = 0, cellSigE = 0;
 579   for (Int_t i = 0; i<ncells; ++i ) {
 580     Double_t cellE = cells->GetAmplitude(i);
 581     cellMeanE += cellE;
 582     cellSigE += cellE*cellE;
 583   }
 584   cellMeanE /= ncells;
 585   cellSigE /= ncells;
 586   cellSigE -= cellMeanE*cellMeanE;
 587   if (cellSigE<0)
 588     cellSigE = 0;
 589   cellSigE = TMath::Sqrt(cellSigE / ncells);
 590
 591   Double_t subE = cellMeanE - 7*cellSigE;
 592   if (subE<0)
 593     return;
 594
 595   for (Int_t i = 0; i<ncells; ++i) {
 596     Short_t id=-1;
 597     Double_t amp=0,time=0;
 598     if (!cells->GetCell(i, id, amp, time))
 599       continue;
 600     amp -= cellMeanE;
 601     if (amp<0.001)
 602       amp = 0;
 603     cells->SetCell(i, id, amp, time);
 604   }
 605
 606   TObjArray *clusters = fEsdClusters;
 607   if (!clusters)
 608     clusters = fAodClusters;
 609
 610   if (!clusters)
 611     return;
 612
 613   Int_t nclus = clusters->GetEntries();
 614   for (Int_t i = 0; i<nclus; ++i) {
 615     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
 616     Int_t nc = clus->GetNCells();
 617     Double_t clusE = 0;
 618     Short_t nc2 = 0;
 619     UShort_t ids[100] = {0};
 620     Double_t fra[100] = {0};
 621     for (Int_t j = 0; j<nc; ++j) {
 622       Short_t id = TMath::Abs(clus->GetCellAbsId(j));
 623       Double_t cen = cells->GetCellAmplitude(id);
 624       clusE += cen;
 625       if (cen>0) {
 626         ids[nc] = id;
 627         ++nc;
 628       }
 629     }
 630     if (clusE<=0) {
 631       clusters->RemoveAt(i);
 632       continue;
 633     }
 634     for (Int_t j = 0; j<nc2; ++j) {
 635       Short_t id = ids[j];
 636       Double_t cen = cells->GetCellAmplitude(id);
 637       fra[j] = cen/clusE;
 638     }
 639     clus->SetE(clusE);
 640     AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus);
 641     if (aodclus) {
 642       aodclus->Clear("");
 643       aodclus->SetNCells(nc2);
 644       aodclus->SetCellsAmplitudeFraction(fra);
 645       aodclus->SetCellsAbsId(ids);
 646     }
 647   }
 648   clusters->Compress();
 649 }
 650
 651 //________________________________________________________________________
 652 Double_t  AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(AliVCluster *cluster)
 653 {
 654   // Get maximum energy of attached cell.
 655
 656   Double_t maxe = 0;
 657   Int_t ncells = cluster->GetNCells();
 658   if (fEsdCells) {
 659     for (Int_t i=0; i<ncells; i++) {
 660       Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
 661       if (e>maxe)
 662         maxe = e;
 663     }
 664   } else {
 665     for (Int_t i=0; i<ncells; i++) {
 666       Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
 667       if (e>maxe)
 668         maxe = e;
 669     }
 670   }
 671   return maxe;
 672 }
 673
 674 //________________________________________________________________________
 675 void AliAnalysisTaskEMCALPi0PbPb::GetSigma(AliVCluster *cluster, Double_t& sigmaMax, Double_t &sigmaMin)
 676 {
 677   // Calculate the (E) weighted variance along the longer (eigen) axis.
 678
 679   sigmaMax = 0;               // cluster variance along its longer axis
 680   sigmaMin = 0;               // cluster variance along its shorter axis
 681   Double_t Ec = cluster->E(); // cluster energy
 682   Double_t Xc = 0 ;           // cluster first moment along X
 683   Double_t Yc = 0 ;           // cluster first moment along Y
 684   Double_t Sxx = 0 ;          // cluster second central moment along X (variance_X^2)
 685   Double_t Sxy = 0 ;          // cluster second central moment along Y (variance_Y^2)
 686   Double_t Syy = 0 ;          // cluster covariance^2
 687
 688   AliVCaloCells *cells = fEsdCells;
 689   if (!cells)
 690     cells = fAodCells;
 691
 692   if (!cells)
 693     return;
 694
 695   Int_t ncells = cluster->GetNCells();
 696   if (ncells==1)
 697     return;
 698
 699   TVector3 pos;
 700   for(Int_t j=0; j<ncells; ++j) {
 701     Int_t id = TMath::Abs(cluster->GetCellAbsId(j));
 702     fGeom->GetGlobal(id,pos);
 703     Double_t cellen = cells->GetCellAmplitude(id);
 704     Xc  += cellen*pos.X();
 705     Yc  += cellen*pos.Y();
 706     Sxx += cellen*pos.X()*pos.X();
 707     Syy += cellen*pos.Y()*pos.Y();
 708     Sxy += cellen*pos.X()*pos.Y();
 709   }
 710   Xc  /= Ec;
 711   Yc  /= Ec;
 712   Sxx /= Ec;
 713   Syy /= Ec;
 714   Sxy /= Ec;
 715   Sxx -= Xc*Xc;
 716   Syy -= Yc*Yc;
 717   Sxy -= Xc*Yc;
 718   sigmaMax = (Sxx + Syy + TMath::Sqrt((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy))/2.0;
 719   sigmaMax = TMath::Sqrt(sigmaMax);
 720   sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy))/2.0;
 721   sigmaMin = TMath::Sqrt(sigmaMin);
 722 }
 723