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

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: Boris Polishchuk                                               *
 * Adapted to AOD reading by Gustavo Conesa  *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//---------------------------------------------------------------------------// 
//                                                                           //
// Fill histograms (one per cell) with two-cluster invariant mass            //
// using calibration coefficients of the previous iteration.                 //
// Histogram for a given cell is filled if the most energy of one cluster    //
// is deposited in this cell and the other cluster could be anywherein EMCAL.//
//                                                                           //
//---------------------------------------------------------------------------//

//#include <cstdlib>
//#include <Riostream.h>
// Root 
#include "TLorentzVector.h"
//#include "TVector3.h"
#include "TRefArray.h"
#include "TList.h"
#include "TH1F.h"

// AliRoot
#include "AliAnalysisTaskEMCALPi0CalibSelection.h"
#include "AliAODEvent.h"
#include "AliESDEvent.h"
#include "AliESDCaloCluster.h"
#include "AliESDCaloCells.h"
#include "AliEMCALGeometry.h"
#include "AliAODCaloCluster.h"
#include "AliAODCaloCells.h"
#include "AliEMCALAodCluster.h"
#include "AliEMCALCalibData.h"

ClassImp(AliAnalysisTaskEMCALPi0CalibSelection)

//__________________________________________________
AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection() :
AliAnalysisTaskSE(),fEMCALGeo(0x0),fCalibData(0x0), fEmin(0.), fLogWeight(4.5), fCopyAOD(kFALSE), 
fEMCALGeoName("EMCAL_COMPLETE"), fOldData(kFALSE),fOutputContainer(0x0),fHmgg(0x0)
{
  //Default constructor.

  for(Int_t iMod=0; iMod < 12; iMod++) {
    for(Int_t iX=0; iX<24; iX++) {
      for(Int_t iZ=0; iZ<48; iZ++) {
	fHmpi0[iMod][iX][iZ]=0;
      }
    } 
  }
		
}

//__________________________________________________
AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(const char* name) :
  AliAnalysisTaskSE(name),fEMCALGeo(0x0),fCalibData(0x0), fEmin(0.), fLogWeight(4.5), fCopyAOD(kFALSE), 
  fEMCALGeoName("EMCAL_COMPLETE"), fOldData(kFALSE), fOutputContainer(0x0),fHmgg(0x0)
{
  //Named constructor which should be used.
  
  DefineOutput(1,TList::Class());
  
  for(Int_t iMod=0; iMod < 12; iMod++) {
    for(Int_t iX=0; iX<24; iX++) {
      for(Int_t iZ=0; iZ<48; iZ++) {
	fHmpi0[iMod][iX][iZ]=0;
      }
    } 
  }

}

//__________________________________________________
AliAnalysisTaskEMCALPi0CalibSelection::~AliAnalysisTaskEMCALPi0CalibSelection()
{
  //Destructor.
  
  if(fOutputContainer){
    fOutputContainer->Delete() ; 
    delete fOutputContainer ;
  }
	
  if(fCalibData)  delete fCalibData;
  if(fEMCALGeo)   delete fEMCALGeo;
	
}


//__________________________________________________
void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromAOD()
{
  // Copy AOD header, vertex, CaloClusters and CaloCells to output AOD
  
  AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
  
  // set arrays and pointers
  Float_t posF[3];
  Double_t pos[3];
  
  Double_t covVtx[6];
  
  for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
  
  // Update the header
  AliAODHeader*headerin = aod->GetHeader();
  AliAODHeader* header = AODEvent()->GetHeader();
  header->SetRunNumber(headerin->GetRunNumber());
  header->SetBunchCrossNumber(headerin->GetBunchCrossNumber());
  header->SetOrbitNumber(headerin->GetOrbitNumber());
  header->SetPeriodNumber(headerin->GetPeriodNumber());
  header->SetEventType(headerin->GetEventType());
  header->SetMuonMagFieldScale(headerin->GetMuonMagFieldScale());
  header->SetCentrality(headerin->GetCentrality()); 
  
  header->SetTriggerMask(headerin->GetTriggerMask()); 
  header->SetTriggerCluster(headerin->GetTriggerCluster());
  header->SetMagneticField(headerin->GetMagneticField());
  header->SetZDCN1Energy(headerin->GetZDCN1Energy());
  header->SetZDCP1Energy(headerin->GetZDCP1Energy());
  header->SetZDCN2Energy(headerin->GetZDCN2Energy());
  header->SetZDCP2Energy(headerin->GetZDCP2Energy());
  header->SetZDCEMEnergy(headerin->GetZDCEMEnergy(0),headerin->GetZDCEMEnergy(1));
  Float_t diamxy[2]={aod->GetDiamondX(),aod->GetDiamondY()};
  Float_t diamcov[3]; aod->GetDiamondCovXY(diamcov);
  header->SetDiamond(diamxy,diamcov);
  //
  //
  Int_t nVertices = 1 ;/* = prim. vtx*/;
  Int_t nCaloClus = aod->GetNCaloClusters();
  
  AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
  
  // Access to the AOD container of vertices
  TClonesArray &vertices = *(AODEvent()->GetVertices());
  Int_t jVertices=0;
  
  // Add primary vertex. The primary tracks will be defined
  // after the loops on the composite objects (V0, cascades, kinks)
  const AliAODVertex *vtx = aod->GetPrimaryVertex();
  
  vtx->GetXYZ(pos); // position
  vtx->GetCovMatrix(covVtx); //covariance matrix
  
  AliAODVertex * primary = new(vertices[jVertices++])
    AliAODVertex(pos, covVtx, vtx->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
  primary->SetName(vtx->GetName());
  primary->SetTitle(vtx->GetTitle());
  
  // Access to the AOD container of clusters
  TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
  Int_t jClusters=0;
  
  for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
    
    AliAODCaloCluster * cluster = aod->GetCaloCluster(iClust);
    
    //Check if it is a EMCAL cluster
    if(!cluster->IsEMCALCluster())  continue ;
    
    Int_t id       = cluster->GetID();
    Float_t energy = cluster->E();
    cluster->GetPosition(posF);
    Char_t ttype   = cluster->GetType(); 
    
    AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) 
      AliAODCaloCluster(id,
			0,
			0x0,
			energy,
			posF,
			NULL,
			ttype);
    
    caloCluster->SetCaloCluster(cluster->GetDistToBadChannel(),
				cluster->GetDispersion(),
				cluster->GetM20(), cluster->GetM02(),
				cluster->GetEmcCpvDistance(),  
				cluster->GetNExMax(),cluster->GetTOF()) ;
    
    caloCluster->SetPIDFromESD(cluster->PID());
    caloCluster->SetNCells(cluster->GetNCells());
    caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
	if(!fOldData){
		caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
	}
	else{
		Double_t *fraction = cluster->GetCellsAmplitudeFraction();
		for(Int_t i = 0; i < cluster->GetNCells() ; i++) fraction[i] = 1;
		caloCluster->SetCellsAmplitudeFraction(fraction);
	}
    
  } 

  caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters	 
  // end of loop on calo clusters
  
  // fill EMCAL cell info
  if (aod->GetEMCALCells()) { // protection against missing AOD information
	    AliAODCaloCells &aodinEMcells = *(aod->GetEMCALCells());
    Int_t nEMcell = aodinEMcells.GetNumberOfCells() ;
    
    AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
    aodEMcells.CreateContainer(nEMcell);
    aodEMcells.SetType(AliAODCaloCells::kEMCAL);
	
	Double_t calibFactor = 1;
	if(fOldData) calibFactor = 0.0153;

    for (Int_t iCell = 0; iCell < nEMcell; iCell++) {      
      aodEMcells.SetCell(iCell,aodinEMcells.GetCellNumber(iCell),aodinEMcells.GetAmplitude(iCell)*calibFactor);
    }
    aodEMcells.Sort();
	  
  }
  
}

