[u/mrichter/AliRoot.git] / PWG4 / CaloCalib / AliAnalysisTaskEMCALPi0PbPb.cxx

// $Id$

#include "TChain.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH1I.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TLorentzVector.h"
#include "AliESDVertex.h"

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"

#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliESDInputHandler.h"

#include "AliAnalysisTaskEMCALPi0PbPb.h"
#include "AliCentrality.h"

#include "AliEMCALGeoUtils.h"
// example of an analysis task creating a Neutral detector QA Histogram 

ClassImp(AliAnalysisTaskEMCALPi0PbPb)
//________________________________________________________________________
AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name) 
  : AliAnalysisTaskSE(name),
  fOutputDataESD(0),
  fhNEvents(0),
  fhCentVsClusMultEMC(0x0),
  fhCentVsCellMultEMC(0x0),
  fhCentVsClusMultEMCAll(0),
  fhCentVsCellMultEMCAll(0),
  fhCentVsInMVsPtEMC(0),
  fhCentVsInMVsPhiEMC(0),
  fhAsyVsPt(0),
  fhCentVsColuRowEMC(0x0),
  fhCentVsColuRowEnerEMC(0x0),
  fhMggPtEr(0x0),
  fhPtEnSg(0x0)   
{
  // Constructor

  nModuleEMC = 4,
  // Define input and output slots here
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  // Output slot #0 id reserved by the base class for AOD
  // Output slot #1 writes into a TH1 container
  DefineOutput(1, TList::Class());
  // Initialize the PHOS geometry
}
//________________________________________________________________________
AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb() {

  if(fOutputDataESD) delete fOutputDataESD;
  fOutputDataESD = 0;

}
//________________________________________________________________________
void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects()
{
  // Create histograms
  // Called once
 
  OpenFile(1); 
  fOutputDataESD = new TList();
  fOutputDataESD->SetOwner();
  
  fhNEvents   = new TH1I ("fhNEvents", "Number of events analyzed",1,0.,1.);
  fOutputDataESD-> Add(fhNEvents);

  fhCentVsClusMultEMC = new TH2F *[nModuleEMC]; 
  fhCentVsCellMultEMC = new TH2F *[nModuleEMC];
  fhCentVsColuRowEMC      = new TH3F *[nModuleEMC];
  fhCentVsColuRowEnerEMC  = new TH3F *[nModuleEMC];

  char key[256];
  char title[256];

    fhCentVsClusMultEMCAll = new TH2F("fhCentVsClusMultEMCAll","EMCal Cluster Multiplicity Dis. at centrality",500, 0, 500, 100, 0., 100.);
   fhCentVsCellMultEMCAll = new TH2F ("fhCentVsCellMultEMCAll","EMCal Cell Multiplicity Dis. at centrality",2000, 0, 2000, 100, 0., 100.);
    fOutputDataESD->Add(fhCentVsClusMultEMCAll);
    fOutputDataESD->Add(fhCentVsCellMultEMCAll);

  fhCentVsInMVsPtEMC =  new TH3F("fhCentVsInMVsPtEMC","EMCal Inva Mass Vs Pt Vs central 2 cluster",500, 0, 1, 400, 0, 40, 100, 0, 100);
  fhCentVsInMVsPhiEMC = new TH3F("fhCentVsInMVsPhiEMC","EMCal Inva Mass Vs Phi Vs central 2 cluster", 500, 0, 1, 360, 0, 7, 100, 0, 100); 
  fhAsyVsPt = new TH2F("fhAsyVsPt", "Asymmetry Vs Pt 2 cluster",100, 0, 1, 400, 0, 40);
  
  fhMggPtEr = new TH3F("fhMggPtEr","Invariant mass vs pT vs Energy ratio",100,0,1,100,0,20,100,0,1);
  fhMggPtEr->GetXaxis()->SetTitle("M_{#gamma#gamma} [GeV/c^{2}_{}]");
  fhMggPtEr->GetYaxis()->SetTitle("p_{T} [GeV/c]");
  fhMggPtEr->GetZaxis()->SetTitle("E_{r}");
  fhPtEnSg  = new TH3F("fhPtEnSg","photon Pt vs Energy vs Large axis",100,0,20,100,0,50,100,0,50); 
  fhPtEnSg->GetXaxis()->SetTitle("p_{T} [GeV/c]");
  fhPtEnSg->GetYaxis()->SetTitle("E [GeV]");
  fhPtEnSg->GetZaxis()->SetTitle("E#sigma_{max} [GeV]");

 fOutputDataESD->Add(fhCentVsInMVsPtEMC);
 fOutputDataESD->Add(fhCentVsInMVsPhiEMC);
 fOutputDataESD->Add(fhAsyVsPt);
 fOutputDataESD->Add(fhMggPtEr);
 fOutputDataESD->Add(fhPtEnSg);

  for(Int_t i=0; i<nModuleEMC; i++) { // for  EMCal
    // for cluster
    sprintf(key,"fhCentVsClusMultEMC_Modu%d",i);
    sprintf(title,"EMCal Module %d Cluster Multiplicity Dis. at centrality",i);
    fhCentVsClusMultEMC[i] = new TH2F (key,title, 500, 0,500, 100, 0.,100.);
   // for Cell
   sprintf(key,"fhCentVsCellMultEMC_Modu%d",i);
   sprintf(title,"EMCal Module %d Cell Multiplicity Dis. at centrality",i);
   fhCentVsCellMultEMC[i] = new TH2F (key,title, 2000, 0, 2000, 100, 0., 100.);
  // for occupancy
   sprintf(key, "fhCentVsColuRowEMC_Modu%d",i); 
   sprintf(title,"entries of cell in Module %d", i); 
   fhCentVsColuRowEMC[i] = new TH3F (key,title, 49, 0, 49, 25, 0, 25, 100, 0., 100.);
  //for energy 
   sprintf(key, "fhCentVsColuRowEMC_Modu%d",i);
   sprintf(title,"energy of cell in Module %d", i);
   fhCentVsColuRowEnerEMC[i] = new TH3F (key,title, 49, 0, 49, 25, 0, 25, 100, 0., 100.);
 
  fOutputDataESD->Add(fhCentVsClusMultEMC[i]);
  fOutputDataESD->Add(fhCentVsCellMultEMC[i]);
  fOutputDataESD->Add(fhCentVsColuRowEMC[i]);
  fOutputDataESD->Add(fhCentVsColuRowEnerEMC[i]);
  } 


 PostData(1, fOutputDataESD); 
}

