[u/mrichter/AliRoot.git] / PWG4 / CaloCalib / AliAnalysisTaskPHOSPi0CalibSelection.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 anywhere in PHOS.//
//                                                                           //
//---------------------------------------------------------------------------//

#include <cstdlib>

// Root 
#include "TLorentzVector.h"
#include "TVector3.h"
#include "TRefArray.h"
#include "TList.h"

// AliRoot
#include "AliAnalysisTaskPHOSPi0CalibSelection.h"
#include "AliAODEvent.h"
#include "AliESDEvent.h"
#include "AliESDCaloCluster.h"
#include "AliESDCaloCells.h"
#include "AliPHOSAodCluster.h"
#include "AliPHOSPIDv1.h"

ClassImp(AliAnalysisTaskPHOSPi0CalibSelection)

//__________________________________________________
AliAnalysisTaskPHOSPi0CalibSelection::AliAnalysisTaskPHOSPi0CalibSelection() :
AliAnalysisTaskSE(),fOutputContainer(0x0),fPhosGeo(0x0),fCalibData(0x0),fHmgg(0x0),
fEmin(0.), fLogWeight(4.5), fCopyAOD(kFALSE)
{
  //Default constructor.

  for(Int_t iMod=0; iMod<5; iMod++) {
    for(Int_t iX=0; iX<64; iX++) {
      for(Int_t iZ=0; iZ<56; iZ++) {
	fHmpi0[iMod][iX][iZ]=0;
      }
    } 
  }
		
}

//__________________________________________________
AliAnalysisTaskPHOSPi0CalibSelection::AliAnalysisTaskPHOSPi0CalibSelection(const char* name) :
  AliAnalysisTaskSE(name),fOutputContainer(0x0),fPhosGeo(0x0),fCalibData(0x0),fHmgg(0x0),
  fEmin(0.), fLogWeight(4.5), fCopyAOD(kFALSE)
{
  //Named constructor which should be used.
  
  DefineOutput(1,TList::Class());
  
  for(Int_t iMod=0; iMod<5; iMod++) {
    for(Int_t iX=0; iX<64; iX++) {
      for(Int_t iZ=0; iZ<56; iZ++) {
	fHmpi0[iMod][iX][iZ]=0;
      }
    } 
  }

}

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

//__________________________________________________
void AliAnalysisTaskPHOSPi0CalibSelection::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->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 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 PHOS cluster
    if(!cluster->IsPHOS())  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->GetDistanceToBadChannel(),
				cluster->GetDispersion(),
				cluster->GetM20(), cluster->GetM02(),
				cluster->GetEmcCpvDistance(),  
				cluster->GetNExMax(),cluster->GetTOF()) ;
    
    caloCluster->SetPIDFromESD(cluster->GetPID());
    caloCluster->SetNCells(cluster->GetNCells());
    caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
    caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
    
  } 

  caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters	 
  // end of loop on calo clusters
  
  // fill PHOS cell info
  if (aod->GetPHOSCells()) { // protection against missing AOD information
    AliAODCaloCells &aodinPHcells = *(aod->GetPHOSCells());
    Int_t nPHcell = aodinPHcells.GetNumberOfCells() ;
    
    AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
    aodPHcells.CreateContainer(nPHcell);
    aodPHcells.SetType(AliVCaloCells::kPHOSCell);
    for (Int_t iCell = 0; iCell < nPHcell; iCell++) {      
      aodPHcells.SetCell(iCell,aodinPHcells.GetCellNumber(iCell),aodinPHcells.GetAmplitude(iCell));
    }
    aodPHcells.Sort();
  }
  
  return;
}

