[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / AliCaloTrackAODReader.cxx


/**************************************************************************
 * Copyright(c) 1998-1999, 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.                  *
 **************************************************************************/
/* $Id: $ */

//_________________________________________________________________________
// Class for reading data (AODs) in order to do prompt gamma
//  or other particle identification and correlations.
// Mixing analysis can be done, input AOD with events
// is opened in the AliCaloTrackReader::Init()
// 
//
//*-- Author: Gustavo Conesa (LNF-INFN) 
//////////////////////////////////////////////////////////////////////////////


// --- ROOT system ---
//#include "Riostream.h"

//---- ANALYSIS system ----
#include "AliCaloTrackAODReader.h" 
#include "AliAODCaloCluster.h"
#include "AliAODTrack.h"
#include "AliESDtrack.h"
#include "AliFiducialCut.h"
#include "AliAODInputHandler.h"
#include "AliAnalysisManager.h"

ClassImp(AliCaloTrackAODReader)

//____________________________________________________________________________
AliCaloTrackAODReader::AliCaloTrackAODReader() : 
  AliCaloTrackReader(), fWriteOutputStdAOD(kFALSE)
{
  //Default Ctor
  
  //Initialize parameters
  fDataType=kAOD;
  fReadStack          = kTRUE;
  fReadAODMCParticles = kFALSE;
  //We want tracks fitted in the detectors:
  fTrackStatus=AliESDtrack::kTPCrefit;
  fTrackStatus|=AliESDtrack::kITSrefit;

}

//____________________________________________________________________________
AliCaloTrackAODReader::AliCaloTrackAODReader(const AliCaloTrackAODReader & aodr) :   
  AliCaloTrackReader(aodr), fWriteOutputStdAOD(aodr.fWriteOutputStdAOD)
{
  // cpy ctor
}

//_________________________________________________________________________
//AliCaloTrackAODReader & AliCaloTrackAODReader::operator = (const AliCaloTrackAODReader & source)
//{
//  // assignment operator
//
//  if(&source == this) return *this;
//
//  return *this;
//
//}

//____________________________________________________________________________
void AliCaloTrackAODReader::FillInputCTS() {
  //Return array with Central Tracking System (CTS) tracks
  if(fDebug > 2 ) printf("AliCaloTrackAODReader::FillInputCTS()\n");
  Int_t nTracks   = fInputEvent->GetNumberOfTracks() ;
  Int_t naod = 0;
  Double_t p[3];
  for (Int_t itrack =  0; itrack <  nTracks; itrack++) {////////////// track loop
    AliAODTrack * track = ((AliAODEvent*)fInputEvent)->GetTrack(itrack) ; // retrieve track from esd
    	  
    //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
	if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus)) continue ;
    
    track->GetPxPyPz(p) ;
    TLorentzVector momentum(p[0],p[1],p[2],0);
    
    if(fCTSPtMin < momentum.Pt() && fFiducialCut->IsInFiducialCut(momentum,"CTS")){
      
      if(fDebug > 2 && momentum.Pt() > 0.1) printf("AliCaloTrackAODReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
						  momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
	  
	  if(fWriteOutputStdAOD){
		   AliAODTrack* newtrack = new((*(fOutputEvent->GetTracks()))[naod++]) AliAODTrack(*track);
		   fAODCTS->Add(newtrack); //Use AOD stored in output for references.
	  }
	  else fAODCTS->Add(track);
    }//Pt and Fiducial cut passed. 
  }// track loop
	
  fAODCTSNormalInputEntries = fAODCTS->GetEntriesFast();
  if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputCTS()   - aod entries %d\n", fAODCTSNormalInputEntries);

  //If second input event available, add the clusters.
  if(fSecondInputAODTree && fSecondInputAODEvent){
	  nTracks   = fSecondInputAODEvent->GetNumberOfTracks() ;
	  if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputCTS()   - Add second input tracks, entries %d\n", nTracks) ;
	  for (Int_t itrack =  0; itrack <  nTracks; itrack++) {////////////// track loop
		  AliAODTrack * track = ((AliAODEvent*)fSecondInputAODEvent)->GetTrack(itrack) ; // retrieve track from esd
		  
		  //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
		  if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus)) continue ;
		  
		  track->GetPxPyPz(p) ;
		  TLorentzVector momentum(p[0],p[1],p[2],0);
		  
		  if(fCTSPtMin < momentum.Pt() && fFiducialCut->IsInFiducialCut(momentum,"CTS")){
			  
			  if(fDebug > 2 && momentum.Pt() > 0.1) printf("AliCaloTrackAODReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
														   momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
			  
			  if(fWriteOutputStdAOD){
				  AliAODTrack* newtrack = new((*(fOutputEvent->GetTracks()))[naod++]) AliAODTrack(*track);
				  fAODCTS->Add(newtrack); //Use AOD stored in output for references.
			  }
			  else fAODCTS->Add(track);
			  
		  }//Pt and Fiducial cut passed. 
	  }// track loop
	  
	  if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputCTS()   - aod normal entries %d, after second input %d\n", fAODCTSNormalInputEntries, fAODCTS->GetEntriesFast());
  }	//second input loop
	
}

