[u/mrichter/AliRoot.git] / PWGLF / totEt / AliAnalysisEmEtReconstructed.cxx

//_________________________________________________________________________
//  Utility Class for transverse energy studies
//  Base class for MC analysis
//  - MC output
//  implementation file
//
//*-- Author: Marcelo G. Munhoz (USP)
	      //_________________________________________________________________________

#include "AliAnalysisEmEtReconstructed.h"
#include "AliAnalysisEtCuts.h"
#include "AliESDtrack.h"
#include "AliStack.h"
#include "AliVEvent.h"
#include "AliMCEvent.h"
#include "AliESDEvent.h"
#include "TH2F.h"
#include "TParticle.h"
#include "AliGenHijingEventHeader.h"
#include "AliGenPythiaEventHeader.h"
#include "TList.h"
#include "AliESDCaloCluster.h"
#include "TGeoGlobalMagField.h"
#include "AliMagF.h"
#include "AliEMCALTrack.h"
#include "AliESDtrackCuts.h"
#include "AliEMCALGeometry.h"
#include "AliExternalTrackParam.h"
#include "AliTrackerBase.h"
#include "TGeoManager.h"
#include "AliCentrality.h"

  using namespace std;

ClassImp(AliAnalysisEmEtReconstructed);


// ctor
AliAnalysisEmEtReconstructed::AliAnalysisEmEtReconstructed():AliAnalysisEtReconstructed()
							    ,fResCut(0)
							    ,fAllRectotETDep(0)
							    ,fElectronMatchtotETDep(0)
							    ,fNeutralRectotET(0)
							    ,fTotEMRectotET(0)

							    ,fMuonMatchtotETDep(0), fPionMatchtotETDep(0), fKaonMatchtotETDep(0), fProtonMatchtotETDep(0)
							    ,fTotChargedMatchtotETDep(0)

							    ,fTotalRectotETDep(0)

							    ,fESD(0)
// 							    ,fGeoUt(0)

							    ,fHistAllRecETDep(0) 
							    ,fHistAllRec(0) 
							    ,fHistAllRectotETDep(0) 

							    ,fHistElectronRecETDep(0) 
							    ,fHistElectronRec(0) 
							    ,fHistElectronMatchtotETDep(0) 
							    ,fHistElectronRecdEdxP(0)

							    ,fHistNeutralRectotET(0)  

							    ,fHistTotEMRectotET(0)

							    ,fHistMuonRecETDep(0) 
							    ,fHistMuonRec(0) 
							    ,fHistMuonMatchtotETDep(0) 
							    ,fHistMuonRecdEdxP(0)

							    ,fHistPionRecETDep(0) 
							    ,fHistPionRec(0) 
							    ,fHistPionMatchtotETDep(0) 
							    ,fHistPionRecdEdxP(0)

							    ,fHistKaonRecETDep(0) 
							    ,fHistKaonRec(0) 
							    ,fHistKaonMatchtotETDep(0) 
							    ,fHistKaonRecdEdxP(0)

							    ,fHistProtonRecETDep(0) 
							    ,fHistProtonRec(0) 
							    ,fHistProtonMatchtotETDep(0) 
							    ,fHistProtonRecdEdxP(0)

							    ,fHistTotChargedMatchtotETDep(0)

							    ,fHistTotalRectotETDep(0)

							    ,fHistDeltaRZ(0)
{//constructor
  fHistogramNameSuffix = TString("EmcalRec");
	
  fResCut = 0.02;
  //fResCut = fEmcalTrackDistanceCut;
	
  TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG));
  //TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
  TGeoManager::Import("geometry.root");
}

// dtor
AliAnalysisEmEtReconstructed::~AliAnalysisEmEtReconstructed() 
{//Destructor 
//   delete fGeoUt;

  delete fHistAllRecETDep;
  delete fHistAllRec;
  delete fHistAllRectotETDep;
	
  delete fHistElectronRecETDep;
  delete fHistElectronRec;
  delete fHistElectronMatchtotETDep; 
	
  delete fHistElectronRecdEdxP;

  delete fHistNeutralRectotET;  

  delete fHistTotEMRectotET;

  delete fHistMuonRecETDep;
  delete fHistMuonRec;
  delete fHistMuonMatchtotETDep; 

  delete fHistMuonRecdEdxP;
	
  delete fHistPionRecETDep;
  delete fHistPionRec;
  delete fHistPionMatchtotETDep; 

  delete fHistPionRecdEdxP;

  delete fHistKaonRecETDep;
  delete fHistKaonRec;
  delete fHistKaonMatchtotETDep; 

  delete fHistKaonRecdEdxP;
	
  delete fHistProtonRecETDep;
  delete fHistProtonRec;
  delete fHistProtonMatchtotETDep; 

  delete fHistProtonRecdEdxP;
	
  delete fHistTotChargedMatchtotETDep;
	
  delete fHistTotalRectotETDep;
	
  //few checks
  delete fHistDeltaRZ;
	
}