//__________________________________________________
void AliAnalysisTaskPHOSPi0CalibSelection::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->IsPHOS())  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::kPHOSNeutral);
    
    caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
				cluster->GetDispersion(),
				cluster->GetM20(), cluster->GetM02(),
				cluster->GetEmcCpvDistance(),  
				cluster->GetNExMax(),cluster->GetTOF()) ;
    
    caloCluster->SetPIDFromESD(cluster->GetPID());
    caloCluster->SetNCells(cluster->GetNCells());
    caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
    caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
    
  } 
  caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters	 
  // end of loop on calo clusters
  
  // fill PHOS cell info
  if( esd->GetPHOSCells()) { // protection against missing ESD information
    AliESDCaloCells &esdPHcells = *(esd->GetPHOSCells());
    Int_t nPHcell = esdPHcells.GetNumberOfCells() ;
    
    AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
    aodPHcells.CreateContainer(nPHcell);
    aodPHcells.SetType(AliVCaloCells::kPHOSCell);
    for (Int_t iCell = 0; iCell < nPHcell; iCell++) {      
      aodPHcells.SetCell(iCell,esdPHcells.GetCellNumber(iCell),esdPHcells.GetAmplitude(iCell));
    }
    aodPHcells.Sort();
  }

}

//__________________________________________________
void AliAnalysisTaskPHOSPi0CalibSelection::UserCreateOutputObjects()
{
  //Create output container
  fOutputContainer = new TList();
  
  char hname[128], htitl[128];
  
  for(Int_t iMod=0; iMod<5; iMod++) {
    for(Int_t iX=0; iX<64; iX++) {
      for(Int_t iZ=0; iZ<56; 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 AliPHOSCalibData();
  fPhosGeo =  AliPHOSGeometry::GetInstance("IHEP") ;	

}

//__________________________________________________
void AliAnalysisTaskPHOSPi0CalibSelection::UserExec(Option_t* /* option */)
{
  //Analysis per event.
  if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection <<< 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("AliAnalysisTaskPHOSPi0CalibSelection: 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); //Check
	
  if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",vtx.X(),vtx.Y(),vtx.Z());
 
  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("AliAnalysisTaskPHOSPi0CalibSelection Use ideal geometry, values geometry matrix not kept in AODs.\n");
  }
  else{	
	  if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection Load Misaligned matrices. \n");
	  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
	  for(Int_t mod=0; mod<5; mod++){ 
		if(esd->GetPHOSMatrix(mod))
			fPhosGeo->SetMisalMatrix(esd->GetPHOSMatrix(mod),mod) ;
		}
  }
	
  if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);

  AliPHOSPIDv1 pid;

  Int_t mod1 = -1;
  TLorentzVector p1;
  TLorentzVector p2;
  TLorentzVector p12;
  Int_t relid[4] ;
  Int_t maxId;

  TRefArray * caloClustersArr  = new TRefArray();
  aod->GetPHOSClusters(caloClustersArr);
  
  const Int_t kNumberOfPhosClusters   = caloClustersArr->GetEntries() ;
  if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection CaloClusters: %d\n", kNumberOfPhosClusters);

  // PHOS cells
  AliAODCaloCells *phsCells = aod->GetPHOSCells();
  
  // loop over PHOS clusters
  for(Int_t iClu=0; iClu<kNumberOfPhosClusters; iClu++) {
    
    AliAODCaloCluster *c1 = (AliAODCaloCluster *) caloClustersArr->At(iClu);
    if(!c1->IsPHOS()) continue; // EMCAL cluster!

    Float_t e1i = c1->E();   // cluster energy before correction
    if(e1i<fEmin) continue;

    AliPHOSAodCluster clu1(*c1);
    clu1.Recalibrate(fCalibData, phsCells);
    clu1.EvalAll(fLogWeight,vtx);
    clu1.EnergyCorrection(&pid) ;

    clu1.GetMomentum(p1,v);

    MaxEnergyCellPos(phsCells,&clu1,maxId);
    fPhosGeo->AbsToRelNumbering(maxId, relid);

    mod1 = relid[0]-1;        // module
    Int_t iX = relid[2]-1;    // cluster X-coord
    Int_t iZ = relid[3]-1;    // cluster Z-coord
	Float_t e1ii = clu1.E(); // cluster energy after correction
    
    for (Int_t jClu=iClu; jClu<kNumberOfPhosClusters; jClu++) {
      AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
      if(!c2->IsPHOS())   continue; // EMCAL cluster!
      if(c2->IsEqual(c1)) continue;

      Float_t e2i = c2->E();
      if(e2i<fEmin) continue;
      
      AliPHOSAodCluster clu2(*c2);
      clu2.Recalibrate(fCalibData, phsCells);
      clu2.EvalAll(fLogWeight,vtx);
      clu2.EnergyCorrection(&pid) ;
        
      clu2.GetMomentum(p2,v);
	  Float_t e2ii = clu2.E();

      p12 = p1+p2;

      fHmgg->Fill(p12.M()*1000); 
      fHmpi0[mod1][iX][iZ]->Fill(p12.M()*1000);
      
	  if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection Mass in (mod%d,%d,%d): %.3f GeV  E1_i=%f E1_ii=%f  E2_i=%f E2_ii=%f\n",
	     mod1,iX,iZ,p12.M(),e1i,e1ii,e2i,e2ii);
    }
    
  } // end of loop over PHOS clusters
  
  delete caloClustersArr;
  PostData(1,fOutputContainer);
}

//__________________________________________________
void AliAnalysisTaskPHOSPi0CalibSelection::MaxEnergyCellPos(AliAODCaloCells *cells, AliAODCaloCluster* 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 AliAnalysisTaskPHOSPi0CalibSelection::SetCalibCorrections(AliPHOSCalibData* cdata)
{
  //Set new correction factors (~1) to calibration coefficients, delete previous.

   if(fCalibData) delete fCalibData;
   fCalibData = cdata;
	
}
Commit	Line	Data
e727d2f7	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: Boris Polishchuk *
7a4cf423	5	* Adapted to AOD reading by Gustavo Conesa *
e727d2f7	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 anywhere in PHOS.//
	22	// //
	23	//---------------------------------------------------------------------------//
	24
	25	#include <cstdlib>
	26
	27	// Root
	28	#include "TLorentzVector.h"
	29	#include "TVector3.h"
e727d2f7	30	#include "TRefArray.h"
	31	#include "TList.h"
	32
	33	// AliRoot
	34	#include "AliAnalysisTaskPHOSPi0CalibSelection.h"
e727d2f7	35	#include "AliAODEvent.h"
4b4efabc	36	#include "AliESDEvent.h"
7a4cf423	37	#include "AliESDCaloCluster.h"
7a4cf423	38	#include "AliESDCaloCells.h"
e727d2f7	39	#include "AliPHOSAodCluster.h"
e727d2f7	40	#include "AliPHOSPIDv1.h"
e727d2f7	41
e727d2f7	42	ClassImp(AliAnalysisTaskPHOSPi0CalibSelection)
e727d2f7	43
7a4cf423	44	//__________________________________________________
e727d2f7	45	AliAnalysisTaskPHOSPi0CalibSelection::AliAnalysisTaskPHOSPi0CalibSelection() :
4b4efabc	46	AliAnalysisTaskSE(),fOutputContainer(0x0),fPhosGeo(0x0),fCalibData(0x0),fHmgg(0x0),
7a4cf423	47	fEmin(0.), fLogWeight(4.5), fCopyAOD(kFALSE)
e727d2f7	48	{
	49	//Default constructor.
	50
	51	for(Int_t iMod=0; iMod<5; iMod++) {
	52	for(Int_t iX=0; iX<64; iX++) {
	53	for(Int_t iZ=0; iZ<56; iZ++) {
	54	fHmpi0[iMod][iX][iZ]=0;
	55	}
	56	}
	57	}
	58
	59	}
	60
7a4cf423	61	//__________________________________________________
e727d2f7	62	AliAnalysisTaskPHOSPi0CalibSelection::AliAnalysisTaskPHOSPi0CalibSelection(const char* name) :
4b4efabc	63	AliAnalysisTaskSE(name),fOutputContainer(0x0),fPhosGeo(0x0),fCalibData(0x0),fHmgg(0x0),
7a4cf423	64	fEmin(0.), fLogWeight(4.5), fCopyAOD(kFALSE)
e727d2f7	65	{
	66	//Named constructor which should be used.
	67
	68	DefineOutput(1,TList::Class());
	69
	70	for(Int_t iMod=0; iMod<5; iMod++) {
	71	for(Int_t iX=0; iX<64; iX++) {
	72	for(Int_t iZ=0; iZ<56; iZ++) {
	73	fHmpi0[iMod][iX][iZ]=0;
	74	}
	75	}
	76	}
	77
	78	}
	79
7a4cf423	80	//__________________________________________________
e727d2f7	81	AliAnalysisTaskPHOSPi0CalibSelection::~AliAnalysisTaskPHOSPi0CalibSelection()
	82	{
	83	//Destructor.
	84
	85	if(fOutputContainer){
	86	fOutputContainer->Delete() ;
	87	delete fOutputContainer ;
	88	}
4b4efabc	89
	90	if(fCalibData) delete fCalibData;
	91	if(fPhosGeo) delete fPhosGeo;
	92
e727d2f7	93	}
e727d2f7	94
7a4cf423	95	//__________________________________________________
	96	void AliAnalysisTaskPHOSPi0CalibSelection::CreateAODFromAOD()
	97	{
	98	// Copy AOD header, vertex, CaloClusters and CaloCells to output AOD
	99
	100	AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
	101
	102	// set arrays and pointers
	103	Float_t posF[3];
	104	Double_t pos[3];
	105
	106	Double_t covVtx[6];
	107
	108	for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
	109
	110	// Update the header
29a41d05	111	AliAODHeader*headerin = aod->GetHeader();
7a4cf423	112	AliAODHeader* header = AODEvent()->GetHeader();
	113	header->SetRunNumber(headerin->GetRunNumber());
	114	header->SetBunchCrossNumber(headerin->GetBunchCrossNumber());
	115	header->SetOrbitNumber(headerin->GetOrbitNumber());
	116	header->SetPeriodNumber(headerin->GetPeriodNumber());
	117	header->SetEventType(headerin->GetEventType());
29a41d05	118	header->SetMuonMagFieldScale(headerin->GetMuonMagFieldScale());
29a41d05	119	header->SetCentrality(headerin->GetCentrality());
7a4cf423	120
	121	header->SetTriggerMask(headerin->GetTriggerMask());
	122	header->SetTriggerCluster(headerin->GetTriggerCluster());
	123	header->SetMagneticField(headerin->GetMagneticField());
	124	header->SetZDCN1Energy(headerin->GetZDCN1Energy());
	125	header->SetZDCP1Energy(headerin->GetZDCP1Energy());
	126	header->SetZDCN2Energy(headerin->GetZDCN2Energy());
	127	header->SetZDCP2Energy(headerin->GetZDCP2Energy());
	128	header->SetZDCEMEnergy(headerin->GetZDCEMEnergy(0),headerin->GetZDCEMEnergy(1));
	129	Float_t diamxy[2]={aod->GetDiamondX(),aod->GetDiamondY()};
	130	Float_t diamcov[3]; aod->GetDiamondCovXY(diamcov);
	131	header->SetDiamond(diamxy,diamcov);
	132	//
	133	//
	134	Int_t nVertices = 1 ;/* = prim. vtx*/;
c8fe2783	135	Int_t nCaloClus = aod->GetNumberOfCaloClusters();
7a4cf423	136
	137	AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
	138
	139	// Access to the AOD container of vertices
	140	TClonesArray &vertices = *(AODEvent()->GetVertices());
	141	Int_t jVertices=0;
	142
	143	// Add primary vertex. The primary tracks will be defined
	144	// after the loops on the composite objects (V0, cascades, kinks)
	145	const AliAODVertex *vtx = aod->GetPrimaryVertex();
	146
	147	vtx->GetXYZ(pos); // position
	148	vtx->GetCovMatrix(covVtx); //covariance matrix
	149
	150	AliAODVertex * primary = new(vertices[jVertices++])
	151	AliAODVertex(pos, covVtx, vtx->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
	152	primary->SetName(vtx->GetName());
	153	primary->SetTitle(vtx->GetTitle());
	154
	155	// Access to the AOD container of clusters
	156	TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
	157	Int_t jClusters=0;
	158
	159	for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
	160
	161	AliAODCaloCluster * cluster = aod->GetCaloCluster(iClust);
	162
	163	//Check if it is a PHOS cluster
c8fe2783	164	if(!cluster->IsPHOS()) continue ;
7a4cf423	165
	166	Int_t id = cluster->GetID();
	167	Float_t energy = cluster->E();
	168	cluster->GetPosition(posF);
	169	Char_t ttype = cluster->GetType();
	170
	171	AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
	172	AliAODCaloCluster(id,
	173	0,
	174	0x0,
	175	energy,
	176	posF,
	177	NULL,
	178	ttype);
	179
c8fe2783	180	caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
7a4cf423	181	cluster->GetDispersion(),
	182	cluster->GetM20(), cluster->GetM02(),
	183	cluster->GetEmcCpvDistance(),
	184	cluster->GetNExMax(),cluster->GetTOF()) ;
	185
c8fe2783	186	caloCluster->SetPIDFromESD(cluster->GetPID());
7a4cf423	187	caloCluster->SetNCells(cluster->GetNCells());
	188	caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
	189	caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
	190
	191	}
	192
	193	caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
	194	// end of loop on calo clusters
	195
	196	// fill PHOS cell info
	197	if (aod->GetPHOSCells()) { // protection against missing AOD information
	198	AliAODCaloCells &aodinPHcells = *(aod->GetPHOSCells());
	199	Int_t nPHcell = aodinPHcells.GetNumberOfCells() ;
	200
	201	AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
	202	aodPHcells.CreateContainer(nPHcell);
c8fe2783	203	aodPHcells.SetType(AliVCaloCells::kPHOSCell);
7a4cf423	204	for (Int_t iCell = 0; iCell < nPHcell; iCell++) {
	205	aodPHcells.SetCell(iCell,aodinPHcells.GetCellNumber(iCell),aodinPHcells.GetAmplitude(iCell));
	206	}
	207	aodPHcells.Sort();
	208	}
	209
	210	return;
	211	}
	212
	213	//__________________________________________________
	214	void AliAnalysisTaskPHOSPi0CalibSelection::CreateAODFromESD()
	215	{
	216
	217	// Copy Header, Vertex, CaloClusters and CaloCells from ESDs to AODs
	218
	219	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
	220
	221	// set arrays and pointers
	222	Float_t posF[3];
	223	Double_t pos[3];
	224
	225	Double_t covVtx[6];
	226
	227	for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
	228
	229	// Update the header
	230
	231	AliAODHeader* header = AODEvent()->GetHeader();
	232	header->SetRunNumber(esd->GetRunNumber());
	233	header->SetBunchCrossNumber(esd->GetBunchCrossNumber());
	234	header->SetOrbitNumber(esd->GetOrbitNumber());
	235	header->SetPeriodNumber(esd->GetPeriodNumber());
	236	header->SetEventType(esd->GetEventType());
	237	header->SetMuonMagFieldScale(-999.); // FIXME
	238	header->SetCentrality(-999.); // FIXME
	239
	240
	241	header->SetTriggerMask(esd->GetTriggerMask());
	242	header->SetTriggerCluster(esd->GetTriggerCluster());
	243	header->SetMagneticField(esd->GetMagneticField());
	244	header->SetZDCN1Energy(esd->GetZDCN1Energy());
	245	header->SetZDCP1Energy(esd->GetZDCP1Energy());
	246	header->SetZDCN2Energy(esd->GetZDCN2Energy());
	247	header->SetZDCP2Energy(esd->GetZDCP2Energy());
	248	header->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
	249	Float_t diamxy[2]={esd->GetDiamondX(),esd->GetDiamondY()};
	250	Float_t diamcov[3]; esd->GetDiamondCovXY(diamcov);
	251	header->SetDiamond(diamxy,diamcov);
	252	//
	253	//
	254	Int_t nVertices = 1 ;/* = prim. vtx*/;
	255	Int_t nCaloClus = esd->GetNumberOfCaloClusters();
	256
	257	AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
	258
	259	// Access to the AOD container of vertices
	260	TClonesArray &vertices = *(AODEvent()->GetVertices());
	261	Int_t jVertices=0;
	262
	263	// Add primary vertex. The primary tracks will be defined
	264	// after the loops on the composite objects (V0, cascades, kinks)
	265	const AliESDVertex *vtx = esd->GetPrimaryVertex();
	266
	267	vtx->GetXYZ(pos); // position
268	vtx->GetCovMatrix(covVtx); //covariance matrix
269
270	AliAODVertex * primary = new(vertices[jVertices++])
271	AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
272	primary->SetName(vtx->GetName());
273	primary->SetTitle(vtx->GetTitle());
274
275	// Access to the AOD container of clusters
276	TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
277	Int_t jClusters=0;
278
279	for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
280
281	AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
282
283	//Check which calorimeter information we want to keep.
284	if(!cluster->IsPHOS()) continue ;
285
286	Int_t id = cluster->GetID();
287	Float_t energy = cluster->E();
288	cluster->GetPosition(posF);
289
290	AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
291	AliAODCaloCluster(id,
292	0,
293	0x0,
294	energy,
295	posF,
296	NULL,
297	AliAODCluster::kPHOSNeutral);
298
299	caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
c8fe2783	300	cluster->GetDispersion(),
7a4cf423	301	cluster->GetM20(), cluster->GetM02(),
	302	cluster->GetEmcCpvDistance(),
	303	cluster->GetNExMax(),cluster->GetTOF()) ;
	304
c8fe2783	305	caloCluster->SetPIDFromESD(cluster->GetPID());
7a4cf423	306	caloCluster->SetNCells(cluster->GetNCells());
	307	caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
	308	caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
	309
	310	}
	311	caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
	312	// end of loop on calo clusters
	313
	314	// fill PHOS cell info
	315	if( esd->GetPHOSCells()) { // protection against missing ESD information
	316	AliESDCaloCells &esdPHcells = *(esd->GetPHOSCells());
	317	Int_t nPHcell = esdPHcells.GetNumberOfCells() ;
	318
	319	AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
	320	aodPHcells.CreateContainer(nPHcell);
c8fe2783	321	aodPHcells.SetType(AliVCaloCells::kPHOSCell);
7a4cf423	322	for (Int_t iCell = 0; iCell < nPHcell; iCell++) {
	323	aodPHcells.SetCell(iCell,esdPHcells.GetCellNumber(iCell),esdPHcells.GetAmplitude(iCell));
	324	}
	325	aodPHcells.Sort();
	326	}
	327
	328	}
	329
	330	//__________________________________________________
