[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / UserTasks / AliAnalysisTaskSEPHOSpPbPi0.cxx

/**************************************************************************
 * Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * 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.                  *
 **************************************************************************/

///////////////////////////////////////////////////////////////////////////
//
// AliAnalysisTaskSE for the gamma and pi0 from pPb collision analysis
//
// Author: H-S. Zhu, hongsheng.zhu@cern.ch
//                   hszhu@iopp.ccnu.edu.cn
///////////////////////////////////////////////////////////////////////////

#include <TF1.h>
#include <TH2I.h>
#include <TList.h>
#include <TMath.h>
#include <TArray.h>
#include <TRandom.h>
#include <TVector3.h>
#include <TGeoManager.h>
#include <TClonesArray.h>

#include "AliInputEventHandler.h"
#include "AliAnalysisManager.h"
#include "AliMCEventHandler.h"
#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliVEvent.h"
#include "AliAODEvent.h"
#include "AliESDEvent.h"
#include "AliAODHeader.h"
#include "AliESDHeader.h"
#include "AliVCluster.h"
#include "AliVCaloCells.h"
#include "AliAODCaloCells.h"
#include "AliESDCaloCells.h"
#include "AliPHOSGeometry.h"
#include "AliPHOSGeoUtils.h"
#include "AliPHOSCalibData.h"
#include "AliPHOSAodCluster.h"
#include "AliPHOSEsdCluster.h"
#include "AliOADBContainer.h"
#include "AliPHOSpPbPi0Header.h"
#include "AliCaloClusterInfo.h"
#include "AliAnalysisTaskSEPHOSpPbPi0.h"

ClassImp(AliAnalysisTaskSEPHOSpPbPi0)

const Float_t AliAnalysisTaskSEPHOSpPbPi0::kLogWeight = 4.5;
Bool_t   AliAnalysisTaskSEPHOSpPbPi0::fgRemovePileup    = kFALSE;
Bool_t   AliAnalysisTaskSEPHOSpPbPi0::fgUseFiducialCut  = kFALSE;
Double_t AliAnalysisTaskSEPHOSpPbPi0::fgDecaliWidth     = 0.055;
Double_t AliAnalysisTaskSEPHOSpPbPi0::fgCuts[5]         = { 0.3, 2., 0.2, 2.5, 7e-8 };

//________________________________________________________________________
AliAnalysisTaskSEPHOSpPbPi0::AliAnalysisTaskSEPHOSpPbPi0():
  AliAnalysisTaskSE(), fIsMC(kFALSE), fCentralityBin(10), fBufferSize(10), fRunNumber(-1),
  fPHOSGeo(0), fPHOSCalibData(0), fListQA(0), fListRD(0), fListMC(0), fHeader(0), fCaloClArr(0)
{
  //
  // Default constructor
  //
  for (Int_t i=0; i<10; i++) { for (Int_t j=0; j<10; j++) fEventList[i][j] = 0; }

  // Init bad channel map
  for (Int_t i=0; i<5; i++)
    fPHOSBadMap[i] = new TH2I(Form("PHOS_BadMap_mod%d", i), Form("PHOS_BadMap_mod%d", i), 64, 0., 64., 56, 0., 56.);
}
//________________________________________________________________________
AliAnalysisTaskSEPHOSpPbPi0::AliAnalysisTaskSEPHOSpPbPi0(const char *name):
  AliAnalysisTaskSE(name), fIsMC(kFALSE), fCentralityBin(10), fBufferSize(10), fRunNumber(-1),
  fPHOSGeo(0), fPHOSCalibData(0), fListQA(0), fListRD(0), fListMC(0), fHeader(0), fCaloClArr(0)
{
  // Constructor
  for (Int_t i=0; i<10; i++) { for (Int_t j=0; j<10; j++) fEventList[i][j] = 0; }

  // Init bad channel map
  for (Int_t i=0; i<5; i++)
    fPHOSBadMap[i] = new TH2I(Form("PHOS_BadMap_mod%d", i), Form("PHOS_BadMap_mod%d", i), 64, 0., 64., 56, 0., 56.);

  DefineOutput(1, TList::Class());
  DefineOutput(2, TList::Class());
  DefineOutput(3, TList::Class());
}

//_____________________________________________________________________________
AliAnalysisTaskSEPHOSpPbPi0::~AliAnalysisTaskSEPHOSpPbPi0()
{
  //
  // Default destructor
  //
  if (fListQA)      { delete fListQA;      fListQA     = NULL; }
  if (fListRD)      { delete fListRD;      fListRD     = NULL; }
  if (fListMC)      { delete fListMC;      fListMC     = NULL; }
  if (fHeader)      { delete fHeader;      fHeader     = NULL; }
  if (fCaloClArr)   { delete fCaloClArr;   fCaloClArr  = NULL; }

}
//________________________________________________________________________
void AliAnalysisTaskSEPHOSpPbPi0::UserCreateOutputObjects()
{
  // Create the output container

  AliPHOSpPbPi0Header::SetIsMC(fIsMC);
  AliPHOSpPbPi0Header::SetUseFiducialCut(fgUseFiducialCut);

  if (!fHeader)     fHeader     = new AliPHOSpPbPi0Header();
  if (!fCaloClArr)  fCaloClArr  = new TClonesArray("AliCaloClusterInfo", 0);
  if (!fListQA)     fListQA     = new TList();
  if (!fListRD)     fListRD     = new TList();
  if (!fListMC)     fListMC     = new TList();

  fHeader->SetNCent(fCentralityBin.GetSize()-1);
  fHeader->CreateHistograms(fListQA, fListRD, fListMC);

  // Post output data.
  PostData(1, fListQA);
  PostData(2, fListRD);
  if (fIsMC) PostData(3, fListMC);

  return;
}

//________________________________________________________________________
void AliAnalysisTaskSEPHOSpPbPi0::UserExec(Option_t *) 
{
  // Main loop
  // Called for each event

  if (fIsMC) {
    if (MCEvent()) {
      if (MCEvent()->GetNumberOfTracks()<=0)
           { AliError("MC event not found. Nothing done!"); return; }
    } else { AliError("MC event not found. Nothing done!"); return; }
  }

  AliAODEvent *aod   = 0x0;
  AliESDEvent *esd   = 0x0;

  if (((TString)fInputEvent->IsA()->GetName())=="AliAODEvent") {
    aod = dynamic_cast<AliAODEvent*>(fInputEvent);
    if (!aod) { AliError("AOD event not found. Nothing done!");   return; }
    if (!fIsMC && (aod->GetHeader()->GetEventType()!=7))          return; // check event type; should be PHYSICS = 7 for data and 0 for MC
  } else {
    esd = dynamic_cast<AliESDEvent*>(fInputEvent);
    if (!esd) { AliError("ESD event not found. Nothing done!");   return; }
    if (!fIsMC && (esd->GetHeader()->GetEventType()!=7))          return; // check event type; should be PHYSICS = 7 for data and 0 for MC
  }

  // Fill Event info
  fHeader->SetEventInfo(fInputHandler);
  fHeader->FillHistosEvent(fListQA);

  // PHOS Geometry and Misalignment initialization at the first time it runs
  if(fRunNumber != fInputEvent->GetRunNumber()) {
    fRunNumber = fInputEvent->GetRunNumber();
    fHeader->SetIspARun(fRunNumber>195344 && fRunNumber<197388);  // flag for pA collisions
    PHOSInitialize(esd);
  }

  // Event Selection
  if (!fHeader->IsSelected())                     return;
  if (fgRemovePileup && fHeader->IsPileupSPD())   return;

  // Fill PHOS cells QA histograms
  fHeader->FillHistosCaloCellsQA(fListQA, fInputEvent->GetPHOSCells(), fPHOSGeo);

  // Fill PHOS cluster Clones Array
  FillCaloClusterInfo(aod, esd);
  aod = 0x0;

  if (!fCaloClArr->GetEntriesFast()) return;

  // vertex bining and centrality bining
  Int_t zvtx = (Int_t)((fHeader->Vz() + 10.)/2.);   if (zvtx<0) zvtx = 0;   if (zvtx>9) zvtx = 9;
  Int_t cent = TMath::BinarySearch<Float_t>(fCentralityBin.GetSize()-1, fCentralityBin.GetArray(), fHeader->Centrality());
  if (!fEventList[zvtx][cent]) fEventList[zvtx][cent] = new TList(); 
  TList *eventList = fEventList[zvtx][cent];

  // Fill cluster histograms
  fHeader->FillHistosCaloCluster(fListQA, fCaloClArr, cent);

  // Fill pi0 histograms
  fHeader->FillHistosPi0(fListRD, fCaloClArr, cent);

  // Fill mixed pi0 histograms
  fHeader->FillHistosMixPi0(fListRD, fCaloClArr, eventList, cent);

  // Fill MC info
  AliStack *stack = 0x0;
  if (fIsMC) {
    if (esd) {
      if (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()) {
        if (static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent())
          stack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
      }
      fHeader->FillHistosMC(fListMC, stack, fPHOSGeo, cent);
    } else
      fHeader->FillHistosMC(fListMC, MCEvent(), fPHOSGeo, cent);
  }
  esd = 0x0;

  // Fill event list for mixing
  if (fCaloClArr->GetEntriesFast()>0) {
    eventList->AddFirst(fCaloClArr);   fCaloClArr = 0x0;
    if (eventList->GetSize()>fBufferSize[cent]) { // Remove redundant events
      TClonesArray *tmp = static_cast<TClonesArray*>(eventList->Last());
      eventList->RemoveLast();
      delete tmp;
    }
  }

  return;
}

//________________________________________________________________________
void AliAnalysisTaskSEPHOSpPbPi0::Terminate(Option_t *) 
{
  // Terminate analysis

  // add the correction matrix

  return;
}

//________________________________________________________________________
void AliAnalysisTaskSEPHOSpPbPi0::PHOSInitialize(AliESDEvent* const esd)
{
  // Initialize PHOS Geometry ,misalignment and calibration 
//TGeoManager::Import("$ALICE_ROOT/test/QA/geometry.root");

  fPHOSGeo =  AliPHOSGeometry::GetInstance("IHEP");
  AliOADBContainer geomContainer("phosGeo");
  geomContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSGeometry.root","PHOSRotationMatrixes");
  TObjArray *matrixes = (TObjArray*)geomContainer.GetObject(fRunNumber,"PHOSRotationMatrixes");
  for (Int_t mod=0; mod<5; mod++) {
    if (!matrixes->At(mod)) continue;
    else fPHOSGeo->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod);
  }

  // sets the PHOS Misalignment vertex if ESD
  if (esd) {
    for (Int_t mod=0; mod<5; mod++) {
      const TGeoHMatrix* modMatrix = fInputEvent->GetPHOSMatrix(mod);
      if (!modMatrix) continue;
      else fPHOSGeo->SetMisalMatrix(modMatrix, mod);
    }
  }

  TRandom random;   random.SetSeed(0); // the seed is set to the current  machine clock
  fPHOSCalibData = new AliPHOSCalibData();   fPHOSCalibData->SetName("PHOSCalibData");
  for(Int_t iMod=1; iMod<6; iMod++) {
    for(Int_t iCol=1; iCol<57; iCol++) {
      for(Int_t iRow=1; iRow<65; iRow++) {
        Float_t value = fIsMC ? random.Gaus(1., fgDecaliWidth) : 0.;  // ADC channel Emc
        fPHOSCalibData->SetADCchannelEmc(iMod, iCol, iRow, value);
      }
    }
  }
}