//__________________________________________________
void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromESD()
{
  
  // Copy Header, Vertex, CaloClusters and CaloCells from ESDs to AODs
  
  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
  
  // set arrays and pointers
  Float_t posF[3];
  Double_t pos[3];
  
  Double_t covVtx[6];
  
  for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
  
  // Update the header
  
  AliAODHeader* header = AODEvent()->GetHeader();
  header->SetRunNumber(esd->GetRunNumber());
  header->SetBunchCrossNumber(esd->GetBunchCrossNumber());
  header->SetOrbitNumber(esd->GetOrbitNumber());
  header->SetPeriodNumber(esd->GetPeriodNumber());
  header->SetEventType(esd->GetEventType());
  header->SetMuonMagFieldScale(-999.); // FIXME
  header->SetCentrality(-999.);        // FIXME
  
  
  header->SetTriggerMask(esd->GetTriggerMask()); 
  header->SetTriggerCluster(esd->GetTriggerCluster());
  header->SetMagneticField(esd->GetMagneticField());
  header->SetZDCN1Energy(esd->GetZDCN1Energy());
  header->SetZDCP1Energy(esd->GetZDCP1Energy());
  header->SetZDCN2Energy(esd->GetZDCN2Energy());
  header->SetZDCP2Energy(esd->GetZDCP2Energy());
  header->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
  Float_t diamxy[2]={esd->GetDiamondX(),esd->GetDiamondY()};
  Float_t diamcov[3]; esd->GetDiamondCovXY(diamcov);
  header->SetDiamond(diamxy,diamcov);
  //
  //
  Int_t nVertices = 1 ;/* = prim. vtx*/;
  Int_t nCaloClus = esd->GetNumberOfCaloClusters();
  
  AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
  
  // Access to the AOD container of vertices
  TClonesArray &vertices = *(AODEvent()->GetVertices());
  Int_t jVertices=0;
  
  // Add primary vertex. The primary tracks will be defined
  // after the loops on the composite objects (V0, cascades, kinks)
  const AliESDVertex *vtx = esd->GetPrimaryVertex();
  
  vtx->GetXYZ(pos); // position
  vtx->GetCovMatrix(covVtx); //covariance matrix
  
  AliAODVertex * primary = new(vertices[jVertices++])
    AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
  primary->SetName(vtx->GetName());
  primary->SetTitle(vtx->GetTitle());
  
  // Access to the AOD container of clusters
  TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
  Int_t jClusters=0;
  
  for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
    
    AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
    
    //Check which calorimeter information we want to keep.
    if(!cluster->IsEMCAL())  continue ;
    
    Int_t id       = cluster->GetID();
    Float_t energy = cluster->E();
    cluster->GetPosition(posF);
    
    AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) 
      AliAODCaloCluster(id,
			0,
			0x0,
			energy,
			posF,
			NULL,
			AliAODCluster::kEMCALClusterv1);
    
    caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
				cluster->GetClusterDisp(),
				cluster->GetM20(), cluster->GetM02(),
				cluster->GetEmcCpvDistance(),  
				cluster->GetNExMax(),cluster->GetTOF()) ;
    
    caloCluster->SetPIDFromESD(cluster->GetPid());
    caloCluster->SetNCells(cluster->GetNCells());
    caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
	if(!fOldData){
		  caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
	}
	else{
		  Double_t *fraction = cluster->GetCellsAmplitudeFraction();
		  for(Int_t i = 0; i < cluster->GetNCells() ; i++) fraction[i] = 1;
		  caloCluster->SetCellsAmplitudeFraction(fraction);
	  }    
  } 
  caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters	 
  // end of loop on calo clusters
  
  // fill EMCAL cell info

  if( esd->GetEMCALCells()) { // protection against missing ESD information
    AliESDCaloCells &esdEMcells = *(esd->GetEMCALCells());
	Int_t nEMcell = esdEMcells.GetNumberOfCells() ;
    
    AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
    aodEMcells.CreateContainer(nEMcell);
    aodEMcells.SetType(AliAODCaloCells::kEMCAL);  
	  
	Double_t calibFactor = 1;
	if(fOldData) calibFactor = 0.0153;
   
	for (Int_t iCell = 0; iCell < nEMcell; iCell++) {      
      aodEMcells.SetCell(iCell,esdEMcells.GetCellNumber(iCell),esdEMcells.GetAmplitude(iCell)*calibFactor);
    }
    aodEMcells.Sort();
	  
  }

}