Int_t AliAnalysisEmEtReconstructed::AnalyseEvent(AliVEvent* ev)
{ // analyse MC and real event info
  if(!ev){
    AliError("ERROR: Event does not exist");   
    return 0;
  }
	
  fESD = dynamic_cast<AliESDEvent*>(ev);
  if (!fESD) {
    AliError("ERROR: Could not retrieve event");
    return 0;
  }

//   if(!fGeoUt){
//     fGeoUt = AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEARV1");//new AliEMCALGeometry("EMCAL_FIRSTYEAR","EMCAL");
//     //fGeoUt = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
//     AliInfo("Creating new AliEMCALGeometry");
//   }
//   //fGeoUt = new AliEMCALGeometry("EMCAL_COMPLETE1","EMCAL");
//   if(!fGeoUt){
//     AliInfo("No fGeoUt!");
//   }
//   else{
//     if(!fESD->GetEMCALMatrix(0)){
//       AliInfo("No matrix!");
//     }
//     else{
//       fGeoUt->SetMisalMatrix(fESD->GetEMCALMatrix(0),0);
//     }
//   }

  fCentBin= -1;
  if(fDataSet==20100){//If this is Pb+Pb
    AliCentrality *centrality = ev->GetCentrality();
    if(fNCentBins<21) fCentBin= centrality->GetCentralityClass10(fCentralityMethod);
    else{ fCentBin= centrality->GetCentralityClass5(fCentralityMethod);}
  }
	
  ResetEventValues();
	
  // get all emcal clusters
  TRefArray* caloClusters = new TRefArray();
  fESD->GetEMCALClusters( caloClusters );
	
  Int_t nCluster = caloClusters->GetEntries();
	
  Float_t pos[3] = {0};
  TVector3 caloPos(0,0,0);
  TVector3 trackPos(0,0,0);
//   Double_t res=0, delta_eta=0, delta_phi=0, maxPid=-99;
  Double_t res=0, maxPid=-99;
  Double_t xCluster[4]={0}, xCharged[7]={0};
	
  AliESDtrack *track = 0;
    
  // loop the clusters
  for (int iCluster = 0; iCluster < nCluster; iCluster++ ) 
    {		
      // Retrieve calo cluster information
      AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
      Float_t caloE = caloCluster->E();
      caloCluster->GetPosition(pos);		
      caloPos.SetXYZ(pos[0],pos[1],pos[2]);
		
      // look for track that matches calo cluster  
      //track = FindMatch(caloCluster, res); // Marcelo's matching	

      // *********************
      // tender's matching
//       delta_eta = caloCluster->GetTrackDz(); 
//       delta_phi = caloCluster->GetTrackDx(); 

      if (caloCluster->GetTrackMatchedIndex() > 0) // tender's matching
          track = fESD->GetTrack(caloCluster->GetTrackMatchedIndex());
        
      //if (track)
      //   if ( !fEsdtrackCutsITSTPC->IsSelected(track) )
      //       track = 0;
      // *********************
        
      // Retrieve track PID
      if (track) {
	maxPid = GetTrackPID(track);
	xCharged[0] = track->Eta();
	xCharged[1] = track->Pt();      }
      else {
	maxPid = -99;
	xCharged[0] = -99;
	xCharged[1] = -99;
      }

      // calculate ET
      Double_t etDep = CorrectForReconstructionEfficiency(*caloCluster);
		
      // All clusters
      //fHistAllRecEtaEDepETDep->Fill(caloE,caloPos.Eta(),etDep);
      //fHistAllRecEtaETDep->Fill(etDep,caloPos.Eta());
		
      xCluster[0] = caloE;
      xCluster[1] = caloPos.Eta();
      xCluster[2] = TMath::RadToDeg()*caloPos.Phi();
      xCluster[3] = caloCluster->GetNCells();
      fAllRectotETDep += etDep;		
		
	
      if(fMakeSparse){
	fHistAllRecETDep->Fill(xCluster,etDep);
	fHistAllRec->Fill(xCluster);
      }

      xCharged[2] = caloE;
      xCharged[3] = caloPos.Eta();
      xCharged[4] = TMath::RadToDeg()*caloPos.Phi();
      xCharged[5] = caloCluster->GetNCells();
      xCharged[6] = res;
		
      Bool_t isCharged = kFALSE;
		
      if (maxPid == AliPID::kProton)
	{
			
	  if(fMakeSparse){
	    fHistProtonRecETDep->Fill(xCharged,etDep);
	    fHistProtonRec->Fill(xCharged);
	  }

	  if (track) fHistProtonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
			
	  if ((res>0.) && (res<fResCut))
	    {
	      fProtonMatchtotETDep += etDep;

	      isCharged = kTRUE;
	    }
	}
      else if (maxPid == AliPID::kPion)
	{
			
	  if(fMakeSparse){
	    fHistPionRecETDep->Fill(xCharged,etDep);
	    fHistPionRec->Fill(xCharged);
	  }

	  if (track) fHistPionRecdEdxP->Fill(track->P(),track->GetTPCsignal());
			
	  if ((res>0.) && (res<fResCut))
	    {
	      fPionMatchtotETDep += etDep;
	      isCharged = kTRUE;
	    }
	}
      else if (maxPid == AliPID::kKaon)
	{
			
	  if(fMakeSparse){
	    fHistKaonRecETDep->Fill(xCharged,etDep);
	    fHistKaonRec->Fill(xCharged);
	  }

	  if (track) fHistKaonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
			
	  if ((res>0.) && (res<fResCut))
	    {

	      fKaonMatchtotETDep += etDep;
	      isCharged = kTRUE;
	    }
	}
      else if (maxPid == AliPID::kMuon)
	{
			
	  if(fMakeSparse){
	    fHistMuonRecETDep->Fill(xCharged,etDep);
	    fHistMuonRec->Fill(xCharged);
	  }
	  if (track) fHistMuonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
			
	  if ((res>0.) && (res<fResCut))
	    {

	      fMuonMatchtotETDep += etDep;						
	      isCharged = kTRUE;
	    }
	}
      else if (maxPid == AliPID::kElectron)
	{
			
	  if(fMakeSparse){
	    fHistElectronRecETDep->Fill(xCharged,etDep);
	    fHistElectronRec->Fill(xCharged);
	  }

	  if (track) fHistElectronRecdEdxP->Fill(track->P(),track->GetTPCsignal());
			
	  if ((res>0.) && (res<fResCut))
	    {
	      fElectronMatchtotETDep += etDep;
	      isCharged = kTRUE;
	    }
	}
		
      if (!isCharged)
	{
	  fNeutralRectotET += etDep;			
	}
		
    } // end of loop over clusters	
	
  fTotEMRectotET = fElectronMatchtotETDep + fNeutralRectotET;
  fTotChargedMatchtotETDep = fMuonMatchtotETDep + fPionMatchtotETDep + fKaonMatchtotETDep + fProtonMatchtotETDep;
  fTotalRectotETDep = fTotEMRectotET + fTotChargedMatchtotETDep;
	
  fHistAllRectotETDep->Fill(fAllRectotETDep);
	
  fHistElectronMatchtotETDep->Fill(fElectronMatchtotETDep); 
  fHistNeutralRectotET->Fill(fNeutralRectotET);
	
  fHistTotEMRectotET->Fill(fTotEMRectotET);
	
  fHistMuonMatchtotETDep->Fill(fMuonMatchtotETDep); 
  fHistPionMatchtotETDep->Fill(fPionMatchtotETDep); 
  fHistKaonMatchtotETDep->Fill(fKaonMatchtotETDep); 
  fHistProtonMatchtotETDep->Fill(fProtonMatchtotETDep); 
  fHistTotChargedMatchtotETDep->Fill(fTotChargedMatchtotETDep);
	
  fHistTotalRectotETDep->Fill(fTotalRectotETDep);
	
  delete caloClusters;
	
  return 0;    
}