e727d2f7	331	void AliAnalysisTaskPHOSPi0CalibSelection::UserCreateOutputObjects()
	332	{
	333	//Create output container
	334	fOutputContainer = new TList();
	335
	336	char hname[128], htitl[128];
	337
	338	for(Int_t iMod=0; iMod<5; iMod++) {
	339	for(Int_t iX=0; iX<64; iX++) {
	340	for(Int_t iZ=0; iZ<56; iZ++) {
	341	sprintf(hname,"%d_%d_%d",iMod,iX,iZ);
	342	sprintf(htitl,"Two-gamma inv. mass for mod %d, cell (%d,%d)",iMod,iX,iZ);
	343	fHmpi0[iMod][iX][iZ] = new TH1F(hname,htitl,100,0.,300.);
	344	fOutputContainer->Add(fHmpi0[iMod][iX][iZ]);
	345	}
	346	}
	347	}
	348
	349	fHmgg = new TH1F("hmgg","2-cluster invariant mass",100,0.,300.);
	350	fOutputContainer->Add(fHmgg);
	351
f6b9f90c	352	fCalibData = new AliPHOSCalibData();
	353	fPhosGeo = AliPHOSGeometry::GetInstance("IHEP") ;
	354
e727d2f7	355	}
e727d2f7	356
7a4cf423	357	//__________________________________________________
e727d2f7	358	void AliAnalysisTaskPHOSPi0CalibSelection::UserExec(Option_t* /* option */)
	359	{
	360	//Analysis per event.
29a41d05	361	if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection <<< Event %d >>>\n",(Int_t)Entry());
7a4cf423	362
e727d2f7	363	AliAODEvent* aod = 0x0;
7a4cf423	364	if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) {
	365	//Input are ESDs
	366	aod = dynamic_cast<AliAODEvent*>(InputEvent());
	367	// Create new AOD with only CaloClusters and CaloCells
	368	if(fCopyAOD) CreateAODFromAOD();
	369	}
	370	else if(!strcmp(InputEvent()->GetName(),"AliESDEvent")) {
	371	//Input are ESDs
	372	aod = AODEvent();
	373	// Create AOD with CaloClusters and use it as input.
	374	// If filtering task is already executed, this is not needed.
	375	if(fCopyAOD) CreateAODFromESD();
	376	}