//__________________________________________________
void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
{
  //Create output container
  fOutputContainer = new TList();
  
  char hname[128], htitl[128];
  
	for(Int_t iMod=0; iMod < AliEMCALGeoParams::fgkEMCALModules; iMod++) {
    for(Int_t iX=0; iX<AliEMCALGeoParams::fgkEMCALCols; iX++) {
      for(Int_t iZ=0; iZ<AliEMCALGeoParams::fgkEMCALRows; iZ++) {
	sprintf(hname,"%d_%d_%d",iMod,iX,iZ);
	sprintf(htitl,"Two-gamma inv. mass for mod %d, cell (%d,%d)",iMod,iX,iZ);
	fHmpi0[iMod][iX][iZ] = new TH1F(hname,htitl,100,0.,300.);
	fOutputContainer->Add(fHmpi0[iMod][iX][iZ]);
      }
    }
  }

  fHmgg = new TH1F("hmgg","2-cluster invariant mass",100,0.,300.);
  fOutputContainer->Add(fHmgg);
	
  fCalibData = new AliEMCALCalibData();
	
  if(!fCalibData)
    AliFatal("Calibration parameters not found in CDB!");
  
  printf("Get Geometry : %s\n", fEMCALGeoName.Data());
  fEMCALGeo =  AliEMCALGeometry::GetInstance(fEMCALGeoName) ;	

}