//________________________________________________________________________
void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *) 
{
  // Main loop
  // Called for each event
      AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
      TRefArray * clusterListESD = new TRefArray();
      event->GetEMCALClusters(clusterListESD); 
      TClonesArray * clusterListAOD = new TClonesArray();
      if(AODEvent()){clusterListAOD = dynamic_cast<TClonesArray*> (AODEvent()->FindListObject("newEMCALClusters"));}

/*
  const AliESDVertex *vertex = event->GetPrimaryVertex();
    if (!vertex)   return;  
    if(! vertex->GetStatus()) return;
    if(TMath::Abs(vertex->GetZv())>10) return;

  Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);// | AliVEvent::kMBNoTRD);

  TString trigClasses = event->GetFiredTriggerClasses();

  if(trigClasses.Contains("CINT1") || trigClasses.Contains("CSH")  || trigClasses.Contains("CMU") || trigClasses.Contains("CBEAM") ) {
     printf("You are analyzing the REAL DATA!!! \n");
     if(!isSelected) return;

  }
*/  

  fhNEvents->Fill(0);

  // centrality
  AliCentrality *centrality =(AliCentrality *) event->GetCentrality();
   Float_t opt = (Float_t)centrality->GetCentralityPercentile("V0M");


  AliEMCALGeoUtils * fGeom   = new AliEMCALGeoUtils("EMCAL_FIRSTYEARV1","EMCAL");
    fGeom->GetEMCGeometry();

  //_____________ for EMCal  cell loop________________________________________________
  AliESDCaloCells *CellEMCAL   = event->GetEMCALCells();
   Int_t multCells = CellEMCAL ->GetNumberOfCells();
   fhCentVsCellMultEMCAll->Fill(multCells,opt);
   Int_t CellModuCount[4]={0,0,0,0};
  for(Int_t icell=0; icell<multCells; icell++){
    Int_t AbsID = CellEMCAL->GetCellNumber(icell); 
//    cout<<"EMCalAbsId =   "<<AbsID<<endl;  
    Int_t iSM1=-1, iTower1=-1, Iphi1=-1, Ieta1=-1;
    fGeom->GetCellIndex(AbsID, iSM1, iTower1, Iphi1, Ieta1 ); 
    Int_t iPhi1=-1, iEta1=-1;
    fGeom->GetCellPhiEtaIndexInSModule(iSM1, iTower1, Iphi1, Ieta1, iPhi1,iEta1);
    for(Int_t i=0; i<nModuleEMC; i++){
      if(i==iSM1){
        CellModuCount[i]++; 
        fhCentVsColuRowEMC[i]->Fill(iEta1,iPhi1,opt,1);// if hit on or not
        fhCentVsColuRowEnerEMC[i]->Fill(iEta1,iPhi1,opt,CellEMCAL->GetCellAmplitude(AbsID));// if hit on or not
      }
    }
  }//cell
    for(Int_t i=0; i<nModuleEMC; i++){
   fhCentVsCellMultEMC[i]->Fill(CellModuCount[i],opt); 
    }  

//__________for EMCal  cluster loop______________________________________________
  Int_t multCaloClust = 0;
  if(!AODEvent()) { multCaloClust = clusterListESD->GetEntriesFast(); 
  }else{            multCaloClust = clusterListAOD->GetEntriesFast(); }
  Double_t vtx0[3] = {0,0,0};
//  Double_t ClustEnergy = 0.;
  Int_t multEMCaClust =0;
  Int_t ClustModuCount[4] ={0,0,0,0};
  
  if(multCaloClust>1){  
  for(Int_t i1=0; i1<multCaloClust-1; i1++){
    
      multEMCaClust++ ; 
      TLorentzVector pemc1;
   if(!AODEvent()){ ((AliESDCaloCluster*) clusterListESD->At(i1))->GetMomentum(pemc1,vtx0);
   }else{           ((AliAODCaloCluster*) clusterListAOD->At(i1))->GetMomentum(pemc1,vtx0);}
// cluster showershape
    Double_t qua1;
    if(!AODEvent()){
      AliESDCaloCluster * ESDcluster = (AliESDCaloCluster*) clusterListESD->At(i1);
      Double_t sigmaMax1 = AliAnalysisTaskEMCALPi0PbPb::GetSigmaMax(ESDcluster, event);
      fhPtEnSg->Fill(pemc1.Pt(), pemc1.E(), sigmaMax1*pemc1.E());
      for(Int_t j=0; j<ESDcluster->GetNCells();j++){
       qua1<event->GetEMCALCells()->GetCellAmplitude(ESDcluster->GetCellAbsId(j))/ESDcluster->E() ? qua1 = event->GetEMCALCells()->GetCellAmplitude(ESDcluster->GetCellAbsId(j))/ESDcluster->E() : qua1 = qua1;
      }
    }else{
      AliAODCaloCluster * AODcluster = (AliAODCaloCluster*) clusterListAOD->At(i1);
      Double_t sigmaMax1 = GetSigmaMax(AODcluster, AODEvent());
      fhPtEnSg->Fill(pemc1.Pt(), pemc1.E(), sigmaMax1*pemc1.E());
      for(Int_t j=0; j<AODcluster->GetNCells();j++){
       qua1<AODEvent()->GetEMCALCells()->GetCellAmplitude(AODcluster->GetCellAbsId(j))/AODcluster->E() ? qua1 = AODEvent()->GetEMCALCells()->GetCellAmplitude(AODcluster->GetCellAbsId(j))/AODcluster->E() : qua1 = qua1;
      }
    }

    Int_t cellAbsId1=-1;
    fGeom->GetAbsCellIdFromEtaPhi(pemc1.Eta(), pemc1.Phi(), cellAbsId1);
    Int_t iSM1=-1, iTower1=-1, Iphi1=-1, Ieta1=-1;
    fGeom->GetCellIndex(cellAbsId1, iSM1, iTower1, Iphi1, Ieta1 );
      for(Int_t i=0; i<nModuleEMC; i++){
        if(i==iSM1){
          ClustModuCount[i]++;
        }
      }
  for(Int_t i2=i1; i2<multCaloClust; i2++){
   TLorentzVector pemc2;
   if(!AODEvent()){ ((AliESDCaloCluster*)clusterListESD->At(i2))->GetMomentum(pemc2,vtx0);
   }else{           ((AliAODCaloCluster*)clusterListAOD->At(i2))->GetMomentum(pemc2,vtx0);}
// cluster showershape
    Double_t qua2;
    if(!AODEvent()){
      AliESDCaloCluster * ESDcluster = (AliESDCaloCluster*) clusterListESD->At(i2);
      for(Int_t j=0; j<ESDcluster->GetNCells();j++){
       qua2<event->GetEMCALCells()->GetCellAmplitude(ESDcluster->GetCellAbsId(j))/ESDcluster->E() ? qua2 = event->GetEMCALCells()->GetCellAmplitude(ESDcluster->GetCellAbsId(j))/ESDcluster->E() : qua2 = qua2;
      }
    }else{
      AliAODCaloCluster * AODcluster = (AliAODCaloCluster*) clusterListAOD->At(i2);
      for(Int_t j=0; j<AODcluster->GetNCells();j++){
       qua2<AODEvent()->GetEMCALCells()->GetCellAmplitude(AODcluster->GetCellAbsId(j))/AODcluster->E() ? qua2 = AODEvent()->GetEMCALCells()->GetCellAmplitude(AODcluster->GetCellAbsId(j))/AODcluster->E() : qua2 = qua2;
      }
    }

   TLorentzVector pemc12  = pemc1  + pemc2;
   Double_t asy  = TMath::Abs((pemc1.E()-pemc2.E())/(pemc1.E()+pemc2.E()));
   Double_t qua = TMath::Max(qua1,qua2);
//   fhEMCALAsymPt->Fill(asy, p12.Pt());
   fhCentVsInMVsPtEMC->Fill(pemc12.M(), pemc12.Pt(), opt);
   fhCentVsInMVsPhiEMC->Fill(pemc12.M(), pemc12.Phi(), opt); 
   fhAsyVsPt->Fill(asy, pemc12.Pt());
   fhMggPtEr->Fill(pemc12.M(), pemc12.Pt(), qua);

  }//loop second cluster
   }// EMCal first cluster loop
  } //multCaloCluster>1

     for(Int_t i=0; i<nModuleEMC; i++){ // fill Module
          fhCentVsClusMultEMC[i]->Fill(ClustModuCount[i],opt);
      }
    fhCentVsClusMultEMCAll->Fill(multEMCaClust,opt);

   
// Post output data.
  PostData(1, fOutputDataESD);
}      