e727d2f7	377	else {
7a4cf423	378	printf("AliAnalysisTaskPHOSPi0CalibSelection: Unknown event type, STOP!\n");
7a4cf423	379	abort();
e727d2f7	380	}
	381
	382	Double_t v[] = {aod->GetVertex(0)->GetX(),aod->GetVertex(0)->GetY(),aod->GetVertex(0)->GetZ()}; //to check!!
	383	//aod->GetVertex()->GetXYZ(v) ;
	384	TVector3 vtx(v); //Check
	385
f6b9f90c	386	if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",vtx.X(),vtx.Y(),vtx.Z());
e727d2f7	387
e727d2f7	388	Int_t runNum = aod->GetRunNumber();
f6b9f90c	389	if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
4b4efabc	390
	391	//Get the matrix with geometry information
	392	//Still not implemented in AOD, just a workaround to be able to work at least with ESDs
f6b9f90c	393	if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) {
	394	if(DebugLevel() > 1)
	395	printf("AliAnalysisTaskPHOSPi0CalibSelection Use ideal geometry, values geometry matrix not kept in AODs.\n");
	396	}
4b4efabc	397	else{
f6b9f90c	398	if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection Load Misaligned matrices. \n");
4b4efabc	399	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
f6b9f90c	400	for(Int_t mod=0; mod<5; mod++){
	401	if(esd->GetPHOSMatrix(mod))
	402	fPhosGeo->SetMisalMatrix(esd->GetPHOSMatrix(mod),mod) ;
	403	}
e727d2f7	404	}
4b4efabc	405
f6b9f90c	406	if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
e727d2f7	407
	408	AliPHOSPIDv1 pid;
	409
	410	Int_t mod1 = -1;
	411	TLorentzVector p1;
	412	TLorentzVector p2;
	413	TLorentzVector p12;
	414	Int_t relid[4] ;
	415	Int_t maxId;
	416
	417	TRefArray * caloClustersArr = new TRefArray();
	418	aod->GetPHOSClusters(caloClustersArr);
	419
	420	const Int_t kNumberOfPhosClusters = caloClustersArr->GetEntries() ;