//__________________________________________________
void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
{
  //Analysis per event.
  if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection <<< Event %d >>>\n",(Int_t)Entry());
  
  AliAODEvent* aod = 0x0;
  if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) {
    //Input are ESDs
    aod = dynamic_cast<AliAODEvent*>(InputEvent());
    // Create new AOD with only CaloClusters and CaloCells
    if(fCopyAOD) CreateAODFromAOD();
  }
  else  if(!strcmp(InputEvent()->GetName(),"AliESDEvent")) {
    //Input are ESDs
    aod = AODEvent();
    // Create AOD with CaloClusters and use it as input.
    // If filtering task is already executed, this is not needed.
    if(fCopyAOD) CreateAODFromESD();
  }
  else {
    printf("AliAnalysisTaskEMCALPi0CalibSelection: Unknown event type, STOP!\n");
    abort();
  }	
  
  Double_t v[] = {aod->GetVertex(0)->GetX(),aod->GetVertex(0)->GetY(),aod->GetVertex(0)->GetZ()}; //to check!!
  //aod->GetVertex()->GetXYZ(v) ;
  //TVector3 vtx(v); 
  
  //if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",vtx.X(),vtx.Y(),vtx.Z());
  if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",v[0],v[1],v[2]);
  
  Int_t runNum = aod->GetRunNumber();
  if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
  
  //Get the matrix with geometry information
  //Still not implemented in AOD, just a workaround to be able to work at least with ESDs	
  if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) {
    if(DebugLevel() > 1) 
      printf("AliAnalysisTaskEMCALPi0CalibSelection Use ideal geometry, values geometry matrix not kept in AODs.\n");
  }
  else{	
    if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Load Misaligned matrices. \n");
    AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
    for(Int_t mod=0; mod < 12; mod++){ 
      if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
    }
  }
  
  if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
  
  //AliEMCALPID pid;
  
  Int_t maxId;
  Int_t iSupMod = -1;
  Int_t iTower  = -1;
  Int_t iIphi   = -1;
  Int_t iIeta   = -1;
  Int_t iphi   = -1;
  Int_t ieta   = -1;
  
  TLorentzVector p1;
  TLorentzVector p2;
  TLorentzVector p12;
  
  
  TRefArray * caloClustersArr  = new TRefArray();
  aod->GetEMCALClusters(caloClustersArr);
  
  const Int_t kNumberOfEMCALClusters   = caloClustersArr->GetEntries() ;
  if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection CaloClusters: %d\n", kNumberOfEMCALClusters);
  
  // EMCAL cells
  AliAODCaloCells *emCells = aod->GetEMCALCells();
	
  // Check if is old data not modified for calibration and change the fraction factor and calibrate amplitudes
  // Only for aliroot older than release 17 tag 3
  if(fOldData && !fCopyAOD){
	  //Calibrate amplitudes
	  AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
	  for (Int_t iCell = 0; iCell < aodEMcells.GetNumberOfCells(); iCell++) {      
			aodEMcells.SetCell(iCell,aodEMcells.GetCellNumber(iCell),aodEMcells.GetAmplitude(iCell)*0.0153);
	  }
	  
	  //Change fraction
	  for(Int_t iClu=0; iClu<kNumberOfEMCALClusters; iClu++) {
		  AliAODCaloCluster *cl = (AliAODCaloCluster *) caloClustersArr->At(iClu);
		  if(!cl->IsEMCALCluster()) continue; // EMCAL cluster!
		  Double_t *fraction = cl->GetCellsAmplitudeFraction();
		  for(Int_t i = 0; i < cl->GetNCells() ; i++) fraction[i] = 1;
		  cl->SetCellsAmplitudeFraction(fraction);	  
	  }
  }
	
  // loop over EMCAL clusters
  for(Int_t iClu=0; iClu<kNumberOfEMCALClusters; iClu++) {
    
    AliAODCaloCluster *c1 = (AliAODCaloCluster *) caloClustersArr->At(iClu);
    if(!c1->IsEMCALCluster()) continue; // EMCAL cluster!
		  
	Float_t e1i = c1->E();   // cluster energy before correction   
	if(e1i<fEmin) continue;

	if(DebugLevel() > 2){ 
		printf("Std  : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClu,e1i, c1->GetDispersion(),c1->GetM02(),c1->GetM20());
		Double_t pos[]={0,0,0};
		c1->GetPosition(pos);
		printf("Std  : i %d, x %f, y %f, z %f\n",iClu, pos[0], pos[1], pos[2]);
	}
	  
    AliEMCALAodCluster clu1(*c1);
    clu1.Recalibrate(fCalibData, emCells, fEMCALGeoName);
	  printf("recalibrated, now calculate the rest\n");
    clu1.EvalAll(fLogWeight, fEMCALGeoName);
    //clu1.EnergyCorrection(&pid) ;
	  
	if(DebugLevel() > 2){ 
		printf("Recal: i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClu,clu1.E(), clu1.GetDispersion(),clu1.GetM02(),clu1.GetM20());    
		Double_t pos2[]={0,0,0};
		clu1.GetPosition(pos2);
		printf("Recal: i %d, x %f, y %f, z %f\n",iClu, pos2[0], pos2[1], pos2[2]);
	}
	  
	clu1.GetMomentum(p1,v);

	MaxEnergyCellPos(emCells,&clu1,maxId);
    
    //Get from the absid the supermodule, tower and eta/phi numbers
    fEMCALGeo->GetCellIndex(maxId,iSupMod,iTower,iIphi,iIeta); 
    //Gives SuperModule and Tower numbers
    fEMCALGeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
					   iIphi, iIeta,iphi,ieta);    
    
    Float_t e1ii = clu1.E(); // cluster energy after correction
    
    for (Int_t jClu=iClu; jClu<kNumberOfEMCALClusters; jClu++) {
      AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
      if(!c2->IsEMCALCluster())   continue; // EMCAL cluster!
      if(c2->IsEqual(c1)) continue;
      
      Float_t e2i = c2->E();
      if(e2i<fEmin) continue;
      
      AliEMCALAodCluster clu2(*c2);
      //printf("i2 %d, E %f\n",iClu,e2i);
      clu2.Recalibrate(fCalibData, emCells,fEMCALGeoName);
      clu2.EvalAll(fLogWeight,fEMCALGeoName);
      //clu2.EnergyCorrection(&pid) ;
     // printf("i2 %d, Erc %f\n",iClu,clu2.E());

      clu2.GetMomentum(p2,v);
      Float_t e2ii = clu2.E();
      
      p12 = p1+p2;
      Float_t invmass = p12.M()*1000; 

      fHmgg->Fill(invmass); 

      //printf("iSM %d, ieta %d, iphi %d,  mass %f  \n",iSupMod, ieta, iphi, invmass);
      if(invmass < 300)
	  fHmpi0[iSupMod][ieta][iphi]->Fill(invmass);
      
      
      if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Mass in (SM%d,%d,%d): %.3f GeV  E1_i=%f E1_ii=%f  E2_i=%f E2_ii=%f\n",
				  iSupMod,iphi,ieta,p12.M(),e1i,e1ii,e2i,e2ii);
    }
    
  } // end of loop over EMCAL clusters
  
  delete caloClustersArr;
  PostData(1,fOutputContainer);
}

//__________________________________________________
void AliAnalysisTaskEMCALPi0CalibSelection::MaxEnergyCellPos(AliAODCaloCells* const cells, AliAODCaloCluster* const clu, Int_t& maxId)
{
  //For a given CaloCluster calculates the absId of the cell 
  //with maximum energy deposit.
  
  Double_t eMax = -111;
  
  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
    Int_t cellAbsId = clu->GetCellAbsId(iDig);
    Double_t eCell  = cells->GetCellAmplitude(cellAbsId)*clu->GetCellAmplitudeFraction(iDig);
    if(eCell>eMax)  { 
      eMax = eCell; 
      maxId = cellAbsId;
    }
  }

}