//________________________________________________________________________
void AliAnalysisTaskSEPHOSpPbPi0::FillCaloClusterInfo(AliAODEvent* const aod, AliESDEvent* const esd)
{
  // Fill calo cluster info
  if (fCaloClArr) fCaloClArr->Clear();
  else fCaloClArr = new TClonesArray("AliCaloClusterInfo", 0);

  Int_t    nclsts       = fInputEvent->GetNumberOfCaloClusters();
  Int_t    countN       = fCaloClArr->GetEntriesFast();
  Int_t    relID[4]     = {0,0,0,0 }; // module = relID[0]; cellX = relID[2]; cellZ = relID[3];
  Float_t  position[3]  = {0.,0.,0.};
  Double_t vtx[3]       = {0,0,0};   fHeader->GetXYZ(vtx);   TVector3 vtxVector(vtx); 
  Double_t magfield     = fInputEvent->GetMagneticField();
  TF1 *nonLinCorr = new TF1("Non-linear", "0.0241+1.0504*x+0.000249*x*x", 0., 50.); // GCB's non-linear correction function

  TClonesArray &caloRef = *fCaloClArr;
  TLorentzVector momentum;
  AliVCluster        *clust       = 0x0;
  AliCaloClusterInfo *caloCluster = 0x0;
  for (Int_t iclst=0; iclst<nclsts; iclst++) {  // loop over all clusters
    clust = fInputEvent->GetCaloCluster(iclst);
    if (!(clust && clust->IsPHOS() && clust->E()>fgCuts[0]))                 { clust=0; continue; }
    if (!(clust->GetNCells()>(Int_t)fgCuts[1] && clust->GetM02()>fgCuts[2])) { clust=0; continue; } // To remove exotic clusters
    if (!(clust->GetDistanceToBadChannel()>fgCuts[3]))                       { clust=0; continue; }
//  if (!(TMath::Abs(clust->GetTOF()<fgCuts[4])))                            { clust=0; continue; } // remove clusters from pileup or same benches
    clust->GetPosition(position); TVector3 global(position); fPHOSGeo->GlobalPos2RelId(global,relID);
    if (relID[0] == 2)                                                       { clust=0; continue; } // !remove module 2
    if (!IsGoodCaloCluster(relID[0], relID[2], relID[3]))                    { clust=0; continue; } // calo cluster selection 

    caloCluster = new AliCaloClusterInfo(clust, esd, relID, magfield);
    if (fgUseFiducialCut && !caloCluster->IsInFiducialRegion(relID[2], relID[3])) { delete caloCluster; caloCluster=0; clust=0; continue; } // Fiducial cut 

    if (aod) { // TODO recalibration for AOD
      if (fIsMC) {
        AliPHOSAodCluster aodClust(*(AliAODCaloCluster*) (clust));
        aodClust.Recalibrate(fPHOSCalibData, aod->GetPHOSCells());
        aodClust.EvalAll(kLogWeight, vtxVector);        // recalculate all cluster parameters
        aodClust.SetE(nonLinCorr->Eval(aodClust.E()));  // non-linear correction
        aodClust.GetMomentum(momentum, vtx);
      } else clust->GetMomentum(momentum, vtx);
      caloCluster->SetLorentzVector(momentum);
    } else { // for ESD
      AliPHOSEsdCluster esdClust(*(AliESDCaloCluster*) (clust));
      esdClust.Recalibrate(fPHOSCalibData, esd->GetPHOSCells());
      esdClust.EvalAll(kLogWeight, vtxVector);       // recalculate all cluster parameters
      esdClust.SetE(nonLinCorr->Eval(esdClust.E()));  // non-linear correction
      esdClust.GetMomentum(momentum, vtx);
      caloCluster->SetLorentzVector(momentum);
    }

    if (fIsMC && !(momentum.E()>fgCuts[0])) { delete caloCluster; caloCluster=0; clust=0; continue; } // check for MC

    Double_t disp = caloCluster->TestDisp();
    if (disp < 2.5*2.5) caloCluster->SetPIDBit(BIT(1)); // set Disp1 bit
    if (disp < 1.5*1.5) caloCluster->SetPIDBit(BIT(3)); // set Disp2 bit

    clust = 0;

    new(caloRef[countN++]) AliCaloClusterInfo(*caloCluster);
    delete caloCluster; caloCluster=0;
  }  // end loop of all clusters
}