f6b9f90c	421	if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection CaloClusters: %d\n", kNumberOfPhosClusters);
e727d2f7	422
	423	// PHOS cells
	424	AliAODCaloCells *phsCells = aod->GetPHOSCells();
	425
	426	// loop over PHOS clusters
	427	for(Int_t iClu=0; iClu<kNumberOfPhosClusters; iClu++) {
	428
	429	AliAODCaloCluster c1 = (AliAODCaloCluster ) caloClustersArr->At(iClu);
c8fe2783	430	if(!c1->IsPHOS()) continue; // EMCAL cluster!
e727d2f7	431
	432	Float_t e1i = c1->E(); // cluster energy before correction
	433	if(e1i<fEmin) continue;
	434
	435	AliPHOSAodCluster clu1(*c1);
4b4efabc	436	clu1.Recalibrate(fCalibData, phsCells);
4b4efabc	437	clu1.EvalAll(fLogWeight,vtx);
e727d2f7	438	clu1.EnergyCorrection(&pid) ;
	439
	440	clu1.GetMomentum(p1,v);
	441
	442	MaxEnergyCellPos(phsCells,&clu1,maxId);
	443	fPhosGeo->AbsToRelNumbering(maxId, relid);
	444
	445	mod1 = relid[0]-1; // module
	446	Int_t iX = relid[2]-1; // cluster X-coord
	447	Int_t iZ = relid[3]-1; // cluster Z-coord
	448	Float_t e1ii = clu1.E(); // cluster energy after correction
	449
	450	for (Int_t jClu=iClu; jClu<kNumberOfPhosClusters; jClu++) {
	451	AliAODCaloCluster c2 = (AliAODCaloCluster ) caloClustersArr->At(jClu);
c8fe2783	452	if(!c2->IsPHOS()) continue; // EMCAL cluster!
e727d2f7	453	if(c2->IsEqual(c1)) continue;
	454
	455	Float_t e2i = c2->E();
	456	if(e2i<fEmin) continue;
	457
	458	AliPHOSAodCluster clu2(*c2);
4b4efabc	459	clu2.Recalibrate(fCalibData, phsCells);
4b4efabc	460	clu2.EvalAll(fLogWeight,vtx);
e727d2f7	461	clu2.EnergyCorrection(&pid) ;
	462
	463	clu2.GetMomentum(p2,v);
	464	Float_t e2ii = clu2.E();
	465
	466	p12 = p1+p2;
	467
	468	fHmgg->Fill(p12.M()*1000);
	469	fHmpi0[mod1][iX][iZ]->Fill(p12.M()*1000);
	470
f6b9f90c	471	if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection Mass in (mod%d,%d,%d): %.3f GeV E1_i=%f E1_ii=%f E2_i=%f E2_ii=%f\n",
e727d2f7	472	mod1,iX,iZ,p12.M(),e1i,e1ii,e2i,e2ii);
	473	}
	474
	475	} // end of loop over PHOS clusters
	476
	477	delete caloClustersArr;
	478	PostData(1,fOutputContainer);
	479	}
	480
7a4cf423	481	//__________________________________________________
e727d2f7	482	void AliAnalysisTaskPHOSPi0CalibSelection::MaxEnergyCellPos(AliAODCaloCells cells, AliAODCaloCluster clu, Int_t& maxId)
	483	{
	484	//For a given CaloCluster calculates the absId of the cell
	485	//with maximum energy deposit.
	486
	487	Double_t eMax = -111;
	488
	489	for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
	490	Int_t cellAbsId = clu->GetCellAbsId(iDig);
	491	Double_t eCell = cells->GetCellAmplitude(cellAbsId)*clu->GetCellAmplitudeFraction(iDig);
	492	if(eCell>eMax) {
	493	eMax = eCell;
	494	maxId = cellAbsId;
	495	}
	496	}
	497
	498	}
4b4efabc	499
7a4cf423	500	//__________________________________________________
4b4efabc	501	void AliAnalysisTaskPHOSPi0CalibSelection::SetCalibCorrections(AliPHOSCalibData* cdata)
	502	{
	503	//Set new correction factors (~1) to calibration coefficients, delete previous.
	504
	505	if(fCalibData) delete fCalibData;
	506	fCalibData = cdata;
	507
	508	}
	509