void AliAnalysisEmEtReconstructed::Init()
{ // init
  AliAnalysisEt::Init();
}


void AliAnalysisEmEtReconstructed::ResetEventValues()
{ // reset event values
  AliAnalysisEt::ResetEventValues();
	
  // collision geometry defaults for p+p:
  fAllRectotETDep = 0;
	
  fElectronMatchtotETDep = 0;
  fNeutralRectotET = 0;
	
  fTotEMRectotET = 0;
	
  fMuonMatchtotETDep = 0; fPionMatchtotETDep = 0; fKaonMatchtotETDep = 0; fProtonMatchtotETDep = 0;
  fTotChargedMatchtotETDep = 0;
	
  fTotalRectotETDep = 0;
}


void AliAnalysisEmEtReconstructed::CreateHistograms()
{ // histogram related additions
  //AliAnalysisEt::CreateHistograms();
	
  if(fMakeSparse){
    fHistAllRecETDep = CreateClusterHistoSparse("fHistAllRecETDep_","E_{T}, all particles");
    fHistAllRec = CreateClusterHistoSparse("fHistAllRec_","counts, all particles");
  }
  TString histname = "fHistAllRectotETDep_" + fHistogramNameSuffix;
  fHistAllRectotETDep = new TH1F(histname.Data(),"total ET, all particles",fgNumOfEBins, fgEAxis);


  if(fMakeSparse){
    fHistElectronRecETDep = CreateChargedPartHistoSparse("fHistElectronRecETDep_","E_{T}, electrons");
    fHistElectronRec = CreateChargedPartHistoSparse("fHistElectronRec_","counts, electrons");
  }
  histname = "fHistElectronMatchtotETDep_" + fHistogramNameSuffix;
  fHistElectronMatchtotETDep = new TH1F(histname.Data(),"total ET, MC primary Electrons",fgNumOfEBins, fgEAxis);

  histname = "fHistElectronRecdEdxP_" + fHistogramNameSuffix;
  fHistElectronRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);

  histname = "fHistNeutralRectotET_" + fHistogramNameSuffix;
  fHistNeutralRectotET = new TH1F(histname.Data(),"total ET, neutral particles",fgNumOfEBins, fgEAxis);
	
  histname = "fHistTotEMRectotET_" + fHistogramNameSuffix;
  fHistTotEMRectotET = new TH1F(histname.Data(),"total electromagnetic ET",fgNumOfEBins, fgEAxis);

  if(fMakeSparse){
    fHistMuonRecETDep = CreateChargedPartHistoSparse("fHistMuonRecETDep_","E_{T}, muons");
    fHistMuonRec = CreateChargedPartHistoSparse("fHistMuonRec_","counts, muons");
  }
  histname = "fHistMuonMatchtotETDep_" + fHistogramNameSuffix;
  fHistMuonMatchtotETDep = new TH1F(histname.Data(),"total ET, Muons",fgNumOfEBins, fgEAxis);

  histname = "fHistMuonRecdEdxP_" + fHistogramNameSuffix;
  fHistMuonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
	
  if(fMakeSparse){
    fHistPionRecETDep = CreateChargedPartHistoSparse("fHistPionRecETDep_","E_{T}, pions");
    fHistPionRec = CreateChargedPartHistoSparse("fHistPionRec_","counts, pions");
  }
  histname = "fHistPionMatchtotETDep_" + fHistogramNameSuffix;
  fHistPionMatchtotETDep = new TH1F(histname.Data(),"total ET, Pions",fgNumOfEBins, fgEAxis);
  histname = "fHistPionRecdEdxP_" + fHistogramNameSuffix;
  fHistPionRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
	
  if(fMakeSparse){
    fHistKaonRecETDep = CreateChargedPartHistoSparse("fHistKaonRecETDep_","E_{T}, kaons");
    fHistKaonRec = CreateChargedPartHistoSparse("fHistKaonRec_","counts, kaons");
  }
  histname = "fHistKaonMatchtotETDep_" + fHistogramNameSuffix;
  fHistKaonMatchtotETDep = new TH1F(histname.Data(),"total ET, Kaons",fgNumOfEBins, fgEAxis);

  histname = "fHistKaonRecdEdxP_" + fHistogramNameSuffix;
  fHistKaonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
	
  if(fMakeSparse){
    fHistProtonRecETDep = CreateChargedPartHistoSparse("fHistProtonRecETDep_","E_{T}, protons");
    fHistProtonRec = CreateChargedPartHistoSparse("fHistProtonRec_","counts, protons");
  }
  histname = "fHistProtonMatchtotETDep_" + fHistogramNameSuffix;
  fHistProtonMatchtotETDep = new TH1F(histname.Data(),"total ET, Protons",fgNumOfEBins, fgEAxis);

  histname = "fHistProtonRecdEdxP_" + fHistogramNameSuffix;
  fHistProtonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
	
  histname = "fHistTotChargedMatchtotETDep_" + fHistogramNameSuffix;
  fHistTotChargedMatchtotETDep = new TH1F(histname.Data(),"total ET, charged particles",fgNumOfEBins, fgEAxis);
	
  histname = "fHistTotalRectotETDep_" + fHistogramNameSuffix;
  fHistTotalRectotETDep = new TH1F(histname.Data(),"total ET, all particles",fgNumOfEBins, fgEAxis);		
	
  histname = "fHistDeltaRZ_" + fHistogramNameSuffix;
  fHistDeltaRZ = new TH2F(histname,"#Delta#phi vs #Delta#eta (track projection - cluster position)",200,-0.1,0.1,200,-0.1,0.1);
}