//________________________________________________________________________
Bool_t AliAnalysisTaskSEPHOSpPbPi0::IsGoodCaloCluster(Int_t iMod, Int_t cellX, Int_t cellZ)
{
  // !Check whether this cluster is not in bad channel

  if (!(iMod>0 && iMod<5 && fPHOSBadMap[iMod]))          return kTRUE;   //No bad maps for this Module
  if (fPHOSBadMap[iMod]->GetBinContent(cellX,cellZ)>0)   return kFALSE;

  return kTRUE;
}
Commit	Line	Data
a01cd2c9	1	/**************************************************************************
	2	* Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	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	// AliAnalysisTaskSE for the gamma and pi0 from pPb collision analysis
	19	//
	20	// Author: H-S. Zhu, hongsheng.zhu@cern.ch
	21	// hszhu@iopp.ccnu.edu.cn
	22	///////////////////////////////////////////////////////////////////////////
	23
e824e25b	24	#include <TF1.h>
a01cd2c9	25	#include <TH2I.h>
	26	#include <TList.h>
	27	#include <TMath.h>
	28	#include <TArray.h>
e824e25b	29	#include <TRandom.h>
	30	#include <TVector3.h>
	31	#include <TGeoManager.h>
a01cd2c9	32	#include <TClonesArray.h>
a01cd2c9	33
e824e25b	34	#include "AliInputEventHandler.h"
a01cd2c9	35	#include "AliAnalysisManager.h"
	36	#include "AliMCEventHandler.h"
	37	#include "AliMCEvent.h"
	38	#include "AliStack.h"
	39	#include "AliVEvent.h"
a01cd2c9	40	#include "AliAODEvent.h"
e824e25b	41	#include "AliESDEvent.h"
a01cd2c9	42	#include "AliAODHeader.h"
e824e25b	43	#include "AliESDHeader.h"
a01cd2c9	44	#include "AliVCluster.h"
a01cd2c9	45	#include "AliVCaloCells.h"
e824e25b	46	#include "AliAODCaloCells.h"
e824e25b	47	#include "AliESDCaloCells.h"
a01cd2c9	48	#include "AliPHOSGeometry.h"
e824e25b	49	#include "AliPHOSGeoUtils.h"
	50	#include "AliPHOSCalibData.h"
	51	#include "AliPHOSAodCluster.h"
	52	#include "AliPHOSEsdCluster.h"
a01cd2c9	53	#include "AliOADBContainer.h"
a01cd2c9	54	#include "AliPHOSpPbPi0Header.h"
	55	#include "AliCaloClusterInfo.h"
	56	#include "AliAnalysisTaskSEPHOSpPbPi0.h"
	57
	58	ClassImp(AliAnalysisTaskSEPHOSpPbPi0)
	59
a91370e6	60	const Float_t AliAnalysisTaskSEPHOSpPbPi0::kLogWeight = 4.5;
e824e25b	61	Bool_t AliAnalysisTaskSEPHOSpPbPi0::fgRemovePileup = kFALSE;
	62	Bool_t AliAnalysisTaskSEPHOSpPbPi0::fgUseFiducialCut = kFALSE;
	63	Double_t AliAnalysisTaskSEPHOSpPbPi0::fgDecaliWidth = 0.055;
	64	Double_t AliAnalysisTaskSEPHOSpPbPi0::fgCuts[5] = { 0.3, 2., 0.2, 2.5, 7e-8 };
a01cd2c9	65
	66	//________________________________________________________________________
	67	AliAnalysisTaskSEPHOSpPbPi0::AliAnalysisTaskSEPHOSpPbPi0():
e824e25b	68	AliAnalysisTaskSE(), fIsMC(kFALSE), fCentralityBin(10), fBufferSize(10), fRunNumber(-1),
e824e25b	69	fPHOSGeo(0), fPHOSCalibData(0), fListQA(0), fListRD(0), fListMC(0), fHeader(0), fCaloClArr(0)
a01cd2c9	70	{
	71	//
	72	// Default constructor
	73	//
	74	for (Int_t i=0; i<10; i++) { for (Int_t j=0; j<10; j++) fEventList[i][j] = 0; }
	75
e824e25b	76	// Init bad channel map
a01cd2c9	77	for (Int_t i=0; i<5; i++)
	78	fPHOSBadMap[i] = new TH2I(Form("PHOS_BadMap_mod%d", i), Form("PHOS_BadMap_mod%d", i), 64, 0., 64., 56, 0., 56.);
	79	}
	80	//________________________________________________________________________
	81	AliAnalysisTaskSEPHOSpPbPi0::AliAnalysisTaskSEPHOSpPbPi0(const char *name):
e824e25b	82	AliAnalysisTaskSE(name), fIsMC(kFALSE), fCentralityBin(10), fBufferSize(10), fRunNumber(-1),
e824e25b	83	fPHOSGeo(0), fPHOSCalibData(0), fListQA(0), fListRD(0), fListMC(0), fHeader(0), fCaloClArr(0)
a01cd2c9	84	{
	85	// Constructor
	86	for (Int_t i=0; i<10; i++) { for (Int_t j=0; j<10; j++) fEventList[i][j] = 0; }
	87
e824e25b	88	// Init bad channel map
a01cd2c9	89	for (Int_t i=0; i<5; i++)
	90	fPHOSBadMap[i] = new TH2I(Form("PHOS_BadMap_mod%d", i), Form("PHOS_BadMap_mod%d", i), 64, 0., 64., 56, 0., 56.);
	91
	92	DefineOutput(1, TList::Class());
0f8c33c1	93	DefineOutput(2, TList::Class());
0f8c33c1	94	DefineOutput(3, TList::Class());
a01cd2c9	95	}
	96
	97	//_____________________________________________________________________________
	98	AliAnalysisTaskSEPHOSpPbPi0::~AliAnalysisTaskSEPHOSpPbPi0()
	99	{
	100	//
	101	// Default destructor
	102	//
e824e25b	103	if (fListQA) { delete fListQA; fListQA = NULL; }
e824e25b	104	if (fListRD) { delete fListRD; fListRD = NULL; }
0f8c33c1	105	if (fListMC) { delete fListMC; fListMC = NULL; }
	106	if (fHeader) { delete fHeader; fHeader = NULL; }
	107	if (fCaloClArr) { delete fCaloClArr; fCaloClArr = NULL; }
a01cd2c9	108
e824e25b	109	}
a01cd2c9	110	//________________________________________________________________________
	111	void AliAnalysisTaskSEPHOSpPbPi0::UserCreateOutputObjects()
	112	{
	113	// Create the output container
e824e25b	114
	115	AliPHOSpPbPi0Header::SetIsMC(fIsMC);
	116	AliPHOSpPbPi0Header::SetUseFiducialCut(fgUseFiducialCut);
a01cd2c9	117
0f8c33c1	118	if (!fHeader) fHeader = new AliPHOSpPbPi0Header();
0f8c33c1	119	if (!fCaloClArr) fCaloClArr = new TClonesArray("AliCaloClusterInfo", 0);
e824e25b	120	if (!fListQA) fListQA = new TList();
e824e25b	121	if (!fListRD) fListRD = new TList();
0f8c33c1	122	if (!fListMC) fListMC = new TList();
a01cd2c9	123
0f8c33c1	124	fHeader->SetNCent(fCentralityBin.GetSize()-1);
e824e25b	125	fHeader->CreateHistograms(fListQA, fListRD, fListMC);
a01cd2c9	126
a01cd2c9	127	// Post output data.
e824e25b	128	PostData(1, fListQA);
	129	PostData(2, fListRD);
	130	if (fIsMC) PostData(3, fListMC);
0f8c33c1	131
a01cd2c9	132	return;
	133	}
	134
	135	//________________________________________________________________________
	136	void AliAnalysisTaskSEPHOSpPbPi0::UserExec(Option_t *)
	137	{
	138	// Main loop
	139	// Called for each event
	140
	141	if (fIsMC) {
	142	if (MCEvent()) {
	143	if (MCEvent()->GetNumberOfTracks()<=0)
	144	{ AliError("MC event not found. Nothing done!"); return; }
	145	} else { AliError("MC event not found. Nothing done!"); return; }
	146	}
	147
e824e25b	148	AliAODEvent *aod = 0x0;
e824e25b	149	AliESDEvent *esd = 0x0;
a01cd2c9	150
	151	if (((TString)fInputEvent->IsA()->GetName())=="AliAODEvent") {
	152	aod = dynamic_cast<AliAODEvent*>(fInputEvent);
	153	if (!aod) { AliError("AOD event not found. Nothing done!"); return; }
	154	if (!fIsMC && (aod->GetHeader()->GetEventType()!=7)) return; // check event type; should be PHYSICS = 7 for data and 0 for MC
a01cd2c9	155	} else {
	156	esd = dynamic_cast<AliESDEvent*>(fInputEvent);
	157	if (!esd) { AliError("ESD event not found. Nothing done!"); return; }
	158	if (!fIsMC && (esd->GetHeader()->GetEventType()!=7)) return; // check event type; should be PHYSICS = 7 for data and 0 for MC
a01cd2c9	159	}
	160
	161	// Fill Event info
	162	fHeader->SetEventInfo(fInputHandler);
e824e25b	163	fHeader->FillHistosEvent(fListQA);
a01cd2c9	164
	165	// PHOS Geometry and Misalignment initialization at the first time it runs
	166	if(fRunNumber != fInputEvent->GetRunNumber()) {
	167	fRunNumber = fInputEvent->GetRunNumber();
e824e25b	168	fHeader->SetIspARun(fRunNumber>195344 && fRunNumber<197388); // flag for pA collisions
a01cd2c9	169	PHOSInitialize(esd);
	170	}
	171
e824e25b	172	// Event Selection
	173	if (!fHeader->IsSelected()) return;
	174	if (fgRemovePileup && fHeader->IsPileupSPD()) return;
	175
a01cd2c9	176	// Fill PHOS cells QA histograms
e824e25b	177	fHeader->FillHistosCaloCellsQA(fListQA, fInputEvent->GetPHOSCells(), fPHOSGeo);
a01cd2c9	178
a01cd2c9	179	// Fill PHOS cluster Clones Array
e824e25b	180	FillCaloClusterInfo(aod, esd);
e824e25b	181	aod = 0x0;
a01cd2c9	182
a01cd2c9	183	if (!fCaloClArr->GetEntriesFast()) return;
0f8c33c1	184
	185	// vertex bining and centrality bining
	186	Int_t zvtx = (Int_t)((fHeader->Vz() + 10.)/2.); if (zvtx<0) zvtx = 0; if (zvtx>9) zvtx = 9;
	187	Int_t cent = TMath::BinarySearch<Float_t>(fCentralityBin.GetSize()-1, fCentralityBin.GetArray(), fHeader->Centrality());
a01cd2c9	188	if (!fEventList[zvtx][cent]) fEventList[zvtx][cent] = new TList();
	189	TList *eventList = fEventList[zvtx][cent];
	190
	191	// Fill cluster histograms
e824e25b	192	fHeader->FillHistosCaloCluster(fListQA, fCaloClArr, cent);
a01cd2c9	193
a01cd2c9	194	// Fill pi0 histograms
e824e25b	195	fHeader->FillHistosPi0(fListRD, fCaloClArr, cent);
a01cd2c9	196
a01cd2c9	197	// Fill mixed pi0 histograms
e824e25b	198	fHeader->FillHistosMixPi0(fListRD, fCaloClArr, eventList, cent);
a01cd2c9	199
a01cd2c9	200	// Fill MC info
0f8c33c1	201	AliStack *stack = 0x0;
a01cd2c9	202	if (fIsMC) {
a01cd2c9	203	if (esd) {
0f8c33c1	204	if (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()) {
e824e25b	205	if (static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent())
e824e25b	206	stack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
a01cd2c9	207	}
0f8c33c1	208	fHeader->FillHistosMC(fListMC, stack, fPHOSGeo, cent);
a01cd2c9	209	} else
0f8c33c1	210	fHeader->FillHistosMC(fListMC, MCEvent(), fPHOSGeo, cent);
a01cd2c9	211	}
e824e25b	212	esd = 0x0;
a01cd2c9	213
a01cd2c9	214	// Fill event list for mixing
0f8c33c1	215	if (fCaloClArr->GetEntriesFast()>0) {
e824e25b	216	eventList->AddFirst(fCaloClArr); fCaloClArr = 0x0;
0f8c33c1	217	if (eventList->GetSize()>fBufferSize[cent]) { // Remove redundant events
0f8c33c1	218	TClonesArray tmp = static_cast<TClonesArray>(eventList->Last());
a01cd2c9	219	eventList->RemoveLast();
0f8c33c1	220	delete tmp;
a01cd2c9	221	}
	222	}
	223
	224	return;
0f8c33c1	225	}
a01cd2c9	226
	227	//________________________________________________________________________
	228	void AliAnalysisTaskSEPHOSpPbPi0::Terminate(Option_t *)
	229	{
	230	// Terminate analysis
	231
	232	// add the correction matrix
	233
	234	return;
	235	}
	236
	237	//________________________________________________________________________
	238	void AliAnalysisTaskSEPHOSpPbPi0::PHOSInitialize(AliESDEvent* const esd)
	239	{
e824e25b	240	// Initialize PHOS Geometry ,misalignment and calibration
	241	//TGeoManager::Import("$ALICE_ROOT/test/QA/geometry.root");
	242
	243	fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP");
a01cd2c9	244	AliOADBContainer geomContainer("phosGeo");
	245	geomContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSGeometry.root","PHOSRotationMatrixes");
	246	TObjArray matrixes = (TObjArray)geomContainer.GetObject(fRunNumber,"PHOSRotationMatrixes");
a01cd2c9	247	for (Int_t mod=0; mod<5; mod++) {
	248	if (!matrixes->At(mod)) continue;
	249	else fPHOSGeo->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod);
	250	}
	251
	252	// sets the PHOS Misalignment vertex if ESD
	253	if (esd) {
	254	for (Int_t mod=0; mod<5; mod++) {
	255	const TGeoHMatrix* modMatrix = fInputEvent->GetPHOSMatrix(mod);
	256	if (!modMatrix) continue;
	257	else fPHOSGeo->SetMisalMatrix(modMatrix, mod);
	258	}
	259	}
e824e25b	260
	261	TRandom random; random.SetSeed(0); // the seed is set to the current machine clock
	262	fPHOSCalibData = new AliPHOSCalibData(); fPHOSCalibData->SetName("PHOSCalibData");
	263	for(Int_t iMod=1; iMod<6; iMod++) {
	264	for(Int_t iCol=1; iCol<57; iCol++) {
	265	for(Int_t iRow=1; iRow<65; iRow++) {
	266	Float_t value = fIsMC ? random.Gaus(1., fgDecaliWidth) : 0.; // ADC channel Emc
	267	fPHOSCalibData->SetADCchannelEmc(iMod, iCol, iRow, value);
	268	}
	269	}
	270	}
a01cd2c9	271	}
	272
	273	//________________________________________________________________________
e824e25b	274	void AliAnalysisTaskSEPHOSpPbPi0::FillCaloClusterInfo(AliAODEvent* const aod, AliESDEvent* const esd)
a01cd2c9	275	{
	276	// Fill calo cluster info
	277	if (fCaloClArr) fCaloClArr->Clear();
	278	else fCaloClArr = new TClonesArray("AliCaloClusterInfo", 0);
	279
0f8c33c1	280	Int_t nclsts = fInputEvent->GetNumberOfCaloClusters();
e824e25b	281	Int_t countN = fCaloClArr->GetEntriesFast();
	282	Int_t relID[4] = {0,0,0,0 }; // module = relID[0]; cellX = relID[2]; cellZ = relID[3];
	283	Float_t position[3] = {0.,0.,0.};
	284	Double_t vtx[3] = {0,0,0}; fHeader->GetXYZ(vtx); TVector3 vtxVector(vtx);
	285	Double_t magfield = fInputEvent->GetMagneticField();
	286	TF1 nonLinCorr = new TF1("Non-linear", "0.0241+1.0504x+0.000249xx", 0., 50.); // GCB's non-linear correction function
	287
0f8c33c1	288	TClonesArray &caloRef = *fCaloClArr;
0f8c33c1	289	TLorentzVector momentum;
e824e25b	290	AliVCluster *clust = 0x0;
e824e25b	291	AliCaloClusterInfo *caloCluster = 0x0;
a01cd2c9	292	for (Int_t iclst=0; iclst<nclsts; iclst++) { // loop over all clusters
a01cd2c9	293	clust = fInputEvent->GetCaloCluster(iclst);
e824e25b	294	if (!(clust && clust->IsPHOS() && clust->E()>fgCuts[0])) { clust=0; continue; }
	295	if (!(clust->GetNCells()>(Int_t)fgCuts[1] && clust->GetM02()>fgCuts[2])) { clust=0; continue; } // To remove exotic clusters
	296	if (!(clust->GetDistanceToBadChannel()>fgCuts[3])) { clust=0; continue; }
	297	// if (!(TMath::Abs(clust->GetTOF()<fgCuts[4]))) { clust=0; continue; } // remove clusters from pileup or same benches
	298	clust->GetPosition(position); TVector3 global(position); fPHOSGeo->GlobalPos2RelId(global,relID);
	299	if (relID[0] == 2) { clust=0; continue; } // !remove module 2
	300	if (!IsGoodCaloCluster(relID[0], relID[2], relID[3])) { clust=0; continue; } // calo cluster selection
	301
	302	caloCluster = new AliCaloClusterInfo(clust, esd, relID, magfield);
	303	if (fgUseFiducialCut && !caloCluster->IsInFiducialRegion(relID[2], relID[3])) { delete caloCluster; caloCluster=0; clust=0; continue; } // Fiducial cut
	304
	305	if (aod) { // TODO recalibration for AOD
	306	if (fIsMC) {
	307	AliPHOSAodCluster aodClust((AliAODCaloCluster) (clust));
	308	aodClust.Recalibrate(fPHOSCalibData, aod->GetPHOSCells());
	309	aodClust.EvalAll(kLogWeight, vtxVector); // recalculate all cluster parameters
	310	aodClust.SetE(nonLinCorr->Eval(aodClust.E())); // non-linear correction
	311	aodClust.GetMomentum(momentum, vtx);
	312	} else clust->GetMomentum(momentum, vtx);
	313	caloCluster->SetLorentzVector(momentum);
	314	} else { // for ESD
	315	AliPHOSEsdCluster esdClust((AliESDCaloCluster) (clust));
	316	esdClust.Recalibrate(fPHOSCalibData, esd->GetPHOSCells());
	317	esdClust.EvalAll(kLogWeight, vtxVector); // recalculate all cluster parameters
	318	esdClust.SetE(nonLinCorr->Eval(esdClust.E())); // non-linear correction
	319	esdClust.GetMomentum(momentum, vtx);
	320	caloCluster->SetLorentzVector(momentum);
	321	}
	322
	323	if (fIsMC && !(momentum.E()>fgCuts[0])) { delete caloCluster; caloCluster=0; clust=0; continue; } // check for MC
	324
	325	Double_t disp = caloCluster->TestDisp();
	326	if (disp < 2.5*2.5) caloCluster->SetPIDBit(BIT(1)); // set Disp1 bit
	327	if (disp < 1.5*1.5) caloCluster->SetPIDBit(BIT(3)); // set Disp2 bit
a01cd2c9	328
	329	clust = 0;
	330
	331	new(caloRef[countN++]) AliCaloClusterInfo(*caloCluster);
a01cd2c9	332	delete caloCluster; caloCluster=0;
	333	} // end loop of all clusters
	334	}
	335
	336	//________________________________________________________________________
	337	Bool_t AliAnalysisTaskSEPHOSpPbPi0::IsGoodCaloCluster(Int_t iMod, Int_t cellX, Int_t cellZ)
	338	{
	339	// !Check whether this cluster is not in bad channel
	340
	341	if (!(iMod>0 && iMod<5 && fPHOSBadMap[iMod])) return kTRUE; //No bad maps for this Module
	342	if (fPHOSBadMap[iMod]->GetBinContent(cellX,cellZ)>0) return kFALSE;
	343
	344	return kTRUE;
	345	}