//____________________________________________________________________________
void AliCaloTrackAODReader::FillInputEMCAL() {
  //Return array with EMCAL clusters in aod format
  if(fDebug > 2 ) printf("AliCaloTrackAODReader::FillInputEMCAL()\n");
	
  //Get vertex for momentum calculation  
  Double_t v[3] ; //vertex ;
  GetVertex(v);

  Int_t naod =  (fOutputEvent->GetCaloClusters())->GetEntriesFast();
  //Loop to select clusters in fiducial cut and fill container with aodClusters
  Int_t nclusters = ((AliAODEvent*)fInputEvent)->GetNCaloClusters();
  for (Int_t iclus =  0; iclus <  nclusters; iclus++) {
    AliAODCaloCluster * clus = 0;
    if ( (clus = ((AliAODEvent*)fInputEvent)->GetCaloCluster(iclus)) ) {
      if (clus->IsEMCALCluster()){
		  
	//Check if the cluster contains any bad channel
	if(ClusterContainsBadChannel("EMCAL",clus->GetCellsAbsId(), clus->GetNCells())) continue;
		  
	TLorentzVector momentum ;
	clus->GetMomentum(momentum, v);      
	
	if(fEMCALPtMin < momentum.Pt() && fFiducialCut->IsInFiducialCut(momentum,"EMCAL")){
    
	  if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackAODReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
						     momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
	  
	  if(fWriteOutputStdAOD){
		AliAODCaloCluster * newclus = new((*(fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(*clus);
		fAODEMCAL->Add(newclus);	
	  }
	  else fAODEMCAL->Add(clus);	
	}//Pt and Fiducial cut passed.
      }//EMCAL cluster
    }// cluster exists
  }// cluster loop
	
  fAODEMCALNormalInputEntries = fAODEMCAL->GetEntriesFast();
  if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputEMCAL() - aod entries %d\n", fAODEMCALNormalInputEntries);

  //If second input event available, add the clusters.
  if(fSecondInputAODTree && fSecondInputAODEvent){
	  GetSecondInputAODVertex(v);
	  nclusters = ((AliAODEvent*)fSecondInputAODEvent)->GetNCaloClusters();
	  if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputEMCAL() - Add second input clusters, entries %d\n", nclusters) ;
		for (Int_t iclus =  0; iclus < nclusters; iclus++) {
			AliAODCaloCluster * clus = 0;
			if ( (clus = ((AliAODEvent*)fSecondInputAODEvent)->GetCaloCluster(iclus)) ) {
				if (clus->IsEMCALCluster()){
					TLorentzVector momentum ;
					clus->GetMomentum(momentum, v);      
					
					if(fEMCALPtMin < momentum.Pt() && fFiducialCut->IsInFiducialCut(momentum,"EMCAL")){
						
						if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackAODReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
																	momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
						if(fWriteOutputStdAOD){
							AliAODCaloCluster * newclus = new((*(fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(*clus);
						    fAODEMCAL->Add(newclus);	
						}
						else fAODEMCAL->Add(clus);	
					}//Pt and Fiducial cut passed.
				}//EMCAL cluster
			}// cluster exists
		}// cluster loop
		
	  if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputEMCAL() - aod normal entries %d, after second input %d\n", fAODEMCALNormalInputEntries, fAODEMCAL->GetEntriesFast());

	} //second input loop
}

//____________________________________________________________________________
void AliCaloTrackAODReader::FillInputPHOS() {
  //Return array with PHOS clusters in aod format
  if(fDebug > 2 ) printf("AliCaloTrackAODReader::FillInputPHOS()\n");
	
  //Get vertex for momentum calculation  
  Double_t v[3] ; //vertex ;
  GetVertex(v);
	
  Int_t naod =  (fOutputEvent->GetCaloClusters())->GetEntriesFast();
  //Loop to select clusters in fiducial cut and fill container with aodClusters
  Int_t nclusters = ((AliAODEvent*)fInputEvent)->GetNCaloClusters();
  for (Int_t iclus =  0; iclus < nclusters; iclus++) {
    AliAODCaloCluster * clus = 0;
    if ( (clus = ((AliAODEvent*)fInputEvent)->GetCaloCluster(iclus)) ) {
      if (clus->IsPHOSCluster()){
		  
	//Check if the cluster contains any bad channel
	if(ClusterContainsBadChannel("PHOS",clus->GetCellsAbsId(), clus->GetNCells())) continue;
		  
	TLorentzVector momentum ;
	clus->GetMomentum(momentum, v);      
	
	if(fPHOSPtMin < momentum.Pt() && fFiducialCut->IsInFiducialCut(momentum,"PHOS")){
	  
	  if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackAODReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
						     momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());

		if(fWriteOutputStdAOD){
			AliAODCaloCluster * newclus = new((*(fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(*clus);
			fAODPHOS->Add(newclus);	
		}
		else fAODPHOS->Add(clus);	
	}//Pt and Fiducial cut passed.
      }//PHOS cluster
    }//cluster exists
  }//esd cluster loop
	
  fAODPHOSNormalInputEntries = fAODPHOS->GetEntriesFast() ;
  if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputPHOS()  - aod entries %d\n", fAODPHOSNormalInputEntries);

  //If second input event available, add the clusters.
  if(fSecondInputAODTree && fSecondInputAODEvent){  
	  GetSecondInputAODVertex(v);
	  nclusters = ((AliAODEvent*)fSecondInputAODEvent)->GetNCaloClusters();
	  if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputPHOS()  - Add second input clusters, entries %d\n", nclusters);
		for (Int_t iclus =  0; iclus < nclusters; iclus++) {
			AliAODCaloCluster * clus = 0;
			if ( (clus = ((AliAODEvent*)fSecondInputAODEvent)->GetCaloCluster(iclus)) ) {
				if (clus->IsPHOSCluster()){
					TLorentzVector momentum ;
					clus->GetMomentum(momentum, v);      
					
					if(fPHOSPtMin < momentum.Pt() && fFiducialCut->IsInFiducialCut(momentum,"PHOS")){
						
						if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackAODReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
																	momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
						if(fWriteOutputStdAOD){
							AliAODCaloCluster * newclus = new((*(fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(*clus);
						    fAODPHOS->Add(newclus);	
						}
						else fAODPHOS->Add(clus);	
					}//Pt and Fiducial cut passed.
				}//PHOS cluster
			}// cluster exists
		}// cluster loop
		if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputPHOS()  - aod normal entries %d, after second input %d\n", fAODPHOSNormalInputEntries, fAODPHOS->GetEntriesFast());
  }	//second input loop

}

//____________________________________________________________________________
void AliCaloTrackAODReader::FillInputEMCALCells() {
  //Return array with EMCAL cells in aod format

  fEMCALCells = (TNamed*) ((AliAODEvent*)fInputEvent)->GetEMCALCells(); 

}

//____________________________________________________________________________
void AliCaloTrackAODReader::FillInputPHOSCells() {
  //Return array with PHOS cells in aod format

  fPHOSCells = (TNamed*) ((AliAODEvent*)fInputEvent)->GetPHOSCells(); 

}

//____________________________________________________________________________
void AliCaloTrackAODReader::GetVertex(Double_t  v[3]) const {
  //Return vertex position

 ((AliAODEvent*)fInputEvent)->GetPrimaryVertex()->GetXYZ(v);
	
}

//____________________________________________________________________________
void AliCaloTrackAODReader::GetSecondInputAODVertex(Double_t  v[3]) const {
	//Return vertex position of second AOD input
	
	fSecondInputAODEvent->GetPrimaryVertex()->GetXYZ(v);

}

//____________________________________________________________________________
Double_t AliCaloTrackAODReader::GetBField() const {
  //Return magnetic field

  Double_t bfield =  ((AliAODEvent*)fInputEvent)->GetMagneticField();

  return bfield;

}

//____________________________________________________________________________________________________________________________________________________
Int_t AliCaloTrackAODReader::GetModuleNumber(AliAODCaloCluster * cluster) const
{
	//Get the EMCAL/PHOS module number that corresponds to this cluster
	TLorentzVector lv;
	Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
	cluster->GetMomentum(lv,v);
	Float_t phi = lv.Phi();
	if(phi < 0) phi+=TMath::TwoPi();	
	Int_t absId = -1;
	if(cluster->IsEMCALCluster()){
		GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
		if(GetDebug() > 2) 
			printf("AliCaloTrackAODReader::GetModuleNumber(ESD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
				   lv.Eta(), phi*TMath::RadToDeg(),absId, GetEMCALGeometry()->GetSuperModuleNumber(absId));
		return GetEMCALGeometry()->GetSuperModuleNumber(absId) ;
	}//EMCAL
	else if(cluster->IsPHOSCluster()) {
		Int_t    relId[4];
		if ( cluster->GetNCells() > 0) {
			absId = cluster->GetCellAbsId(0);
			if(GetDebug() > 2) 
				printf("AliCaloTrackAODReader::GetModuleNumber(AOD) - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
					   lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
		}
		else return -1;
		
		if ( absId >= 0) {
			GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
			if(GetDebug() > 2) 
				printf("AliCaloTrackAODReader::GetModuleNumber(AOD) - PHOS: Module %d\n",relId[0]-1);
			return relId[0]-1;
		}
		else return -1;
	}//PHOS
	
	return -1;
}


//____________________________________________________________________________
void AliCaloTrackAODReader::SetInputOutputMCEvent(AliVEvent* input, AliAODEvent* aod, AliMCEvent* mc) {
  // Connect the data pointers
  // If input is AOD, do analysis with input, if not, do analysis with the output aod.

  if(!strcmp(input->GetName(),"AliESDEvent"))   {
    SetInputEvent(aod);
    SetOutputEvent(aod);
  }
  else if(!strcmp(input->GetName(),"AliAODEvent")){
	  
	  AliAODInputHandler* aodIH = static_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
	  //printf("AODInputHandler %p, MergeEvents %d \n",aodIH, aodIH->GetMergeEvents());
	  if (aodIH && aodIH->GetMergeEvents()) {
		  //Merged events, use output AOD.
		  SetInputEvent(aod);
		  SetOutputEvent(aod);
	  }
	  else{
		  SetInputEvent(input);
		  SetOutputEvent(aod);
	  }
  }
  else{ 
    printf("AliCaloTrackAODReader::SetInputOutputMCEvent() - STOP : Wrong data format: %s\n",input->GetName());
    abort();
  }
  
  SetMC(mc);
  
}

//________________________________________________________________
void AliCaloTrackAODReader::Print(const Option_t * opt) const
{
	
	//Print some relevant parameters set for the analysis
	AliCaloTrackReader::Print(opt);
	
	printf("Write std AOD       =     %d\n", fWriteOutputStdAOD) ;
	
	printf("    \n") ;
}
Commit	Line	Data
1c5acb87	1
	2	/**************************************************************************
	3	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	4	* *
	5	* Author: The ALICE Off-line Project. *
	6	* Contributors are mentioned in the code where appropriate. *
	7	* *
	8	* Permission to use, copy, modify and distribute this software and its *
	9	* documentation strictly for non-commercial purposes is hereby granted *
	10	* without fee, provided that the above copyright notice appears in all *
	11	* copies and that both the copyright notice and this permission notice *
	12	* appear in the supporting documentation. The authors make no claims *
	13	* about the suitability of this software for any purpose. It is *
	14	* provided "as is" without express or implied warranty. *
	15	**************************************************************************/
	16	/* $Id: $ */
	17
	18	//_________________________________________________________________________
	19	// Class for reading data (AODs) in order to do prompt gamma
591cc579	20	// or other particle identification and correlations.
	21	// Mixing analysis can be done, input AOD with events
	22	// is opened in the AliCaloTrackReader::Init()
1c5acb87	23	//
	24	//
	25	//*-- Author: Gustavo Conesa (LNF-INFN)
	26	//////////////////////////////////////////////////////////////////////////////
	27
	28
	29	// --- ROOT system ---
	30	//#include "Riostream.h"
	31
	32	//---- ANALYSIS system ----
	33	#include "AliCaloTrackAODReader.h"
1c5acb87	34	#include "AliAODCaloCluster.h"
1c5acb87	35	#include "AliAODTrack.h"
591cc579	36	#include "AliESDtrack.h"
ff45398a	37	#include "AliFiducialCut.h"
8a587055	38	#include "AliAODInputHandler.h"
8a587055	39	#include "AliAnalysisManager.h"
1c5acb87	40
	41	ClassImp(AliCaloTrackAODReader)
	42
	43	//____________________________________________________________________________
	44	AliCaloTrackAODReader::AliCaloTrackAODReader() :
08a064bc	45	AliCaloTrackReader(), fWriteOutputStdAOD(kFALSE)
1c5acb87	46	{
	47	//Default Ctor
	48
	49	//Initialize parameters
	50	fDataType=kAOD;
591cc579	51	fReadStack = kTRUE;
	52	fReadAODMCParticles = kFALSE;
	53	//We want tracks fitted in the detectors:
	54	fTrackStatus=AliESDtrack::kTPCrefit;
	55	fTrackStatus\|=AliESDtrack::kITSrefit;
	56
1c5acb87	57	}
	58
	59	//____________________________________________________________________________
a8a55c9d	60	AliCaloTrackAODReader::AliCaloTrackAODReader(const AliCaloTrackAODReader & aodr) :
08a064bc	61	AliCaloTrackReader(aodr), fWriteOutputStdAOD(aodr.fWriteOutputStdAOD)
1c5acb87	62	{
	63	// cpy ctor
	64	}
	65
	66	//_________________________________________________________________________
	67	//AliCaloTrackAODReader & AliCaloTrackAODReader::operator = (const AliCaloTrackAODReader & source)
	68	//{
	69	// // assignment operator
	70	//
	71	// if(&source == this) return *this;
	72	//
	73	// return *this;
	74	//
	75	//}
	76
	77	//____________________________________________________________________________
	78	void AliCaloTrackAODReader::FillInputCTS() {
591cc579	79	//Return array with Central Tracking System (CTS) tracks
be1f5fa4	80	if(fDebug > 2 ) printf("AliCaloTrackAODReader::FillInputCTS()\n");
477d6cee	81	Int_t nTracks = fInputEvent->GetNumberOfTracks() ;
a8a55c9d	82	Int_t naod = 0;
1c5acb87	83	Double_t p[3];
1c5acb87	84	for (Int_t itrack = 0; itrack < nTracks; itrack++) {////////////// track loop
477d6cee	85	AliAODTrack * track = ((AliAODEvent*)fInputEvent)->GetTrack(itrack) ; // retrieve track from esd
591cc579	86
	87	//Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
	88	if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus)) continue ;
1c5acb87	89
	90	track->GetPxPyPz(p) ;
	91	TLorentzVector momentum(p[0],p[1],p[2],0);
	92
ff45398a	93	if(fCTSPtMin < momentum.Pt() && fFiducialCut->IsInFiducialCut(momentum,"CTS")){
1c5acb87	94
a3aebfff	95	if(fDebug > 2 && momentum.Pt() > 0.1) printf("AliCaloTrackAODReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
1c5acb87	96	momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
a8a55c9d	97
42dc8e7d	98	if(fWriteOutputStdAOD){
a8a55c9d	99	AliAODTrack* newtrack = new(((fOutputEvent->GetTracks()))[naod++]) AliAODTrack(track);
	100	fAODCTS->Add(newtrack); //Use AOD stored in output for references.
	101	}
	102	else fAODCTS->Add(track);
ff45398a	103	}//Pt and Fiducial cut passed.
1c5acb87	104	}// track loop
591cc579	105
	106	fAODCTSNormalInputEntries = fAODCTS->GetEntriesFast();
	107	if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputCTS() - aod entries %d\n", fAODCTSNormalInputEntries);
	108
	109	//If second input event available, add the clusters.
	110	if(fSecondInputAODTree && fSecondInputAODEvent){
	111	nTracks = fSecondInputAODEvent->GetNumberOfTracks() ;
	112	if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputCTS() - Add second input tracks, entries %d\n", nTracks) ;
	113	for (Int_t itrack = 0; itrack < nTracks; itrack++) {////////////// track loop
	114	AliAODTrack * track = ((AliAODEvent*)fSecondInputAODEvent)->GetTrack(itrack) ; // retrieve track from esd
	115
	116	//Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
	117	if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus)) continue ;
	118
	119	track->GetPxPyPz(p) ;
	120	TLorentzVector momentum(p[0],p[1],p[2],0);
	121
ff45398a	122	if(fCTSPtMin < momentum.Pt() && fFiducialCut->IsInFiducialCut(momentum,"CTS")){
591cc579	123
	124	if(fDebug > 2 && momentum.Pt() > 0.1) printf("AliCaloTrackAODReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
	125	momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
	126
42dc8e7d	127	if(fWriteOutputStdAOD){
a8a55c9d	128	AliAODTrack* newtrack = new(((fOutputEvent->GetTracks()))[naod++]) AliAODTrack(track);
	129	fAODCTS->Add(newtrack); //Use AOD stored in output for references.
	130	}
	131	else fAODCTS->Add(track);
591cc579	132
ff45398a	133	}//Pt and Fiducial cut passed.
591cc579	134	}// track loop
	135
	136	if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputCTS() - aod normal entries %d, after second input %d\n", fAODCTSNormalInputEntries, fAODCTS->GetEntriesFast());
	137	} //second input loop
	138
1c5acb87	139	}
	140
	141	//____________________________________________________________________________
	142	void AliCaloTrackAODReader::FillInputEMCAL() {
	143	//Return array with EMCAL clusters in aod format
be1f5fa4	144	if(fDebug > 2 ) printf("AliCaloTrackAODReader::FillInputEMCAL()\n");
be1f5fa4	145
1c5acb87	146	//Get vertex for momentum calculation
	147	Double_t v[3] ; //vertex ;
	148	GetVertex(v);
a8a55c9d	149
a8a55c9d	150	Int_t naod = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
ff45398a	151	//Loop to select clusters in fiducial cut and fill container with aodClusters
591cc579	152	Int_t nclusters = ((AliAODEvent*)fInputEvent)->GetNCaloClusters();
591cc579	153	for (Int_t iclus = 0; iclus < nclusters; iclus++) {
477d6cee	154	AliAODCaloCluster * clus = 0;
	155	if ( (clus = ((AliAODEvent*)fInputEvent)->GetCaloCluster(iclus)) ) {
	156	if (clus->IsEMCALCluster()){
c1ac3823	157
	158	//Check if the cluster contains any bad channel
	159	if(ClusterContainsBadChannel("EMCAL",clus->GetCellsAbsId(), clus->GetNCells())) continue;
	160
477d6cee	161	TLorentzVector momentum ;
	162	clus->GetMomentum(momentum, v);
	163
ff45398a	164	if(fEMCALPtMin < momentum.Pt() && fFiducialCut->IsInFiducialCut(momentum,"EMCAL")){
1c5acb87	165
a3aebfff	166	if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackAODReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
477d6cee	167	momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
477d6cee	168
42dc8e7d	169	if(fWriteOutputStdAOD){
a8a55c9d	170	AliAODCaloCluster * newclus = new(((fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(clus);
	171	fAODEMCAL->Add(newclus);
	172	}
	173	else fAODEMCAL->Add(clus);
ff45398a	174	}//Pt and Fiducial cut passed.
477d6cee	175	}//EMCAL cluster
477d6cee	176	}// cluster exists
591cc579	177	}// cluster loop
	178
	179	fAODEMCALNormalInputEntries = fAODEMCAL->GetEntriesFast();
	180	if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputEMCAL() - aod entries %d\n", fAODEMCALNormalInputEntries);
	181
	182	//If second input event available, add the clusters.
	183	if(fSecondInputAODTree && fSecondInputAODEvent){
233e0df8	184	GetSecondInputAODVertex(v);
591cc579	185	nclusters = ((AliAODEvent*)fSecondInputAODEvent)->GetNCaloClusters();
	186	if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputEMCAL() - Add second input clusters, entries %d\n", nclusters) ;
	187	for (Int_t iclus = 0; iclus < nclusters; iclus++) {
	188	AliAODCaloCluster * clus = 0;
	189	if ( (clus = ((AliAODEvent*)fSecondInputAODEvent)->GetCaloCluster(iclus)) ) {
	190	if (clus->IsEMCALCluster()){
	191	TLorentzVector momentum ;
	192	clus->GetMomentum(momentum, v);
	193
ff45398a	194	if(fEMCALPtMin < momentum.Pt() && fFiducialCut->IsInFiducialCut(momentum,"EMCAL")){
591cc579	195
	196	if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackAODReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
	197	momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
42dc8e7d	198	if(fWriteOutputStdAOD){
a8a55c9d	199	AliAODCaloCluster * newclus = new(((fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(clus);
	200	fAODEMCAL->Add(newclus);
	201	}
	202	else fAODEMCAL->Add(clus);
ff45398a	203	}//Pt and Fiducial cut passed.
591cc579	204	}//EMCAL cluster
	205	}// cluster exists
	206	}// cluster loop
	207
	208	if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputEMCAL() - aod normal entries %d, after second input %d\n", fAODEMCALNormalInputEntries, fAODEMCAL->GetEntriesFast());
1c5acb87	209
591cc579	210	} //second input loop
1c5acb87	211	}
	212
	213	//____________________________________________________________________________
	214	void AliCaloTrackAODReader::FillInputPHOS() {
	215	//Return array with PHOS clusters in aod format
be1f5fa4	216	if(fDebug > 2 ) printf("AliCaloTrackAODReader::FillInputPHOS()\n");
be1f5fa4	217
1c5acb87	218	//Get vertex for momentum calculation
	219	Double_t v[3] ; //vertex ;
	220	GetVertex(v);
a8a55c9d	221
a8a55c9d	222	Int_t naod = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
ff45398a	223	//Loop to select clusters in fiducial cut and fill container with aodClusters
591cc579	224	Int_t nclusters = ((AliAODEvent*)fInputEvent)->GetNCaloClusters();
591cc579	225	for (Int_t iclus = 0; iclus < nclusters; iclus++) {
477d6cee	226	AliAODCaloCluster * clus = 0;
	227	if ( (clus = ((AliAODEvent*)fInputEvent)->GetCaloCluster(iclus)) ) {
	228	if (clus->IsPHOSCluster()){
c1ac3823	229
	230	//Check if the cluster contains any bad channel
	231	if(ClusterContainsBadChannel("PHOS",clus->GetCellsAbsId(), clus->GetNCells())) continue;
	232
477d6cee	233	TLorentzVector momentum ;
	234	clus->GetMomentum(momentum, v);
	235
ff45398a	236	if(fPHOSPtMin < momentum.Pt() && fFiducialCut->IsInFiducialCut(momentum,"PHOS")){
477d6cee	237
a3aebfff	238	if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackAODReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
477d6cee	239	momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
591cc579	240
42dc8e7d	241	if(fWriteOutputStdAOD){
a8a55c9d	242	AliAODCaloCluster * newclus = new(((fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(clus);
	243	fAODPHOS->Add(newclus);
	244	}
	245	else fAODPHOS->Add(clus);
ff45398a	246	}//Pt and Fiducial cut passed.
477d6cee	247	}//PHOS cluster
477d6cee	248	}//cluster exists
1c5acb87	249	}//esd cluster loop
591cc579	250
	251	fAODPHOSNormalInputEntries = fAODPHOS->GetEntriesFast() ;
	252	if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputPHOS() - aod entries %d\n", fAODPHOSNormalInputEntries);
1c5acb87	253
591cc579	254	//If second input event available, add the clusters.
233e0df8	255	if(fSecondInputAODTree && fSecondInputAODEvent){
233e0df8	256	GetSecondInputAODVertex(v);
591cc579	257	nclusters = ((AliAODEvent*)fSecondInputAODEvent)->GetNCaloClusters();
	258	if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputPHOS() - Add second input clusters, entries %d\n", nclusters);
	259	for (Int_t iclus = 0; iclus < nclusters; iclus++) {
	260	AliAODCaloCluster * clus = 0;
	261	if ( (clus = ((AliAODEvent*)fSecondInputAODEvent)->GetCaloCluster(iclus)) ) {
	262	if (clus->IsPHOSCluster()){
	263	TLorentzVector momentum ;
	264	clus->GetMomentum(momentum, v);
	265
ff45398a	266	if(fPHOSPtMin < momentum.Pt() && fFiducialCut->IsInFiducialCut(momentum,"PHOS")){
591cc579	267
	268	if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackAODReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
	269	momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
42dc8e7d	270	if(fWriteOutputStdAOD){
a8a55c9d	271	AliAODCaloCluster * newclus = new(((fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(clus);
	272	fAODPHOS->Add(newclus);
	273	}
	274	else fAODPHOS->Add(clus);
ff45398a	275	}//Pt and Fiducial cut passed.
591cc579	276	}//PHOS cluster
	277	}// cluster exists
	278	}// cluster loop
	279	if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputPHOS() - aod normal entries %d, after second input %d\n", fAODPHOSNormalInputEntries, fAODPHOS->GetEntriesFast());
	280	} //second input loop
1c5acb87	281
	282	}
	283
	284	//____________________________________________________________________________
	285	void AliCaloTrackAODReader::FillInputEMCALCells() {
	286	//Return array with EMCAL cells in aod format
	287
477d6cee	288	fEMCALCells = (TNamed) ((AliAODEvent)fInputEvent)->GetEMCALCells();
1c5acb87	289
	290	}
	291
	292	//____________________________________________________________________________
	293	void AliCaloTrackAODReader::FillInputPHOSCells() {
	294	//Return array with PHOS cells in aod format
	295
477d6cee	296	fPHOSCells = (TNamed) ((AliAODEvent)fInputEvent)->GetPHOSCells();
1c5acb87	297
	298	}
	299
	300	//____________________________________________________________________________
	301	void AliCaloTrackAODReader::GetVertex(Double_t v[3]) const {
	302	//Return vertex position
	303
8a587055	304	((AliAODEvent*)fInputEvent)->GetPrimaryVertex()->GetXYZ(v);
8a587055	305
1c5acb87	306	}
1c5acb87	307
233e0df8	308	//____________________________________________________________________________
	309	void AliCaloTrackAODReader::GetSecondInputAODVertex(Double_t v[3]) const {
	310	//Return vertex position of second AOD input
	311
8a587055	312	fSecondInputAODEvent->GetPrimaryVertex()->GetXYZ(v);
8a587055	313
233e0df8	314	}
233e0df8	315
8dacfd76	316	//____________________________________________________________________________
	317	Double_t AliCaloTrackAODReader::GetBField() const {
	318	//Return magnetic field
	319
	320	Double_t bfield = ((AliAODEvent*)fInputEvent)->GetMagneticField();
	321
	322	return bfield;
	323
	324	}
1c5acb87	325
e6e63685	326	//____________________________________________________________________________________________________________________________________________________
	327	Int_t AliCaloTrackAODReader::GetModuleNumber(AliAODCaloCluster * cluster) const
	328	{
	329	//Get the EMCAL/PHOS module number that corresponds to this cluster
	330	TLorentzVector lv;
	331	Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
	332	cluster->GetMomentum(lv,v);
	333	Float_t phi = lv.Phi();
	334	if(phi < 0) phi+=TMath::TwoPi();
	335	Int_t absId = -1;
	336	if(cluster->IsEMCALCluster()){
	337	GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
	338	if(GetDebug() > 2)
	339	printf("AliCaloTrackAODReader::GetModuleNumber(ESD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
	340	lv.Eta(), phi*TMath::RadToDeg(),absId, GetEMCALGeometry()->GetSuperModuleNumber(absId));
	341	return GetEMCALGeometry()->GetSuperModuleNumber(absId) ;
	342	}//EMCAL
	343	else if(cluster->IsPHOSCluster()) {
	344	Int_t relId[4];
	345	if ( cluster->GetNCells() > 0) {
	346	absId = cluster->GetCellAbsId(0);
	347	if(GetDebug() > 2)
	348	printf("AliCaloTrackAODReader::GetModuleNumber(AOD) - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
	349	lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
	350	}
	351	else return -1;
	352
	353	if ( absId >= 0) {
	354	GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
	355	if(GetDebug() > 2)
	356	printf("AliCaloTrackAODReader::GetModuleNumber(AOD) - PHOS: Module %d\n",relId[0]-1);
	357	return relId[0]-1;
	358	}
	359	else return -1;
	360	}//PHOS
	361
	362	return -1;
	363	}
	364
	365
1c5acb87	366	//____________________________________________________________________________
477d6cee	367	void AliCaloTrackAODReader::SetInputOutputMCEvent(AliVEvent* input, AliAODEvent* aod, AliMCEvent* mc) {
1c5acb87	368	// Connect the data pointers
477d6cee	369	// If input is AOD, do analysis with input, if not, do analysis with the output aod.
8a587055	370
477d6cee	371	if(!strcmp(input->GetName(),"AliESDEvent")) {
	372	SetInputEvent(aod);
	373	SetOutputEvent(aod);
	374	}
	375	else if(!strcmp(input->GetName(),"AliAODEvent")){
8a587055	376
	377	AliAODInputHandler* aodIH = static_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
	378	//printf("AODInputHandler %p, MergeEvents %d \n",aodIH, aodIH->GetMergeEvents());
	379	if (aodIH && aodIH->GetMergeEvents()) {
	380	//Merged events, use output AOD.
	381	SetInputEvent(aod);
	382	SetOutputEvent(aod);
	383	}
	384	else{
	385	SetInputEvent(input);
	386	SetOutputEvent(aod);
	387	}
477d6cee	388	}
477d6cee	389	else{
591cc579	390	printf("AliCaloTrackAODReader::SetInputOutputMCEvent() - STOP : Wrong data format: %s\n",input->GetName());
477d6cee	391	abort();
	392	}
	393
	394	SetMC(mc);
1c5acb87	395
1c5acb87	396	}
08a064bc	397
	398	//________________________________________________________________
	399	void AliCaloTrackAODReader::Print(const Option_t * opt) const
	400	{
	401
	402	//Print some relevant parameters set for the analysis
	403	AliCaloTrackReader::Print(opt);
	404
	405	printf("Write std AOD = %d\n", fWriteOutputStdAOD) ;
	406
	407	printf(" \n") ;
	408	}
	409