void AliAnalysisEmEtReconstructed::FillOutputList(TList *list)
{//Function for filling the output list
  //AliAnalysisEt::FillOutputList(list);
	
	
  if(fMakeSparse){
    list->Add(fHistAllRecETDep); 
    list->Add(fHistAllRec); 
  }
  list->Add(fHistAllRectotETDep); 

  if(fMakeSparse){
    list->Add(fHistElectronRecETDep); 
    list->Add(fHistElectronRec); 
  }
  list->Add(fHistElectronMatchtotETDep); 
  list->Add(fHistElectronRecdEdxP);
	
	
  list->Add(fHistTotEMRectotET); 

  list->Add(fHistMuonRec); 
  list->Add(fHistMuonRecdEdxP);
  list->Add(fHistMuonMatchtotETDep); 

  if(fMakeSparse){
    list->Add(fHistPionRecETDep); 
    list->Add(fHistPionRec); 
  }
  list->Add(fHistPionMatchtotETDep); 
  list->Add(fHistPionRecdEdxP);
	
  if(fMakeSparse){
    list->Add(fHistKaonRecETDep); 
    list->Add(fHistKaonRec); 
  }
  list->Add(fHistKaonMatchtotETDep); 
  list->Add(fHistKaonRecdEdxP);

  if(fMakeSparse){
    list->Add(fHistProtonRecETDep); 
    list->Add(fHistProtonRec); 
  }
  list->Add(fHistProtonMatchtotETDep); 
  list->Add(fHistProtonRecdEdxP);
	
  list->Add(fHistTotChargedMatchtotETDep); 
  list->Add(fHistTotalRectotETDep); 
	
  list->Add(fHistDeltaRZ);
}

//________________________________________________________________________	
Double_t AliAnalysisEmEtReconstructed::GetTrackPID(const AliESDtrack *track) const
{//Get the default track ID
  const Double_t *pidWeights = track->PID();
  Int_t maxpid = -1;
  Double_t maxpidweight = 0;
	
  if (pidWeights)
    {
      for (Int_t p =0; p < AliPID::kSPECIES; p++)
	{
	  if (pidWeights[p] > maxpidweight)
	    {
	      maxpidweight = pidWeights[p];
	      maxpid = p;
	    }
	}
    }
	
  return maxpid;
}
Commit	Line	Data
2e12941b	1	//_________________________________________________________________________
	2	// Utility Class for transverse energy studies
	3	// Base class for MC analysis
	4	// - MC output
	5	// implementation file
	6	//
	7	//*-- Author: Marcelo G. Munhoz (USP)
f2adae2b	8	//_________________________________________________________________________
2e12941b	9
	10	#include "AliAnalysisEmEtReconstructed.h"
	11	#include "AliAnalysisEtCuts.h"
	12	#include "AliESDtrack.h"
	13	#include "AliStack.h"
	14	#include "AliVEvent.h"
	15	#include "AliMCEvent.h"
	16	#include "AliESDEvent.h"
	17	#include "TH2F.h"
	18	#include "TParticle.h"
	19	#include "AliGenHijingEventHeader.h"
	20	#include "AliGenPythiaEventHeader.h"
	21	#include "TList.h"
	22	#include "AliESDCaloCluster.h"
	23	#include "TGeoGlobalMagField.h"
	24	#include "AliMagF.h"
	25	#include "AliEMCALTrack.h"
	26	#include "AliESDtrackCuts.h"
	27	#include "AliEMCALGeometry.h"
	28	#include "AliExternalTrackParam.h"
	29	#include "AliTrackerBase.h"
	30	#include "TGeoManager.h"
2c70dc4d	31	#include "AliCentrality.h"
2e12941b	32
f2adae2b	33	using namespace std;
2e12941b	34
	35	ClassImp(AliAnalysisEmEtReconstructed);
	36
	37
	38	// ctor
	39	AliAnalysisEmEtReconstructed::AliAnalysisEmEtReconstructed():AliAnalysisEtReconstructed()
0af25ca1	40	,fResCut(0)
f2adae2b	41	,fAllRectotETDep(0)
	42	,fElectronMatchtotETDep(0)
	43	,fNeutralRectotET(0)
	44	,fTotEMRectotET(0)
	45
	46	,fMuonMatchtotETDep(0), fPionMatchtotETDep(0), fKaonMatchtotETDep(0), fProtonMatchtotETDep(0)
	47	,fTotChargedMatchtotETDep(0)
	48
	49	,fTotalRectotETDep(0)
	50
	51	,fESD(0)