//__________________________________________________
void AliAnalysisTaskEMCALPi0CalibSelection::SetCalibCorrections(AliEMCALCalibData* const cdata)
{
  //Set new correction factors (~1) to calibration coefficients, delete previous.

   if(fCalibData) delete fCalibData;
   fCalibData = cdata;
	
}
Commit	Line	Data
375cec9b	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: Boris Polishchuk *
	5	* Adapted to AOD reading by Gustavo Conesa *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//---------------------------------------------------------------------------//
	17	// //
	18	// Fill histograms (one per cell) with two-cluster invariant mass //
	19	// using calibration coefficients of the previous iteration. //
	20	// Histogram for a given cell is filled if the most energy of one cluster //
	21	// is deposited in this cell and the other cluster could be anywherein EMCAL.//
	22	// //
	23	//---------------------------------------------------------------------------//
	24
	25	//#include <cstdlib>
	26	//#include <Riostream.h>
	27	// Root
	28	#include "TLorentzVector.h"
	29	//#include "TVector3.h"
	30	#include "TRefArray.h"
	31	#include "TList.h"
	32	#include "TH1F.h"
	33
	34	// AliRoot
	35	#include "AliAnalysisTaskEMCALPi0CalibSelection.h"
	36	#include "AliAODEvent.h"
	37	#include "AliESDEvent.h"
	38	#include "AliESDCaloCluster.h"
	39	#include "AliESDCaloCells.h"
	40	#include "AliEMCALGeometry.h"
	41	#include "AliAODCaloCluster.h"
	42	#include "AliAODCaloCells.h"
	43	#include "AliEMCALAodCluster.h"
	44	#include "AliEMCALCalibData.h"
	45
	46	ClassImp(AliAnalysisTaskEMCALPi0CalibSelection)
	47
	48	//__________________________________________________
	49	AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection() :
	50	AliAnalysisTaskSE(),fEMCALGeo(0x0),fCalibData(0x0), fEmin(0.), fLogWeight(4.5), fCopyAOD(kFALSE),
	51	fEMCALGeoName("EMCAL_COMPLETE"), fOldData(kFALSE),fOutputContainer(0x0),fHmgg(0x0)
	52	{
	53	//Default constructor.
	54
bdd2a262	55	for(Int_t iMod=0; iMod < 12; iMod++) {
	56	for(Int_t iX=0; iX<24; iX++) {
	57	for(Int_t iZ=0; iZ<48; iZ++) {
375cec9b	58	fHmpi0[iMod][iX][iZ]=0;
	59	}
	60	}
	61	}
	62
	63	}
	64
	65	//__________________________________________________
	66	AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(const char* name) :
	67	AliAnalysisTaskSE(name),fEMCALGeo(0x0),fCalibData(0x0), fEmin(0.), fLogWeight(4.5), fCopyAOD(kFALSE),
	68	fEMCALGeoName("EMCAL_COMPLETE"), fOldData(kFALSE), fOutputContainer(0x0),fHmgg(0x0)
	69	{
	70	//Named constructor which should be used.
	71
	72	DefineOutput(1,TList::Class());
	73
bdd2a262	74	for(Int_t iMod=0; iMod < 12; iMod++) {
	75	for(Int_t iX=0; iX<24; iX++) {
	76	for(Int_t iZ=0; iZ<48; iZ++) {
375cec9b	77	fHmpi0[iMod][iX][iZ]=0;
	78	}
	79	}
	80	}
	81
	82	}
	83
	84	//__________________________________________________
	85	AliAnalysisTaskEMCALPi0CalibSelection::~AliAnalysisTaskEMCALPi0CalibSelection()
	86	{
	87	//Destructor.
	88
	89	if(fOutputContainer){
	90	fOutputContainer->Delete() ;
	91	delete fOutputContainer ;
	92	}
	93
	94	if(fCalibData) delete fCalibData;
	95	if(fEMCALGeo) delete fEMCALGeo;
	96
	97	}
	98
	99
	100	//__________________________________________________
	101	void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromAOD()
	102	{
	103	// Copy AOD header, vertex, CaloClusters and CaloCells to output AOD
	104
	105	AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
	106
	107	// set arrays and pointers
	108	Float_t posF[3];
	109	Double_t pos[3];
	110
	111	Double_t covVtx[6];
	112
	113	for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
	114
	115	// Update the header
	116	AliAODHeader*headerin = aod->GetHeader();
	117	AliAODHeader* header = AODEvent()->GetHeader();
	118	header->SetRunNumber(headerin->GetRunNumber());
	119	header->SetBunchCrossNumber(headerin->GetBunchCrossNumber());
	120	header->SetOrbitNumber(headerin->GetOrbitNumber());
	121	header->SetPeriodNumber(headerin->GetPeriodNumber());
	122	header->SetEventType(headerin->GetEventType());
	123	header->SetMuonMagFieldScale(headerin->GetMuonMagFieldScale());
	124	header->SetCentrality(headerin->GetCentrality());
	125
	126	header->SetTriggerMask(headerin->GetTriggerMask());
	127	header->SetTriggerCluster(headerin->GetTriggerCluster());
	128	header->SetMagneticField(headerin->GetMagneticField());
	129	header->SetZDCN1Energy(headerin->GetZDCN1Energy());
	130	header->SetZDCP1Energy(headerin->GetZDCP1Energy());
	131	header->SetZDCN2Energy(headerin->GetZDCN2Energy());
	132	header->SetZDCP2Energy(headerin->GetZDCP2Energy());
	133	header->SetZDCEMEnergy(headerin->GetZDCEMEnergy(0),headerin->GetZDCEMEnergy(1));
	134	Float_t diamxy[2]={aod->GetDiamondX(),aod->GetDiamondY()};
	135	Float_t diamcov[3]; aod->GetDiamondCovXY(diamcov);
	136	header->SetDiamond(diamxy,diamcov);
	137	//
	138	//
	139	Int_t nVertices = 1 ;/* = prim. vtx*/;
	140	Int_t nCaloClus = aod->GetNCaloClusters();
141
142	AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
143
144	// Access to the AOD container of vertices
145	TClonesArray &vertices = *(AODEvent()->GetVertices());
146	Int_t jVertices=0;
147
148	// Add primary vertex. The primary tracks will be defined
149	// after the loops on the composite objects (V0, cascades, kinks)
150	const AliAODVertex *vtx = aod->GetPrimaryVertex();
151
152	vtx->GetXYZ(pos); // position
153	vtx->GetCovMatrix(covVtx); //covariance matrix
154
155	AliAODVertex * primary = new(vertices[jVertices++])
156	AliAODVertex(pos, covVtx, vtx->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
157	primary->SetName(vtx->GetName());
158	primary->SetTitle(vtx->GetTitle());
159
160	// Access to the AOD container of clusters
161	TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
162	Int_t jClusters=0;
163
164	for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
165
166	AliAODCaloCluster * cluster = aod->GetCaloCluster(iClust);
167
168	//Check if it is a EMCAL cluster
169	if(!cluster->IsEMCALCluster()) continue ;
170
171	Int_t id = cluster->GetID();
172	Float_t energy = cluster->E();
173	cluster->GetPosition(posF);
174	Char_t ttype = cluster->GetType();
175
176	AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
177	AliAODCaloCluster(id,
178	0,
179	0x0,
180	energy,
181	posF,
182	NULL,
183	ttype);
184
185	caloCluster->SetCaloCluster(cluster->GetDistToBadChannel(),
186	cluster->GetDispersion(),
187	cluster->GetM20(), cluster->GetM02(),
188	cluster->GetEmcCpvDistance(),
189	cluster->GetNExMax(),cluster->GetTOF()) ;
190
191	caloCluster->SetPIDFromESD(cluster->PID());
192	caloCluster->SetNCells(cluster->GetNCells());
193	caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
194	if(!fOldData){
195	caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
196	}
197	else{
198	Double_t *fraction = cluster->GetCellsAmplitudeFraction();
199	for(Int_t i = 0; i < cluster->GetNCells() ; i++) fraction[i] = 1;
200	caloCluster->SetCellsAmplitudeFraction(fraction);
201	}
202
203	}
204
205	caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
206	// end of loop on calo clusters
207
208	// fill EMCAL cell info
209	if (aod->GetEMCALCells()) { // protection against missing AOD information
210	AliAODCaloCells &aodinEMcells = *(aod->GetEMCALCells());
211	Int_t nEMcell = aodinEMcells.GetNumberOfCells() ;
212
213	AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
214	aodEMcells.CreateContainer(nEMcell);
215	aodEMcells.SetType(AliAODCaloCells::kEMCAL);
216
217	Double_t calibFactor = 1;
218	if(fOldData) calibFactor = 0.0153;
219
220	for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
221	aodEMcells.SetCell(iCell,aodinEMcells.GetCellNumber(iCell),aodinEMcells.GetAmplitude(iCell)*calibFactor);
222	}
223	aodEMcells.Sort();
224
225	}
226
227	}
228
229	//__________________________________________________
230	void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromESD()
231	{
232
233	// Copy Header, Vertex, CaloClusters and CaloCells from ESDs to AODs
234
235	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
236
237	// set arrays and pointers
238	Float_t posF[3];
239	Double_t pos[3];
240
241	Double_t covVtx[6];
242
243	for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
244
245	// Update the header
246
247	AliAODHeader* header = AODEvent()->GetHeader();
248	header->SetRunNumber(esd->GetRunNumber());
249	header->SetBunchCrossNumber(esd->GetBunchCrossNumber());
250	header->SetOrbitNumber(esd->GetOrbitNumber());
251	header->SetPeriodNumber(esd->GetPeriodNumber());
252	header->SetEventType(esd->GetEventType());
253	header->SetMuonMagFieldScale(-999.); // FIXME
254	header->SetCentrality(-999.); // FIXME
255
256
257	header->SetTriggerMask(esd->GetTriggerMask());
258	header->SetTriggerCluster(esd->GetTriggerCluster());
259	header->SetMagneticField(esd->GetMagneticField());
260	header->SetZDCN1Energy(esd->GetZDCN1Energy());
261	header->SetZDCP1Energy(esd->GetZDCP1Energy());
262	header->SetZDCN2Energy(esd->GetZDCN2Energy());
263	header->SetZDCP2Energy(esd->GetZDCP2Energy());
264	header->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
265	Float_t diamxy[2]={esd->GetDiamondX(),esd->GetDiamondY()};
266	Float_t diamcov[3]; esd->GetDiamondCovXY(diamcov);
267	header->SetDiamond(diamxy,diamcov);
268	//
269	//
270	Int_t nVertices = 1 ;/* = prim. vtx*/;
271	Int_t nCaloClus = esd->GetNumberOfCaloClusters();
272
273	AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
274
275	// Access to the AOD container of vertices
276	TClonesArray &vertices = *(AODEvent()->GetVertices());
277	Int_t jVertices=0;
278
279	// Add primary vertex. The primary tracks will be defined
280	// after the loops on the composite objects (V0, cascades, kinks)
281	const AliESDVertex *vtx = esd->GetPrimaryVertex();
282
283	vtx->GetXYZ(pos); // position
284	vtx->GetCovMatrix(covVtx); //covariance matrix
285
286	AliAODVertex * primary = new(vertices[jVertices++])
287	AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
288	primary->SetName(vtx->GetName());
289	primary->SetTitle(vtx->GetTitle());
290
291	// Access to the AOD container of clusters
292	TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
293	Int_t jClusters=0;
294
295	for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
296
297	AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
298
299	//Check which calorimeter information we want to keep.
300	if(!cluster->IsEMCAL()) continue ;
301
302	Int_t id = cluster->GetID();
303	Float_t energy = cluster->E();
304	cluster->GetPosition(posF);
305
306	AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
307	AliAODCaloCluster(id,
308	0,
309	0x0,
310	energy,
311	posF,
312	NULL,
313	AliAODCluster::kEMCALClusterv1);
314
315	caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
316	cluster->GetClusterDisp(),
317	cluster->GetM20(), cluster->GetM02(),
318	cluster->GetEmcCpvDistance(),
319	cluster->GetNExMax(),cluster->GetTOF()) ;
320
321	caloCluster->SetPIDFromESD(cluster->GetPid());
322	caloCluster->SetNCells(cluster->GetNCells());
323	caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
324	if(!fOldData){
325	caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
326	}
327	else{
328	Double_t *fraction = cluster->GetCellsAmplitudeFraction();
329	for(Int_t i = 0; i < cluster->GetNCells() ; i++) fraction[i] = 1;
330	caloCluster->SetCellsAmplitudeFraction(fraction);
331	}
332	}
333	caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
334	// end of loop on calo clusters
335
336	// fill EMCAL cell info
337
338	if( esd->GetEMCALCells()) { // protection against missing ESD information
339	AliESDCaloCells &esdEMcells = *(esd->GetEMCALCells());
340	Int_t nEMcell = esdEMcells.GetNumberOfCells() ;
341
342	AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
343	aodEMcells.CreateContainer(nEMcell);
344	aodEMcells.SetType(AliAODCaloCells::kEMCAL);
345
346	Double_t calibFactor = 1;
347	if(fOldData) calibFactor = 0.0153;
348
349	for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
350	aodEMcells.SetCell(iCell,esdEMcells.GetCellNumber(iCell),esdEMcells.GetAmplitude(iCell)*calibFactor);
351	}
352	aodEMcells.Sort();
353
354	}
355
356	}
357
358	//__________________________________________________
359	void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
360	{
361	//Create output container
362	fOutputContainer = new TList();
363
364	char hname[128], htitl[128];
365
366	for(Int_t iMod=0; iMod < AliEMCALGeoParams::fgkEMCALModules; iMod++) {
367	for(Int_t iX=0; iX<AliEMCALGeoParams::fgkEMCALCols; iX++) {
368	for(Int_t iZ=0; iZ<AliEMCALGeoParams::fgkEMCALRows; iZ++) {
369	sprintf(hname,"%d_%d_%d",iMod,iX,iZ);
370	sprintf(htitl,"Two-gamma inv. mass for mod %d, cell (%d,%d)",iMod,iX,iZ);
371	fHmpi0[iMod][iX][iZ] = new TH1F(hname,htitl,100,0.,300.);
372	fOutputContainer->Add(fHmpi0[iMod][iX][iZ]);
373	}
374	}
375	}
376
377	fHmgg = new TH1F("hmgg","2-cluster invariant mass",100,0.,300.);
378	fOutputContainer->Add(fHmgg);
379
380	fCalibData = new AliEMCALCalibData();
381
382	if(!fCalibData)
383	AliFatal("Calibration parameters not found in CDB!");
384
385	printf("Get Geometry : %s\n", fEMCALGeoName.Data());
386	fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
387
388	}
389
390	//__________________________________________________
391	void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
392	{
393	//Analysis per event.
394	if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection <<< Event %d >>>\n",(Int_t)Entry());
395
396	AliAODEvent* aod = 0x0;
397	if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) {
398	//Input are ESDs
399	aod = dynamic_cast<AliAODEvent*>(InputEvent());
400	// Create new AOD with only CaloClusters and CaloCells
401	if(fCopyAOD) CreateAODFromAOD();
402	}
403	else if(!strcmp(InputEvent()->GetName(),"AliESDEvent")) {
404	//Input are ESDs
405	aod = AODEvent();
406	// Create AOD with CaloClusters and use it as input.
407	// If filtering task is already executed, this is not needed.
408	if(fCopyAOD) CreateAODFromESD();
409	}
410	else {
411	printf("AliAnalysisTaskEMCALPi0CalibSelection: Unknown event type, STOP!\n");
412	abort();
413	}
414
415	Double_t v[] = {aod->GetVertex(0)->GetX(),aod->GetVertex(0)->GetY(),aod->GetVertex(0)->GetZ()}; //to check!!
416	//aod->GetVertex()->GetXYZ(v) ;
417	//TVector3 vtx(v);
418
419	//if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",vtx.X(),vtx.Y(),vtx.Z());
420	if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",v[0],v[1],v[2]);
421
422	Int_t runNum = aod->GetRunNumber();
423	if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
424
425	//Get the matrix with geometry information
426	//Still not implemented in AOD, just a workaround to be able to work at least with ESDs
427	if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) {
428	if(DebugLevel() > 1)
429	printf("AliAnalysisTaskEMCALPi0CalibSelection Use ideal geometry, values geometry matrix not kept in AODs.\n");
430	}
431	else{
432	if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Load Misaligned matrices. \n");
433	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
bdd2a262	434	for(Int_t mod=0; mod < 12; mod++){
375cec9b	435	if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
	436	}
	437	}
	438
	439	if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
	440
	441	//AliEMCALPID pid;
	442
	443	Int_t maxId;
	444	Int_t iSupMod = -1;
	445	Int_t iTower = -1;
	446	Int_t iIphi = -1;
	447	Int_t iIeta = -1;
	448	Int_t iphi = -1;
	449	Int_t ieta = -1;
	450
	451	TLorentzVector p1;
	452	TLorentzVector p2;
	453	TLorentzVector p12;
	454
	455
	456	TRefArray * caloClustersArr = new TRefArray();
	457	aod->GetEMCALClusters(caloClustersArr);
	458
	459	const Int_t kNumberOfEMCALClusters = caloClustersArr->GetEntries() ;
	460	if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection CaloClusters: %d\n", kNumberOfEMCALClusters);
	461
	462	// EMCAL cells
	463	AliAODCaloCells *emCells = aod->GetEMCALCells();
	464
	465	// Check if is old data not modified for calibration and change the fraction factor and calibrate amplitudes
	466	// Only for aliroot older than release 17 tag 3
	467	if(fOldData && !fCopyAOD){
	468	//Calibrate amplitudes
	469	AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
	470	for (Int_t iCell = 0; iCell < aodEMcells.GetNumberOfCells(); iCell++) {
	471	aodEMcells.SetCell(iCell,aodEMcells.GetCellNumber(iCell),aodEMcells.GetAmplitude(iCell)*0.0153);
	472	}
	473
	474	//Change fraction
	475	for(Int_t iClu=0; iClu<kNumberOfEMCALClusters; iClu++) {
	476	AliAODCaloCluster cl = (AliAODCaloCluster ) caloClustersArr->At(iClu);
	477	if(!cl->IsEMCALCluster()) continue; // EMCAL cluster!
	478	Double_t *fraction = cl->GetCellsAmplitudeFraction();
	479	for(Int_t i = 0; i < cl->GetNCells() ; i++) fraction[i] = 1;
	480	cl->SetCellsAmplitudeFraction(fraction);
	481	}
	482	}
	483
	484	// loop over EMCAL clusters
	485	for(Int_t iClu=0; iClu<kNumberOfEMCALClusters; iClu++) {
	486
	487	AliAODCaloCluster c1 = (AliAODCaloCluster ) caloClustersArr->At(iClu);
	488	if(!c1->IsEMCALCluster()) continue; // EMCAL cluster!
	489
	490	Float_t e1i = c1->E(); // cluster energy before correction
	491	if(e1i<fEmin) continue;
	492
	493	if(DebugLevel() > 2){
	494	printf("Std : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClu,e1i, c1->GetDispersion(),c1->GetM02(),c1->GetM20());
	495	Double_t pos[]={0,0,0};
	496	c1->GetPosition(pos);
	497	printf("Std : i %d, x %f, y %f, z %f\n",iClu, pos[0], pos[1], pos[2]);
	498	}
499
500	AliEMCALAodCluster clu1(*c1);
501	clu1.Recalibrate(fCalibData, emCells, fEMCALGeoName);
502	printf("recalibrated, now calculate the rest\n");
503	clu1.EvalAll(fLogWeight, fEMCALGeoName);
504	//clu1.EnergyCorrection(&pid) ;
505
506	if(DebugLevel() > 2){
507	printf("Recal: i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClu,clu1.E(), clu1.GetDispersion(),clu1.GetM02(),clu1.GetM20());
508	Double_t pos2[]={0,0,0};
509	clu1.GetPosition(pos2);
510	printf("Recal: i %d, x %f, y %f, z %f\n",iClu, pos2[0], pos2[1], pos2[2]);
511	}
512
513	clu1.GetMomentum(p1,v);
514
515	MaxEnergyCellPos(emCells,&clu1,maxId);
516
517	//Get from the absid the supermodule, tower and eta/phi numbers
518	fEMCALGeo->GetCellIndex(maxId,iSupMod,iTower,iIphi,iIeta);
519	//Gives SuperModule and Tower numbers
520	fEMCALGeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
521	iIphi, iIeta,iphi,ieta);
522
523	Float_t e1ii = clu1.E(); // cluster energy after correction
524
525	for (Int_t jClu=iClu; jClu<kNumberOfEMCALClusters; jClu++) {
526	AliAODCaloCluster c2 = (AliAODCaloCluster ) caloClustersArr->At(jClu);
527	if(!c2->IsEMCALCluster()) continue; // EMCAL cluster!
528	if(c2->IsEqual(c1)) continue;
529
530	Float_t e2i = c2->E();
531	if(e2i<fEmin) continue;
532
533	AliEMCALAodCluster clu2(*c2);
534	//printf("i2 %d, E %f\n",iClu,e2i);
535	clu2.Recalibrate(fCalibData, emCells,fEMCALGeoName);
536	clu2.EvalAll(fLogWeight,fEMCALGeoName);
537	//clu2.EnergyCorrection(&pid) ;
538	// printf("i2 %d, Erc %f\n",iClu,clu2.E());
539
540	clu2.GetMomentum(p2,v);
541	Float_t e2ii = clu2.E();
542
543	p12 = p1+p2;
544	Float_t invmass = p12.M()*1000;
545
546	fHmgg->Fill(invmass);
547
548	//printf("iSM %d, ieta %d, iphi %d, mass %f \n",iSupMod, ieta, iphi, invmass);
549	if(invmass < 300)
550	fHmpi0[iSupMod][ieta][iphi]->Fill(invmass);
551
552
553	if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Mass in (SM%d,%d,%d): %.3f GeV E1_i=%f E1_ii=%f E2_i=%f E2_ii=%f\n",
554	iSupMod,iphi,ieta,p12.M(),e1i,e1ii,e2i,e2ii);
555	}
556
557	} // end of loop over EMCAL clusters
558
559	delete caloClustersArr;
560	PostData(1,fOutputContainer);
561	}
562
563	//__________________________________________________
bdd2a262	564	void AliAnalysisTaskEMCALPi0CalibSelection::MaxEnergyCellPos(AliAODCaloCells* const cells, AliAODCaloCluster* const clu, Int_t& maxId)
375cec9b	565	{
	566	//For a given CaloCluster calculates the absId of the cell
	567	//with maximum energy deposit.
	568
	569	Double_t eMax = -111;
	570
	571	for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
	572	Int_t cellAbsId = clu->GetCellAbsId(iDig);
	573	Double_t eCell = cells->GetCellAmplitude(cellAbsId)*clu->GetCellAmplitudeFraction(iDig);
	574	if(eCell>eMax) {
	575	eMax = eCell;
	576	maxId = cellAbsId;
	577	}
	578	}
	579
	580	}
	581
	582	//__________________________________________________
bdd2a262	583	void AliAnalysisTaskEMCALPi0CalibSelection::SetCalibCorrections(AliEMCALCalibData* const cdata)
375cec9b	584	{
	585	//Set new correction factors (~1) to calibration coefficients, delete previous.
	586
	587	if(fCalibData) delete fCalibData;
	588	fCalibData = cdata;
	589
	590	}
	591