//________________________________________________________________________
void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *) 
{
  // Draw result to the screen
  // Called once at the end of the query
}
//________________________________________________________________________
Double_t AliAnalysisTaskEMCALPi0PbPb::GetSigmaMax(AliESDCaloCluster * cluster, AliESDEvent * event){
 Double_t sigmaMax;
 Double_t Ec = cluster->E();
 Double_t Xc = 0 ;
 Double_t Yc = 0 ;
 Double_t Sxx = 0 ;
 Double_t Sxy = 0 ;
 Double_t Syy = 0 ;
 AliEMCALGeoUtils * fGeom   = new AliEMCALGeoUtils("EMCAL_FIRSTYEARV1","EMCAL");
 fGeom->GetEMCGeometry();
 for(Int_t j=0; j<cluster->GetNCells();j++){
   TVector3 pos;
   fGeom->GetGlobal(cluster->GetCellAbsId(j),pos);
   Xc = Xc + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))*pos.X()/Ec; 
   Yc = Yc + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))*pos.Y()/Ec;    
 }
 for(Int_t j=0; j<cluster->GetNCells();j++){
   TVector3 pos;
   fGeom->GetGlobal(cluster->GetCellAbsId(j),pos);
   Sxx = Sxx + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))*(pos.X()-Xc)*(pos.X()-Xc)/Ec; 
   Syy = Syy + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))*(pos.Y()-Yc)*(pos.Y()-Yc)/Ec; 
   Sxy = Sxy + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))*(pos.X()-Xc)*(pos.Y()-Yc)/Ec; 
 } 
 sigmaMax = (Sxx + Syy + TMath::Sqrt((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy))/2.0;
 sigmaMax = TMath::Sqrt(sigmaMax); 
 return  sigmaMax;
}
//________________________________________________________________________
Double_t AliAnalysisTaskEMCALPi0PbPb::GetSigmaMax(AliAODCaloCluster * cluster, AliAODEvent * event){
 Double_t sigmaMax;
 Double_t Ec = cluster->E();
 Double_t Xc = 0 ;
 Double_t Yc = 0 ;
 Double_t Sxx = 0 ;
 Double_t Sxy = 0 ;
 Double_t Syy = 0 ;
 AliEMCALGeoUtils * fGeom   = new AliEMCALGeoUtils("EMCAL_FIRSTYEARV1","EMCAL");
 fGeom->GetEMCGeometry();
 for(Int_t j=0; j<cluster->GetNCells();j++){
   TVector3 pos;
   fGeom->GetGlobal(cluster->GetCellAbsId(j),pos);
   Xc = Xc + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))*pos.X()/Ec; 
   Yc = Yc + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))*pos.Y()/Ec;    
 }
 for(Int_t j=0; j<cluster->GetNCells();j++){
   TVector3 pos;
   fGeom->GetGlobal(cluster->GetCellAbsId(j),pos);
   Sxx = Sxx + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))*(pos.X()-Xc)*(pos.X()-Xc)/Ec; 
   Syy = Syy + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))*(pos.Y()-Yc)*(pos.Y()-Yc)/Ec; 
   Sxy = Sxy + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))*(pos.X()-Xc)*(pos.Y()-Yc)/Ec; 
 } 
 sigmaMax = (Sxx + Syy + TMath::Sqrt((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy))/2.0;
 sigmaMax = TMath::Sqrt(sigmaMax); 
 return  sigmaMax;

}
Commit	Line	Data
ea3fd2d5	1	// $Id$
	2
	3	#include "TChain.h"
	4	#include "TTree.h"
	5	#include "TH1F.h"
	6	#include "TH1I.h"
	7	#include "TH2F.h"
	8	#include "TCanvas.h"
	9	#include "TLorentzVector.h"
	10	#include "AliESDVertex.h"
	11
	12	#include "AliAnalysisTask.h"
	13	#include "AliAnalysisManager.h"
	14
	15	#include "AliESDEvent.h"
	16	#include "AliAODEvent.h"
	17	#include "AliESDInputHandler.h"
	18
	19	#include "AliAnalysisTaskEMCALPi0PbPb.h"
	20	#include "AliCentrality.h"
	21
	22	#include "AliEMCALGeoUtils.h"
	23	// example of an analysis task creating a Neutral detector QA Histogram
	24
	25	ClassImp(AliAnalysisTaskEMCALPi0PbPb)
	26	//________________________________________________________________________
	27	AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name)
	28	: AliAnalysisTaskSE(name),
	29	fOutputDataESD(0),
	30	fhNEvents(0),
	31	fhCentVsClusMultEMC(0x0),
	32	fhCentVsCellMultEMC(0x0),
	33	fhCentVsClusMultEMCAll(0),
	34	fhCentVsCellMultEMCAll(0),
	35	fhCentVsInMVsPtEMC(0),
	36	fhCentVsInMVsPhiEMC(0),
	37	fhAsyVsPt(0),
	38	fhCentVsColuRowEMC(0x0),
	39	fhCentVsColuRowEnerEMC(0x0),
	40	fhMggPtEr(0x0),
	41	fhPtEnSg(0x0)
	42	{
	43	// Constructor
	44
	45	nModuleEMC = 4,
	46	// Define input and output slots here
	47	// Input slot #0 works with a TChain
	48	DefineInput(0, TChain::Class());
	49	// Output slot #0 id reserved by the base class for AOD
	50	// Output slot #1 writes into a TH1 container
	51	DefineOutput(1, TList::Class());
	52	// Initialize the PHOS geometry
	53	}
	54	//________________________________________________________________________
	55	AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb() {
	56
	57	if(fOutputDataESD) delete fOutputDataESD;
	58	fOutputDataESD = 0;
	59
	60	}
	61	//________________________________________________________________________
	62	void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects()
	63	{
	64	// Create histograms
65	// Called once
66
67	OpenFile(1);
68	fOutputDataESD = new TList();
69	fOutputDataESD->SetOwner();
70
71	fhNEvents = new TH1I ("fhNEvents", "Number of events analyzed",1,0.,1.);
72	fOutputDataESD-> Add(fhNEvents);
73
74	fhCentVsClusMultEMC = new TH2F *[nModuleEMC];
75	fhCentVsCellMultEMC = new TH2F *[nModuleEMC];
76	fhCentVsColuRowEMC = new TH3F *[nModuleEMC];
77	fhCentVsColuRowEnerEMC = new TH3F *[nModuleEMC];
78
79	char key[256];
80	char title[256];
81
82	fhCentVsClusMultEMCAll = new TH2F("fhCentVsClusMultEMCAll","EMCal Cluster Multiplicity Dis. at centrality",500, 0, 500, 100, 0., 100.);
83	fhCentVsCellMultEMCAll = new TH2F ("fhCentVsCellMultEMCAll","EMCal Cell Multiplicity Dis. at centrality",2000, 0, 2000, 100, 0., 100.);
84	fOutputDataESD->Add(fhCentVsClusMultEMCAll);
85	fOutputDataESD->Add(fhCentVsCellMultEMCAll);
86
87	fhCentVsInMVsPtEMC = new TH3F("fhCentVsInMVsPtEMC","EMCal Inva Mass Vs Pt Vs central 2 cluster",500, 0, 1, 400, 0, 40, 100, 0, 100);
88	fhCentVsInMVsPhiEMC = new TH3F("fhCentVsInMVsPhiEMC","EMCal Inva Mass Vs Phi Vs central 2 cluster", 500, 0, 1, 360, 0, 7, 100, 0, 100);
89	fhAsyVsPt = new TH2F("fhAsyVsPt", "Asymmetry Vs Pt 2 cluster",100, 0, 1, 400, 0, 40);
90
91	fhMggPtEr = new TH3F("fhMggPtEr","Invariant mass vs pT vs Energy ratio",100,0,1,100,0,20,100,0,1);
92	fhMggPtEr->GetXaxis()->SetTitle("M_{#gamma#gamma} [GeV/c^{2}_{}]");
93	fhMggPtEr->GetYaxis()->SetTitle("p_{T} [GeV/c]");
94	fhMggPtEr->GetZaxis()->SetTitle("E_{r}");
95	fhPtEnSg = new TH3F("fhPtEnSg","photon Pt vs Energy vs Large axis",100,0,20,100,0,50,100,0,50);
96	fhPtEnSg->GetXaxis()->SetTitle("p_{T} [GeV/c]");
97	fhPtEnSg->GetYaxis()->SetTitle("E [GeV]");
98	fhPtEnSg->GetZaxis()->SetTitle("E#sigma_{max} [GeV]");
99
100	fOutputDataESD->Add(fhCentVsInMVsPtEMC);
101	fOutputDataESD->Add(fhCentVsInMVsPhiEMC);
102	fOutputDataESD->Add(fhAsyVsPt);
103	fOutputDataESD->Add(fhMggPtEr);
104	fOutputDataESD->Add(fhPtEnSg);
105
106	for(Int_t i=0; i<nModuleEMC; i++) { // for EMCal
107	// for cluster
108	sprintf(key,"fhCentVsClusMultEMC_Modu%d",i);
109	sprintf(title,"EMCal Module %d Cluster Multiplicity Dis. at centrality",i);
110	fhCentVsClusMultEMC[i] = new TH2F (key,title, 500, 0,500, 100, 0.,100.);
111	// for Cell
112	sprintf(key,"fhCentVsCellMultEMC_Modu%d",i);
113	sprintf(title,"EMCal Module %d Cell Multiplicity Dis. at centrality",i);
114	fhCentVsCellMultEMC[i] = new TH2F (key,title, 2000, 0, 2000, 100, 0., 100.);
115	// for occupancy
116	sprintf(key, "fhCentVsColuRowEMC_Modu%d",i);
117	sprintf(title,"entries of cell in Module %d", i);
118	fhCentVsColuRowEMC[i] = new TH3F (key,title, 49, 0, 49, 25, 0, 25, 100, 0., 100.);
119	//for energy
120	sprintf(key, "fhCentVsColuRowEMC_Modu%d",i);
121	sprintf(title,"energy of cell in Module %d", i);
122	fhCentVsColuRowEnerEMC[i] = new TH3F (key,title, 49, 0, 49, 25, 0, 25, 100, 0., 100.);
123
124	fOutputDataESD->Add(fhCentVsClusMultEMC[i]);
125	fOutputDataESD->Add(fhCentVsCellMultEMC[i]);
126	fOutputDataESD->Add(fhCentVsColuRowEMC[i]);
127	fOutputDataESD->Add(fhCentVsColuRowEnerEMC[i]);
128	}
129
130
131	PostData(1, fOutputDataESD);
132	}
133
134	//________________________________________________________________________
135	void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *)
136	{
137	// Main loop
138	// Called for each event
139	AliESDEvent event = dynamic_cast<AliESDEvent>(InputEvent());
140	TRefArray * clusterListESD = new TRefArray();
141	event->GetEMCALClusters(clusterListESD);
142	TClonesArray * clusterListAOD = new TClonesArray();
143	if(AODEvent()){clusterListAOD = dynamic_cast<TClonesArray*> (AODEvent()->FindListObject("newEMCALClusters"));}
144
145	/*
146	const AliESDVertex *vertex = event->GetPrimaryVertex();
147	if (!vertex) return;
148	if(! vertex->GetStatus()) return;
149	if(TMath::Abs(vertex->GetZv())>10) return;
150
151	Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);// \| AliVEvent::kMBNoTRD);
152
153	TString trigClasses = event->GetFiredTriggerClasses();
154
155	if(trigClasses.Contains("CINT1") \|\| trigClasses.Contains("CSH") \|\| trigClasses.Contains("CMU") \|\| trigClasses.Contains("CBEAM") ) {
156	printf("You are analyzing the REAL DATA!!! \n");
157	if(!isSelected) return;
158
159	}
160	*/
161
162	fhNEvents->Fill(0);
163
164	// centrality
165	AliCentrality centrality =(AliCentrality ) event->GetCentrality();
166	Float_t opt = (Float_t)centrality->GetCentralityPercentile("V0M");
167
168
169
170
171	AliEMCALGeoUtils * fGeom = new AliEMCALGeoUtils("EMCAL_FIRSTYEARV1","EMCAL");
172	fGeom->GetEMCGeometry();
173
174	//_____________ for EMCal cell loop________________________________________________
175	AliESDCaloCells *CellEMCAL = event->GetEMCALCells();
176	Int_t multCells = CellEMCAL ->GetNumberOfCells();
177	fhCentVsCellMultEMCAll->Fill(multCells,opt);
178	Int_t CellModuCount[4]={0,0,0,0};
179	for(Int_t icell=0; icell<multCells; icell++){
180	Int_t AbsID = CellEMCAL->GetCellNumber(icell);
181	// cout<<"EMCalAbsId = "<<AbsID<<endl;
182	Int_t iSM1=-1, iTower1=-1, Iphi1=-1, Ieta1=-1;
183	fGeom->GetCellIndex(AbsID, iSM1, iTower1, Iphi1, Ieta1 );
184	Int_t iPhi1=-1, iEta1=-1;
185	fGeom->GetCellPhiEtaIndexInSModule(iSM1, iTower1, Iphi1, Ieta1, iPhi1,iEta1);
186	for(Int_t i=0; i<nModuleEMC; i++){
187	if(i==iSM1){
188	CellModuCount[i]++;
189	fhCentVsColuRowEMC[i]->Fill(iEta1,iPhi1,opt,1);// if hit on or not
190	fhCentVsColuRowEnerEMC[i]->Fill(iEta1,iPhi1,opt,CellEMCAL->GetCellAmplitude(AbsID));// if hit on or not
191	}
192	}
193	}//cell
194	for(Int_t i=0; i<nModuleEMC; i++){
195	fhCentVsCellMultEMC[i]->Fill(CellModuCount[i],opt);
196	}
197
198	//__________for EMCal cluster loop______________________________________________
199	Int_t multCaloClust = 0;
200	if(!AODEvent()) { multCaloClust = clusterListESD->GetEntriesFast();
201	}else{ multCaloClust = clusterListAOD->GetEntriesFast(); }
202	Double_t vtx0[3] = {0,0,0};
203	// Double_t ClustEnergy = 0.;
204	Int_t multEMCaClust =0;
205	Int_t ClustModuCount[4] ={0,0,0,0};
206
207	if(multCaloClust>1){
208	for(Int_t i1=0; i1<multCaloClust-1; i1++){
209
210	multEMCaClust++ ;
211	TLorentzVector pemc1;
212	if(!AODEvent()){ ((AliESDCaloCluster*) clusterListESD->At(i1))->GetMomentum(pemc1,vtx0);
213	}else{ ((AliAODCaloCluster*) clusterListAOD->At(i1))->GetMomentum(pemc1,vtx0);}
214	// cluster showershape
215	Double_t qua1;
216	if(!AODEvent()){
217	AliESDCaloCluster * ESDcluster = (AliESDCaloCluster*) clusterListESD->At(i1);
218	Double_t sigmaMax1 = AliAnalysisTaskEMCALPi0PbPb::GetSigmaMax(ESDcluster, event);
219	fhPtEnSg->Fill(pemc1.Pt(), pemc1.E(), sigmaMax1*pemc1.E());
220	for(Int_t j=0; j<ESDcluster->GetNCells();j++){
221	qua1<event->GetEMCALCells()->GetCellAmplitude(ESDcluster->GetCellAbsId(j))/ESDcluster->E() ? qua1 = event->GetEMCALCells()->GetCellAmplitude(ESDcluster->GetCellAbsId(j))/ESDcluster->E() : qua1 = qua1;
222	}
223	}else{
224	AliAODCaloCluster * AODcluster = (AliAODCaloCluster*) clusterListAOD->At(i1);
225	Double_t sigmaMax1 = GetSigmaMax(AODcluster, AODEvent());
226	fhPtEnSg->Fill(pemc1.Pt(), pemc1.E(), sigmaMax1*pemc1.E());
227	for(Int_t j=0; j<AODcluster->GetNCells();j++){
228	qua1<AODEvent()->GetEMCALCells()->GetCellAmplitude(AODcluster->GetCellAbsId(j))/AODcluster->E() ? qua1 = AODEvent()->GetEMCALCells()->GetCellAmplitude(AODcluster->GetCellAbsId(j))/AODcluster->E() : qua1 = qua1;
229	}
230	}
231
232	Int_t cellAbsId1=-1;
233	fGeom->GetAbsCellIdFromEtaPhi(pemc1.Eta(), pemc1.Phi(), cellAbsId1);
234	Int_t iSM1=-1, iTower1=-1, Iphi1=-1, Ieta1=-1;
235	fGeom->GetCellIndex(cellAbsId1, iSM1, iTower1, Iphi1, Ieta1 );
236	for(Int_t i=0; i<nModuleEMC; i++){
237	if(i==iSM1){
238	ClustModuCount[i]++;
239	}
240	}
241	for(Int_t i2=i1; i2<multCaloClust; i2++){
242	TLorentzVector pemc2;
243	if(!AODEvent()){ ((AliESDCaloCluster*)clusterListESD->At(i2))->GetMomentum(pemc2,vtx0);
244	}else{ ((AliAODCaloCluster*)clusterListAOD->At(i2))->GetMomentum(pemc2,vtx0);}
245	// cluster showershape
246	Double_t qua2;
247	if(!AODEvent()){
248	AliESDCaloCluster * ESDcluster = (AliESDCaloCluster*) clusterListESD->At(i2);
249	for(Int_t j=0; j<ESDcluster->GetNCells();j++){
250	qua2<event->GetEMCALCells()->GetCellAmplitude(ESDcluster->GetCellAbsId(j))/ESDcluster->E() ? qua2 = event->GetEMCALCells()->GetCellAmplitude(ESDcluster->GetCellAbsId(j))/ESDcluster->E() : qua2 = qua2;
251	}
252	}else{
253	AliAODCaloCluster * AODcluster = (AliAODCaloCluster*) clusterListAOD->At(i2);
254	for(Int_t j=0; j<AODcluster->GetNCells();j++){
255	qua2<AODEvent()->GetEMCALCells()->GetCellAmplitude(AODcluster->GetCellAbsId(j))/AODcluster->E() ? qua2 = AODEvent()->GetEMCALCells()->GetCellAmplitude(AODcluster->GetCellAbsId(j))/AODcluster->E() : qua2 = qua2;
256	}
257	}
258
259	TLorentzVector pemc12 = pemc1 + pemc2;
260	Double_t asy = TMath::Abs((pemc1.E()-pemc2.E())/(pemc1.E()+pemc2.E()));
261	Double_t qua = TMath::Max(qua1,qua2);
262	// fhEMCALAsymPt->Fill(asy, p12.Pt());
263	fhCentVsInMVsPtEMC->Fill(pemc12.M(), pemc12.Pt(), opt);
264	fhCentVsInMVsPhiEMC->Fill(pemc12.M(), pemc12.Phi(), opt);
265	fhAsyVsPt->Fill(asy, pemc12.Pt());
266	fhMggPtEr->Fill(pemc12.M(), pemc12.Pt(), qua);
267
268	}//loop second cluster
269	}// EMCal first cluster loop
270	} //multCaloCluster>1
271
272	for(Int_t i=0; i<nModuleEMC; i++){ // fill Module
273	fhCentVsClusMultEMC[i]->Fill(ClustModuCount[i],opt);
274	}
275	fhCentVsClusMultEMCAll->Fill(multEMCaClust,opt);
276
277
278	// Post output data.
279	PostData(1, fOutputDataESD);
280	}
281
282	//________________________________________________________________________
283	void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *)
284	{
285	// Draw result to the screen
286	// Called once at the end of the query
287	}
288	//________________________________________________________________________
289	Double_t AliAnalysisTaskEMCALPi0PbPb::GetSigmaMax(AliESDCaloCluster * cluster, AliESDEvent * event){
290	Double_t sigmaMax;
291	Double_t Ec = cluster->E();
292	Double_t Xc = 0 ;
293	Double_t Yc = 0 ;
294	Double_t Sxx = 0 ;
295	Double_t Sxy = 0 ;
296	Double_t Syy = 0 ;
297	AliEMCALGeoUtils * fGeom = new AliEMCALGeoUtils("EMCAL_FIRSTYEARV1","EMCAL");
298	fGeom->GetEMCGeometry();
299	for(Int_t j=0; j<cluster->GetNCells();j++){
300	TVector3 pos;
301	fGeom->GetGlobal(cluster->GetCellAbsId(j),pos);
302	Xc = Xc + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))*pos.X()/Ec;
303	Yc = Yc + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))*pos.Y()/Ec;
304	}
305	for(Int_t j=0; j<cluster->GetNCells();j++){
306	TVector3 pos;
307	fGeom->GetGlobal(cluster->GetCellAbsId(j),pos);
308	Sxx = Sxx + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))(pos.X()-Xc)(pos.X()-Xc)/Ec;
309	Syy = Syy + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))(pos.Y()-Yc)(pos.Y()-Yc)/Ec;
310	Sxy = Sxy + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))(pos.X()-Xc)(pos.Y()-Yc)/Ec;
311	}
312	sigmaMax = (Sxx + Syy + TMath::Sqrt((Sxx-Syy)(Sxx-Syy)+4.0Sxy*Sxy))/2.0;
313	sigmaMax = TMath::Sqrt(sigmaMax);
314	return sigmaMax;
315	}
316	//________________________________________________________________________
317	Double_t AliAnalysisTaskEMCALPi0PbPb::GetSigmaMax(AliAODCaloCluster * cluster, AliAODEvent * event){
318	Double_t sigmaMax;
319	Double_t Ec = cluster->E();
320	Double_t Xc = 0 ;
321	Double_t Yc = 0 ;
322	Double_t Sxx = 0 ;
323	Double_t Sxy = 0 ;
324	Double_t Syy = 0 ;
325	AliEMCALGeoUtils * fGeom = new AliEMCALGeoUtils("EMCAL_FIRSTYEARV1","EMCAL");
326	fGeom->GetEMCGeometry();
327	for(Int_t j=0; j<cluster->GetNCells();j++){
328	TVector3 pos;
329	fGeom->GetGlobal(cluster->GetCellAbsId(j),pos);
330	Xc = Xc + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))*pos.X()/Ec;
331	Yc = Yc + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))*pos.Y()/Ec;
332	}
333	for(Int_t j=0; j<cluster->GetNCells();j++){
334	TVector3 pos;
335	fGeom->GetGlobal(cluster->GetCellAbsId(j),pos);
336	Sxx = Sxx + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))(pos.X()-Xc)(pos.X()-Xc)/Ec;
337	Syy = Syy + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))(pos.Y()-Yc)(pos.Y()-Yc)/Ec;
338	Sxy = Sxy + event->GetEMCALCells()->GetCellAmplitude(cluster->GetCellAbsId(j))(pos.X()-Xc)(pos.Y()-Yc)/Ec;
339	}
340	sigmaMax = (Sxx + Syy + TMath::Sqrt((Sxx-Syy)(Sxx-Syy)+4.0Sxy*Sxy))/2.0;
341	sigmaMax = TMath::Sqrt(sigmaMax);
342	return sigmaMax;
343
344	}