2ed44673	52	// ,fGeoUt(0)
f2adae2b	53
	54	,fHistAllRecETDep(0)
	55	,fHistAllRec(0)
	56	,fHistAllRectotETDep(0)
	57
	58	,fHistElectronRecETDep(0)
	59	,fHistElectronRec(0)
	60	,fHistElectronMatchtotETDep(0)
	61	,fHistElectronRecdEdxP(0)
	62
	63	,fHistNeutralRectotET(0)
	64
	65	,fHistTotEMRectotET(0)
	66
	67	,fHistMuonRecETDep(0)
	68	,fHistMuonRec(0)
	69	,fHistMuonMatchtotETDep(0)
	70	,fHistMuonRecdEdxP(0)
	71
	72	,fHistPionRecETDep(0)
	73	,fHistPionRec(0)
	74	,fHistPionMatchtotETDep(0)
	75	,fHistPionRecdEdxP(0)
	76
	77	,fHistKaonRecETDep(0)
	78	,fHistKaonRec(0)
	79	,fHistKaonMatchtotETDep(0)
	80	,fHistKaonRecdEdxP(0)
	81
	82	,fHistProtonRecETDep(0)
	83	,fHistProtonRec(0)
	84	,fHistProtonMatchtotETDep(0)
	85	,fHistProtonRecdEdxP(0)
	86
	87	,fHistTotChargedMatchtotETDep(0)
	88
	89	,fHistTotalRectotETDep(0)
	90
	91	,fHistDeltaRZ(0)
2e12941b	92	{//constructor
f2adae2b	93	fHistogramNameSuffix = TString("EmcalRec");
2e12941b	94
f2adae2b	95	fResCut = 0.02;
f2adae2b	96	//fResCut = fEmcalTrackDistanceCut;
2e12941b	97
f2adae2b	98	TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG));
	99	//TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
	100	TGeoManager::Import("geometry.root");
2e12941b	101	}
	102
	103	// dtor
	104	AliAnalysisEmEtReconstructed::~AliAnalysisEmEtReconstructed()
