[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"
#include "AliCentrality.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());
  
  if(!aod) {
    printf("AliAnalysisTaskPHOSPi0CalibSelection::CreateAODFromAOD() - This event does not contain AODs?");
    return;
  }
  
  // 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());

  if(headerin->GetCentrality()){
    header->SetCentrality(new AliCentrality(*(headerin->GetCentralityP())));
  } else {
    header->SetCentrality(0);
  }

  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());
  
  if(!esd) {
    printf("AliAnalysisTaskPHOSPi0CalibSelection::CreateAODFromESD() - This event does not contain AODs?");
    return;
  }
  
  // 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(0);            // 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();
  const Int_t buffersize = 255;
  char hname[buffersize], htitl[buffersize];
  
  for(Int_t iMod=0; iMod<5; iMod++) {
    for(Int_t iX=0; iX<64; iX++) {
      for(Int_t iZ=0; iZ<56; iZ++) {
        snprintf(hname,buffersize,"%d_%d_%d",iMod,iX,iZ);
        snprintf(htitl,buffersize,"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());
    if(!aod) {
      printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - This event does not contain AODs?");
      return;
    }    
    
    // Create new AOD with only CaloClusters and CaloCells
    if(fCopyAOD) CreateAODFromAOD();
  }
  else  if(!strcmp(InputEvent()->GetName(),"AliESDEvent")) {
    //Input are ESDs
    aod = AODEvent();
    if(!aod) {
      printf("AliAnalysisTaskPHOSPi0CalibSelection::UserExec() - This event does not contain AODs?");
      return;
    } 
    
    // 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()) ;
    if(!esd) {
       printf("AliAnalysisTaskPHOSPi0CalibSelection::UserExec() - This event does not contain ESDs?");
      return;
    }
	  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);

  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() ;

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