f2adae2b	105	{//Destructor
2ed44673	106	// delete fGeoUt;
a4d91e7c	107
f2adae2b	108	delete fHistAllRecETDep;
	109	delete fHistAllRec;
	110	delete fHistAllRectotETDep;
a4d91e7c	111
f2adae2b	112	delete fHistElectronRecETDep;
	113	delete fHistElectronRec;
	114	delete fHistElectronMatchtotETDep;
a4d91e7c	115
f2adae2b	116	delete fHistElectronRecdEdxP;
a4d91e7c	117
f2adae2b	118	delete fHistNeutralRectotET;
a4d91e7c	119
f2adae2b	120	delete fHistTotEMRectotET;
a4d91e7c	121
f2adae2b	122	delete fHistMuonRecETDep;
	123	delete fHistMuonRec;
	124	delete fHistMuonMatchtotETDep;
a4d91e7c	125
f2adae2b	126	delete fHistMuonRecdEdxP;
a4d91e7c	127
f2adae2b	128	delete fHistPionRecETDep;
	129	delete fHistPionRec;
	130	delete fHistPionMatchtotETDep;
a4d91e7c	131
f2adae2b	132	delete fHistPionRecdEdxP;
a4d91e7c	133
f2adae2b	134	delete fHistKaonRecETDep;
	135	delete fHistKaonRec;
	136	delete fHistKaonMatchtotETDep;
a4d91e7c	137
f2adae2b	138	delete fHistKaonRecdEdxP;
a4d91e7c	139
f2adae2b	140	delete fHistProtonRecETDep;
	141	delete fHistProtonRec;
	142	delete fHistProtonMatchtotETDep;
a4d91e7c	143
f2adae2b	144	delete fHistProtonRecdEdxP;
a4d91e7c	145
f2adae2b	146	delete fHistTotChargedMatchtotETDep;
a4d91e7c	147
f2adae2b	148	delete fHistTotalRectotETDep;
a4d91e7c	149
f2adae2b	150	//few checks
f2adae2b	151	delete fHistDeltaRZ;
a4d91e7c	152
2e12941b	153	}
	154
	155	Int_t AliAnalysisEmEtReconstructed::AnalyseEvent(AliVEvent* ev)
	156	{ // analyse MC and real event info
f2adae2b	157	if(!ev){
	158	AliError("ERROR: Event does not exist");
	159	return 0;
	160	}
	161
	162	fESD = dynamic_cast<AliESDEvent*>(ev);
6481ee88	163	if (!fESD) {
	164	AliError("ERROR: Could not retrieve event");
	165	return 0;
	166	}
f2adae2b	167
2ed44673	168	// if(!fGeoUt){
	169	// fGeoUt = AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEARV1");//new AliEMCALGeometry("EMCAL_FIRSTYEAR","EMCAL");
	170	// //fGeoUt = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
	171	// AliInfo("Creating new AliEMCALGeometry");
	172	// }
	173	// //fGeoUt = new AliEMCALGeometry("EMCAL_COMPLETE1","EMCAL");
	174	// if(!fGeoUt){
	175	// AliInfo("No fGeoUt!");
	176	// }
	177	// else{
	178	// if(!fESD->GetEMCALMatrix(0)){
	179	// AliInfo("No matrix!");
	180	// }
	181	// else{
	182	// fGeoUt->SetMisalMatrix(fESD->GetEMCALMatrix(0),0);
	183	// }
	184	// }
2c70dc4d	185
	186	fCentBin= -1;
	187	if(fDataSet==20100){//If this is Pb+Pb
	188	AliCentrality *centrality = ev->GetCentrality();
	189	if(fNCentBins<21) fCentBin= centrality->GetCentralityClass10(fCentralityMethod);
	190	else{ fCentBin= centrality->GetCentralityClass5(fCentralityMethod);}
	191	}
2e12941b	192
f2adae2b	193	ResetEventValues();
2e12941b	194
f2adae2b	195	// get all emcal clusters
	196	TRefArray* caloClusters = new TRefArray();
	197	fESD->GetEMCALClusters( caloClusters );
2e12941b	198
f2adae2b	199	Int_t nCluster = caloClusters->GetEntries();
2e12941b	200
f2adae2b	201	Float_t pos[3] = {0};
	202	TVector3 caloPos(0,0,0);
	203	TVector3 trackPos(0,0,0);
f6b36c54	204	// Double_t res=0, delta_eta=0, delta_phi=0, maxPid=-99;
f6b36c54	205	Double_t res=0, maxPid=-99;
f2adae2b	206	Double_t xCluster[4]={0}, xCharged[7]={0};
7687e5b3	207
	208	AliESDtrack *track = 0;
	209
f2adae2b	210	// loop the clusters
	211	for (int iCluster = 0; iCluster < nCluster; iCluster++ )
	212	{
	213	// Retrieve calo cluster information
	214	AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
	215	Float_t caloE = caloCluster->E();
	216	caloCluster->GetPosition(pos);
	217	caloPos.SetXYZ(pos[0],pos[1],pos[2]);
2e12941b	218
f2adae2b	219	// look for track that matches calo cluster
7687e5b3	220	//track = FindMatch(caloCluster, res); // Marcelo's matching
	221
	222	// *********************
	223	// tender's matching
f6b36c54	224	// delta_eta = caloCluster->GetTrackDz();
f6b36c54	225	// delta_phi = caloCluster->GetTrackDx();
7687e5b3	226
	227	if (caloCluster->GetTrackMatchedIndex() > 0) // tender's matching
	228	track = fESD->GetTrack(caloCluster->GetTrackMatchedIndex());
	229
	230	//if (track)
	231	// if ( !fEsdtrackCutsITSTPC->IsSelected(track) )
	232	// track = 0;
	233	// *********************
	234
f2adae2b	235	// Retrieve track PID
6481ee88	236	if (track) {
f2adae2b	237	maxPid = GetTrackPID(track);
6481ee88	238	xCharged[0] = track->Eta();
	239	xCharged[1] = track->Pt(); }
	240	else {
f2adae2b	241	maxPid = -99;
6481ee88	242	xCharged[0] = -99;
	243	xCharged[1] = -99;
	244	}
	245
f2adae2b	246	// calculate ET
cfe23ff0	247	Double_t etDep = CorrectForReconstructionEfficiency(*caloCluster);
2e12941b	248
f2adae2b	249	// All clusters
	250	//fHistAllRecEtaEDepETDep->Fill(caloE,caloPos.Eta(),etDep);
	251	//fHistAllRecEtaETDep->Fill(etDep,caloPos.Eta());
2e12941b	252
f2adae2b	253	xCluster[0] = caloE;
	254	xCluster[1] = caloPos.Eta();
	255	xCluster[2] = TMath::RadToDeg()*caloPos.Phi();
	256	xCluster[3] = caloCluster->GetNCells();
	257	fAllRectotETDep += etDep;
2e12941b	258
f2adae2b	259
	260	if(fMakeSparse){
	261	fHistAllRecETDep->Fill(xCluster,etDep);
	262	fHistAllRec->Fill(xCluster);
	263	}
	264
f2adae2b	265	xCharged[2] = caloE;
	266	xCharged[3] = caloPos.Eta();
	267	xCharged[4] = TMath::RadToDeg()*caloPos.Phi();
	268	xCharged[5] = caloCluster->GetNCells();
	269	xCharged[6] = res;
2e12941b	270
f2adae2b	271	Bool_t isCharged = kFALSE;
2e12941b	272
f2adae2b	273	if (maxPid == AliPID::kProton)
f2adae2b	274	{
2e12941b	275
f2adae2b	276	if(fMakeSparse){
	277	fHistProtonRecETDep->Fill(xCharged,etDep);
	278	fHistProtonRec->Fill(xCharged);
	279	}
2e12941b	280
6481ee88	281	if (track) fHistProtonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
2e12941b	282
f2adae2b	283	if ((res>0.) && (res<fResCut))
	284	{
	285	fProtonMatchtotETDep += etDep;
	286
	287	isCharged = kTRUE;
	288	}
	289	}
	290	else if (maxPid == AliPID::kPion)
	291	{
2e12941b	292
f2adae2b	293	if(fMakeSparse){
	294	fHistPionRecETDep->Fill(xCharged,etDep);
	295	fHistPionRec->Fill(xCharged);
	296	}
2e12941b	297
6481ee88	298	if (track) fHistPionRecdEdxP->Fill(track->P(),track->GetTPCsignal());
2e12941b	299
f2adae2b	300	if ((res>0.) && (res<fResCut))
	301	{
	302	fPionMatchtotETDep += etDep;
	303	isCharged = kTRUE;
	304	}
	305	}
	306	else if (maxPid == AliPID::kKaon)
	307	{
2e12941b	308
f2adae2b	309	if(fMakeSparse){
	310	fHistKaonRecETDep->Fill(xCharged,etDep);
	311	fHistKaonRec->Fill(xCharged);
	312	}
2e12941b	313
6481ee88	314	if (track) fHistKaonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
2e12941b	315
f2adae2b	316	if ((res>0.) && (res<fResCut))
f2adae2b	317	{
2e12941b	318
f2adae2b	319	fKaonMatchtotETDep += etDep;
	320	isCharged = kTRUE;
	321	}
	322	}
	323	else if (maxPid == AliPID::kMuon)
	324	{
2e12941b	325
f2adae2b	326	if(fMakeSparse){
	327	fHistMuonRecETDep->Fill(xCharged,etDep);
	328	fHistMuonRec->Fill(xCharged);
	329	}
6481ee88	330	if (track) fHistMuonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
2e12941b	331
f2adae2b	332	if ((res>0.) && (res<fResCut))
f2adae2b	333	{
2e12941b	334
f2adae2b	335	fMuonMatchtotETDep += etDep;
	336	isCharged = kTRUE;
	337	}
	338	}
	339	else if (maxPid == AliPID::kElectron)
	340	{
2e12941b	341
f2adae2b	342	if(fMakeSparse){
	343	fHistElectronRecETDep->Fill(xCharged,etDep);
	344	fHistElectronRec->Fill(xCharged);
	345	}
	346
6481ee88	347	if (track) fHistElectronRecdEdxP->Fill(track->P(),track->GetTPCsignal());
f2adae2b	348
	349	if ((res>0.) && (res<fResCut))
	350	{
	351	fElectronMatchtotETDep += etDep;
	352	isCharged = kTRUE;
	353	}
	354	}
2e12941b	355
f2adae2b	356	if (!isCharged)
	357	{
	358	fNeutralRectotET += etDep;
	359	}
2e12941b	360
f2adae2b	361	} // end of loop over clusters
2e12941b	362
f2adae2b	363	fTotEMRectotET = fElectronMatchtotETDep + fNeutralRectotET;
	364	fTotChargedMatchtotETDep = fMuonMatchtotETDep + fPionMatchtotETDep + fKaonMatchtotETDep + fProtonMatchtotETDep;
	365	fTotalRectotETDep = fTotEMRectotET + fTotChargedMatchtotETDep;
2e12941b	366
f2adae2b	367	fHistAllRectotETDep->Fill(fAllRectotETDep);
2e12941b	368
f2adae2b	369	fHistElectronMatchtotETDep->Fill(fElectronMatchtotETDep);
f2adae2b	370	fHistNeutralRectotET->Fill(fNeutralRectotET);
2e12941b	371
f2adae2b	372	fHistTotEMRectotET->Fill(fTotEMRectotET);
2e12941b	373
f2adae2b	374	fHistMuonMatchtotETDep->Fill(fMuonMatchtotETDep);
	375	fHistPionMatchtotETDep->Fill(fPionMatchtotETDep);
	376	fHistKaonMatchtotETDep->Fill(fKaonMatchtotETDep);
	377	fHistProtonMatchtotETDep->Fill(fProtonMatchtotETDep);
	378	fHistTotChargedMatchtotETDep->Fill(fTotChargedMatchtotETDep);
2e12941b	379
f2adae2b	380	fHistTotalRectotETDep->Fill(fTotalRectotETDep);
2e12941b	381
f2adae2b	382	delete caloClusters;
2e12941b	383
f2adae2b	384	return 0;
2e12941b	385	}
	386
	387	void AliAnalysisEmEtReconstructed::Init()
	388	{ // init
f2adae2b	389	AliAnalysisEt::Init();
2e12941b	390	}
	391
	392
	393	void AliAnalysisEmEtReconstructed::ResetEventValues()
	394	{ // reset event values
f2adae2b	395	AliAnalysisEt::ResetEventValues();
2e12941b	396
f2adae2b	397	// collision geometry defaults for p+p:
f2adae2b	398	fAllRectotETDep = 0;
2e12941b	399
f2adae2b	400	fElectronMatchtotETDep = 0;
f2adae2b	401	fNeutralRectotET = 0;
2e12941b	402
f2adae2b	403	fTotEMRectotET = 0;
2e12941b	404
f2adae2b	405	fMuonMatchtotETDep = 0; fPionMatchtotETDep = 0; fKaonMatchtotETDep = 0; fProtonMatchtotETDep = 0;
f2adae2b	406	fTotChargedMatchtotETDep = 0;
2e12941b	407
f2adae2b	408	fTotalRectotETDep = 0;
2e12941b	409	}
	410
	411
	412	void AliAnalysisEmEtReconstructed::CreateHistograms()
	413	{ // histogram related additions
f2adae2b	414	//AliAnalysisEt::CreateHistograms();
	415
	416	if(fMakeSparse){
	417	fHistAllRecETDep = CreateClusterHistoSparse("fHistAllRecETDep_","E_{T}, all particles");
	418	fHistAllRec = CreateClusterHistoSparse("fHistAllRec_","counts, all particles");
	419	}
	420	TString histname = "fHistAllRectotETDep_" + fHistogramNameSuffix;
	421	fHistAllRectotETDep = new TH1F(histname.Data(),"total ET, all particles",fgNumOfEBins, fgEAxis);
	422
	423
	424	if(fMakeSparse){
	425	fHistElectronRecETDep = CreateChargedPartHistoSparse("fHistElectronRecETDep_","E_{T}, electrons");
	426	fHistElectronRec = CreateChargedPartHistoSparse("fHistElectronRec_","counts, electrons");
	427	}
	428	histname = "fHistElectronMatchtotETDep_" + fHistogramNameSuffix;
	429	fHistElectronMatchtotETDep = new TH1F(histname.Data(),"total ET, MC primary Electrons",fgNumOfEBins, fgEAxis);
	430
	431	histname = "fHistElectronRecdEdxP_" + fHistogramNameSuffix;
	432	fHistElectronRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
	433
	434	histname = "fHistNeutralRectotET_" + fHistogramNameSuffix;
	435	fHistNeutralRectotET = new TH1F(histname.Data(),"total ET, neutral particles",fgNumOfEBins, fgEAxis);
	436
	437	histname = "fHistTotEMRectotET_" + fHistogramNameSuffix;
	438	fHistTotEMRectotET = new TH1F(histname.Data(),"total electromagnetic ET",fgNumOfEBins, fgEAxis);
	439
	440	if(fMakeSparse){
	441	fHistMuonRecETDep = CreateChargedPartHistoSparse("fHistMuonRecETDep_","E_{T}, muons");
	442	fHistMuonRec = CreateChargedPartHistoSparse("fHistMuonRec_","counts, muons");
	443	}
	444	histname = "fHistMuonMatchtotETDep_" + fHistogramNameSuffix;
	445	fHistMuonMatchtotETDep = new TH1F(histname.Data(),"total ET, Muons",fgNumOfEBins, fgEAxis);
	446
	447	histname = "fHistMuonRecdEdxP_" + fHistogramNameSuffix;
	448	fHistMuonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
	449
	450	if(fMakeSparse){
	451	fHistPionRecETDep = CreateChargedPartHistoSparse("fHistPionRecETDep_","E_{T}, pions");
	452	fHistPionRec = CreateChargedPartHistoSparse("fHistPionRec_","counts, pions");
	453	}
	454	histname = "fHistPionMatchtotETDep_" + fHistogramNameSuffix;
	455	fHistPionMatchtotETDep = new TH1F(histname.Data(),"total ET, Pions",fgNumOfEBins, fgEAxis);
	456	histname = "fHistPionRecdEdxP_" + fHistogramNameSuffix;
	457	fHistPionRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
	458
	459	if(fMakeSparse){
	460	fHistKaonRecETDep = CreateChargedPartHistoSparse("fHistKaonRecETDep_","E_{T}, kaons");
	461	fHistKaonRec = CreateChargedPartHistoSparse("fHistKaonRec_","counts, kaons");
	462	}
	463	histname = "fHistKaonMatchtotETDep_" + fHistogramNameSuffix;
	464	fHistKaonMatchtotETDep = new TH1F(histname.Data(),"total ET, Kaons",fgNumOfEBins, fgEAxis);
	465
	466	histname = "fHistKaonRecdEdxP_" + fHistogramNameSuffix;
	467	fHistKaonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
	468
	469	if(fMakeSparse){
	470	fHistProtonRecETDep = CreateChargedPartHistoSparse("fHistProtonRecETDep_","E_{T}, protons");
	471	fHistProtonRec = CreateChargedPartHistoSparse("fHistProtonRec_","counts, protons");
	472	}
	473	histname = "fHistProtonMatchtotETDep_" + fHistogramNameSuffix;
	474	fHistProtonMatchtotETDep = new TH1F(histname.Data(),"total ET, Protons",fgNumOfEBins, fgEAxis);
	475
	476	histname = "fHistProtonRecdEdxP_" + fHistogramNameSuffix;
	477	fHistProtonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
478
479	histname = "fHistTotChargedMatchtotETDep_" + fHistogramNameSuffix;
480	fHistTotChargedMatchtotETDep = new TH1F(histname.Data(),"total ET, charged particles",fgNumOfEBins, fgEAxis);
481
482	histname = "fHistTotalRectotETDep_" + fHistogramNameSuffix;
483	fHistTotalRectotETDep = new TH1F(histname.Data(),"total ET, all particles",fgNumOfEBins, fgEAxis);
484
485	histname = "fHistDeltaRZ_" + fHistogramNameSuffix;
486	fHistDeltaRZ = new TH2F(histname,"#Delta#phi vs #Delta#eta (track projection - cluster position)",200,-0.1,0.1,200,-0.1,0.1);
2e12941b	487	}
	488
	489	void AliAnalysisEmEtReconstructed::FillOutputList(TList *list)
	490	{//Function for filling the output list
f2adae2b	491	//AliAnalysisEt::FillOutputList(list);
	492
	493
	494	if(fMakeSparse){
	495	list->Add(fHistAllRecETDep);
	496	list->Add(fHistAllRec);
	497	}
	498	list->Add(fHistAllRectotETDep);
	499
	500	if(fMakeSparse){
	501	list->Add(fHistElectronRecETDep);
	502	list->Add(fHistElectronRec);
	503	}
	504	list->Add(fHistElectronMatchtotETDep);
	505	list->Add(fHistElectronRecdEdxP);
	506
	507
	508	list->Add(fHistTotEMRectotET);
	509
	510	list->Add(fHistMuonRec);
	511	list->Add(fHistMuonRecdEdxP);
	512	list->Add(fHistMuonMatchtotETDep);
	513
	514	if(fMakeSparse){
	515	list->Add(fHistPionRecETDep);
	516	list->Add(fHistPionRec);
	517	}
	518	list->Add(fHistPionMatchtotETDep);
	519	list->Add(fHistPionRecdEdxP);
	520
	521	if(fMakeSparse){
	522	list->Add(fHistKaonRecETDep);
	523	list->Add(fHistKaonRec);
	524	}
	525	list->Add(fHistKaonMatchtotETDep);
	526	list->Add(fHistKaonRecdEdxP);
	527
	528	if(fMakeSparse){
	529	list->Add(fHistProtonRecETDep);
	530	list->Add(fHistProtonRec);
	531	}
	532	list->Add(fHistProtonMatchtotETDep);
	533	list->Add(fHistProtonRecdEdxP);
	534
	535	list->Add(fHistTotChargedMatchtotETDep);
	536	list->Add(fHistTotalRectotETDep);
	537
	538	list->Add(fHistDeltaRZ);
2e12941b	539	}
2e12941b	540
2e12941b	541	//________________________________________________________________________
	542	Double_t AliAnalysisEmEtReconstructed::GetTrackPID(const AliESDtrack *track) const
	543	{//Get the default track ID
f2adae2b	544	const Double_t *pidWeights = track->PID();
	545	Int_t maxpid = -1;
	546	Double_t maxpidweight = 0;
2e12941b	547
f2adae2b	548	if (pidWeights)
	549	{
	550	for (Int_t p =0; p < AliPID::kSPECIES; p++)
2e12941b	551	{
f2adae2b	552	if (pidWeights[p] > maxpidweight)
	553	{
	554	maxpidweight = pidWeights[p];
	555	maxpid = p;
	556	}
2e12941b	557	}
f2adae2b	558	}
2e12941b	559
f2adae2b	560	return maxpid;
2e12941b	561	}