[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / AliCaloTrackReader.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:  $ */

//_________________________________________________________________________
// Base class for reading data: MonteCarlo, ESD or AOD, of PHOS EMCAL and 
// Central Barrel Tracking detectors (CTS).
// Not all MC particles/tracks/clusters are kept, some kinematical/fidutial restrictions are done.
// Mother class of : AliCaloTrackESDReader: Fills ESD data in 3 TObjArrays (PHOS, EMCAL, CTS)
//                 : AliCaloTrackMCReader: Fills Kinematics data in 3 TObjArrays (PHOS, EMCAL, CTS)
//                 : AliCaloTrackAODReader: Fills AOD data in 3 TObjArrays (PHOS, EMCAL, CTS) 
//                
//-- Author: Gustavo Conesa (LNF-INFN) 
//////////////////////////////////////////////////////////////////////////////


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

//---- ANALYSIS system ----
#include "AliCaloTrackReader.h"
#include "AliFidutialCut.h"
#include "AliMCEvent.h"
#include "AliAODMCHeader.h"
#include "AliGenPythiaEventHeader.h"
#include "AliAODEvent.h"

ClassImp(AliCaloTrackReader)
  
  
//____________________________________________________________________________
  AliCaloTrackReader::AliCaloTrackReader() : 
    TObject(), fEventNumber(-1), fCurrentFileName(""),fDataType(0), fDebug(0), 
    fFidutialCut(0x0), fComparePtHardAndJetPt(kFALSE), fPtHardAndJetPtFactor(7),
    fCTSPtMin(0), fEMCALPtMin(0),fPHOSPtMin(0),
    fAODCTS(new TObjArray()), fAODEMCAL(new TObjArray()), fAODPHOS(new TObjArray()),
    fEMCALCells(0x0), fPHOSCells(0x0),
    fInputEvent(0x0), fOutputEvent(0x0),fMC(0x0),
    fFillCTS(0),fFillEMCAL(0),fFillPHOS(0),
    fFillEMCALCells(0),fFillPHOSCells(0), 
	fSecondInputAODTree(0x0), fSecondInputAODEvent(0x0),
	fSecondInputFileName(""),fSecondInputFirstEvent(0), 
	fAODCTSNormalInputEntries(0), fAODEMCALNormalInputEntries(0), 
    fAODPHOSNormalInputEntries(0), fTrackStatus(0), 
	fReadStack(kFALSE), fReadAODMCParticles(kFALSE), 
	fCleanOutputStdAOD(kFALSE), fDeltaAODFileName("deltaAODPartCorr.root")
{
  //Ctor
  
  //Initialize parameters
  InitParameters();
}

//____________________________________________________________________________
AliCaloTrackReader::AliCaloTrackReader(const AliCaloTrackReader & g) :   
  TObject(g), fEventNumber(g.fEventNumber), fCurrentFileName(g.fCurrentFileName), 
  fDataType(g.fDataType), fDebug(g.fDebug),
  fFidutialCut(g.fFidutialCut),
  fComparePtHardAndJetPt(g.fComparePtHardAndJetPt),
  fPtHardAndJetPtFactor(g.fPtHardAndJetPtFactor),
  fCTSPtMin(g.fCTSPtMin), fEMCALPtMin(g.fEMCALPtMin),fPHOSPtMin(g.fPHOSPtMin), 
  fAODCTS(new TObjArray(*g.fAODCTS)),  
  fAODEMCAL(new TObjArray(*g.fAODEMCAL)),
  fAODPHOS(new TObjArray(*g.fAODPHOS)),
  fEMCALCells(new TNamed(*g.fEMCALCells)),
  fPHOSCells(new TNamed(*g.fPHOSCells)),
  fInputEvent(g.fInputEvent), fOutputEvent(g.fOutputEvent), fMC(g.fMC),
  fFillCTS(g.fFillCTS),fFillEMCAL(g.fFillEMCAL),fFillPHOS(g.fFillPHOS),
  fFillEMCALCells(g.fFillEMCALCells),fFillPHOSCells(g.fFillPHOSCells),
  fSecondInputAODTree(g.fSecondInputAODTree), 
  fSecondInputAODEvent(g.fSecondInputAODEvent),
  fSecondInputFileName(g.fSecondInputFileName), 
  fSecondInputFirstEvent(g.fSecondInputFirstEvent),
  fAODCTSNormalInputEntries(g.fAODCTSNormalInputEntries), 
  fAODEMCALNormalInputEntries(g.fAODEMCALNormalInputEntries), 
  fAODPHOSNormalInputEntries(g.fAODPHOSNormalInputEntries),
  fTrackStatus(g.fTrackStatus),
  fReadStack(g.fReadStack), fReadAODMCParticles(g.fReadAODMCParticles),
  fCleanOutputStdAOD(g.fCleanOutputStdAOD), fDeltaAODFileName(g.fDeltaAODFileName)
{
  // cpy ctor
  
}

//_________________________________________________________________________
AliCaloTrackReader & AliCaloTrackReader::operator = (const AliCaloTrackReader & source)
{
  // assignment operator
  
  if(&source == this) return *this;
  
  fDataType    = source.fDataType ;
  fDebug       = source.fDebug ;
  fEventNumber = source.fEventNumber ;
  fCurrentFileName = source.fCurrentFileName ;
  fFidutialCut = source.fFidutialCut;
	
  fComparePtHardAndJetPt = source.fComparePtHardAndJetPt;
  fPtHardAndJetPtFactor  = source.fPtHardAndJetPtFactor;
	
  fCTSPtMin    = source.fCTSPtMin ;
  fEMCALPtMin  = source.fEMCALPtMin ;
  fPHOSPtMin   = source.fPHOSPtMin ; 
  
  fAODCTS     = new TObjArray(*source.fAODCTS) ;
  fAODEMCAL   = new TObjArray(*source.fAODEMCAL) ;
  fAODPHOS    = new TObjArray(*source.fAODPHOS) ;
  fEMCALCells = new TNamed(*source.fEMCALCells) ;
  fPHOSCells  = new TNamed(*source.fPHOSCells) ;

  fInputEvent  = source.fInputEvent;
  fOutputEvent = source.fOutputEvent;
  fMC          = source.fMC;
  
  fFillCTS        = source.fFillCTS;
  fFillEMCAL      = source.fFillEMCAL;
  fFillPHOS       = source.fFillPHOS;
  fFillEMCALCells = source.fFillEMCALCells;
  fFillPHOSCells  = source.fFillPHOSCells;

  fSecondInputAODTree    = source.fSecondInputAODTree;
  fSecondInputAODEvent   = source.fSecondInputAODEvent;
  fSecondInputFileName   = source.fSecondInputFileName;
  fSecondInputFirstEvent = source.fSecondInputFirstEvent;

  fAODCTSNormalInputEntries   = source.fAODCTSNormalInputEntries; 
  fAODEMCALNormalInputEntries = source.fAODEMCALNormalInputEntries; 
  fAODPHOSNormalInputEntries  = source.fAODPHOSNormalInputEntries;
	
  fTrackStatus        = source.fTrackStatus;
  fReadStack          = source.fReadStack;
  fReadAODMCParticles = source.fReadAODMCParticles;	
	
  fCleanOutputStdAOD  = source.fCleanOutputStdAOD;
  fDeltaAODFileName   = source.fDeltaAODFileName;
	
  return *this;
  
}

//_________________________________
AliCaloTrackReader::~AliCaloTrackReader() {
  //Dtor
  
  if(fFidutialCut) delete fFidutialCut ;
	
  if(fAODCTS){
    fAODCTS->Clear() ; 
    delete fAODCTS ;
  }
  
  if(fAODEMCAL){
    fAODEMCAL->Clear() ; 
    delete fAODEMCAL ;
  }
  
  if(fAODPHOS){
    fAODPHOS->Clear() ; 
    delete fAODPHOS ;
  }
  
  if(fEMCALCells){
    fEMCALCells->Clear() ; 
    delete fEMCALCells ;
  }
  
  if(fPHOSCells){
    fPHOSCells->Clear() ; 
    delete fPHOSCells ;
  }

  if(fInputEvent)  delete fInputEvent ;
  if(fOutputEvent) delete fOutputEvent ;
  if(fMC)          delete fMC ;  
	
  if(fSecondInputAODTree){
	fSecondInputAODTree->Clear();
	delete fSecondInputAODTree;
  }
	
  if(fSecondInputAODEvent) delete fSecondInputAODEvent ;

}


//_________________________________________________________________________
Bool_t AliCaloTrackReader::ComparePtHardAndJetPt(){
	// Check the event, if the requested ptHard is much larger than the jet pT, then there is a problem.
	// Only for PYTHIA.
	if(!fReadStack) return kTRUE; //Information not filtered to AOD
	
	if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader")){
		TParticle * jet =  new TParticle;
		AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
		Int_t nTriggerJets =  pygeh->NTriggerJets();
		Float_t ptHard = pygeh->GetPtHard();
		
		//if(fDebug > 1) printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %f\n",nTriggerJets, ptHard);
	    Float_t tmpjet[]={0,0,0,0};
		for(Int_t ijet = 0; ijet< nTriggerJets; ijet++){
			pygeh->TriggerJet(ijet, tmpjet);
			jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
			//Compare jet pT and pt Hard
			//if(fDebug > 1) printf("AliMCAnalysisUtils:: %d pycell jet pT %f\n",ijet, jet->Pt());
			if(jet->Pt() > fPtHardAndJetPtFactor * ptHard) {
				printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n",
					   nTriggerJets, ptHard, jet->Pt(), fPtHardAndJetPtFactor);
				return kFALSE;
			}
		}
	}
	
	return kTRUE ;
	
}

//____________________________________________________________________________
AliStack* AliCaloTrackReader::GetStack() const {
  //Return pointer to stack
  if(fMC)
    return fMC->Stack();
  else{
    if(fDebug > 1) printf("AliCaloTrackReader::GetStack() - Stack is not available\n"); 
    return 0x0 ;
  }
}

//____________________________________________________________________________
AliHeader* AliCaloTrackReader::GetHeader() const {
  //Return pointer to header
  if(fMC)
    return fMC->Header();
  else{
    printf("AliCaloTrackReader::Header is not available\n"); 
    return 0x0 ;
  }
}
//____________________________________________________________________________
AliGenEventHeader* AliCaloTrackReader::GetGenEventHeader() const {
  //Return pointer to Generated event header
  if(fMC)
    return fMC->GenEventHeader();
  else{
    printf("AliCaloTrackReader::GenEventHeader is not available\n"); 
    return 0x0 ;
  }
}

//____________________________________________________________________________
TClonesArray* AliCaloTrackReader::GetAODMCParticles(Int_t input) const {
	//Return list of particles in AOD. Do it for the corresponding input event.
	if(fDataType == kAOD){
	 //Normal input AOD
	 if(input == 0) return (TClonesArray*)((AliAODEvent*)fInputEvent)->FindListObject("mcparticles");
	  //Second input AOD
	 else if(input == 1 && fSecondInputAODEvent) return (TClonesArray*) fSecondInputAODEvent->FindListObject("mcparticles");	
	 else {
	     printf("AliCaloTrackReader::GetAODMCParticles() - wrong AOD input index? %d, or non existing tree? \n",input); 
		 return 0x0;
	 }
	}
	else {
		printf("AliCaloTrackReader::GetAODMCParticles() - Input are not AODs\n"); 
		return 0x0;
	}
}

//____________________________________________________________________________
AliAODMCHeader* AliCaloTrackReader::GetAODMCHeader(Int_t input) const {
	//Return MC header in AOD. Do it for the corresponding input event.
	if(fDataType == kAOD){
		//Normal input AOD
		if(input == 0) return (AliAODMCHeader*)((AliAODEvent*)fInputEvent)->FindListObject("mcheader");
		//Second input AOD
		else if(input == 1) return  (AliAODMCHeader*) fSecondInputAODEvent->FindListObject("mcheader");	
		else {
			printf("AliCaloTrackReader::GetAODMCHeader() - wrong AOD input index, %d\n",input);
			return 0x0;
		}
	}
	else {
		printf("AliCaloTrackReader::GetAODMCHeader() - Input are not AODs\n");
		return 0x0;
	}
}

//_______________________________________________________________
void AliCaloTrackReader::Init()
{
	//Init reader. Method to be called in AliAnaPartCorrMaker
	
	//Get the file with second input events if the filename is given
	//Get the tree and connect the AODEvent. Only with AODs

	if(fReadStack && fReadAODMCParticles){
		printf("AliCaloTrackReader::Init() - Cannot access stack and mcparticles at the same time, change them \n");
		fReadStack = kFALSE;
		fReadAODMCParticles = kFALSE;
	}
	
	if(fSecondInputFileName!=""){
		if(fDataType == kAOD){
			TFile * input2   = new TFile(fSecondInputFileName,"read");
			printf("AliCaloTrackReader::Init() - Second input file opened: %s, size %d \n", input2->GetName(), (Int_t) input2->GetSize());
			fSecondInputAODTree = (TTree*) input2->Get("aodTree");
			if(fSecondInputAODTree) printf("AliCaloTrackReader::Init() - Second input tree opened: %s, entries %d \n", 
										   fSecondInputAODTree->GetName(), (Int_t) fSecondInputAODTree->GetEntries());
			else{
			 printf("AliCaloTrackReader::Init() - Second input tree not available, STOP \n");
			 abort();
			}
			fSecondInputAODEvent = new AliAODEvent;
			fSecondInputAODEvent->ReadFromTree(fSecondInputAODTree);
			if(fSecondInputFirstEvent >= fSecondInputAODTree->GetEntriesFast()){
				printf("AliCaloTrackReader::Init() - Requested first event of second input %d, is larger than number of events %d, STOP\n", 
					   fSecondInputFirstEvent, (Int_t) fSecondInputAODTree->GetEntriesFast());
				abort();
			}
		}
		else printf("AliCaloTrackReader::Init() - Second input not added, reader is not AOD\n");
	}
	
	
}
//_______________________________________________________________
void AliCaloTrackReader::InitParameters()
{
 
  //Initialize the parameters of the analysis.
  fDataType   = kESD ;
  fCTSPtMin   = 0.2 ;
  fEMCALPtMin = 0.2 ;
  fPHOSPtMin  = 0.2 ;

  //Do not filter the detectors input by default.
  fFillEMCAL      = kFALSE;
  fFillPHOS       = kFALSE;
  fFillCTS        = kFALSE;
  fFillEMCALCells = kFALSE;
  fFillPHOSCells  = kFALSE;

  fFidutialCut           = new AliFidutialCut();
  fSecondInputFileName   = "" ;
  fSecondInputFirstEvent = 0 ;
  fReadStack             = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
  fReadAODMCParticles    = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
  fCleanOutputStdAOD     = kFALSE; // Clean the standard clusters/tracks?
  fDeltaAODFileName      = "deltaAODPartCorr.root";
}


//________________________________________________________________
void AliCaloTrackReader::Print(const Option_t * opt) const
{

  //Print some relevant parameters set for the analysis
  if(! opt)
    return;

  printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
  printf("Data type      : %d\n", fDataType) ;
  printf("CTS Min pT     : %2.1f GeV/c\n", fCTSPtMin) ;
  printf("EMCAL Min pT   : %2.1f GeV/c\n", fEMCALPtMin) ;
  printf("PHOS Min pT    : %2.1f GeV/c\n", fPHOSPtMin) ;
  printf("Use CTS         =     %d\n", fFillCTS) ;
  printf("Use EMCAL       =     %d\n", fFillEMCAL) ;
  printf("Use PHOS        =     %d\n", fFillPHOS) ;
  printf("Use EMCAL Cells =     %d\n", fFillEMCALCells) ;
  printf("Use PHOS  Cells =     %d\n", fFillPHOSCells) ;
  printf("Track status    =     %d\n", (Int_t) fTrackStatus) ;
  if(fComparePtHardAndJetPt)
	  printf("Compare jet pt and pt hard to accept event, factor = %2.2f",fPtHardAndJetPtFactor);
	
  if(fSecondInputFileName!="") {
	  printf("Second Input File Name     =     %s\n", fSecondInputFileName.Data()) ;
	  printf("Second Input First Event   =     %d\n", fSecondInputFirstEvent) ;
  }
	
  printf("Read Kine from, stack? %d, AOD ? %d \n", fReadStack, fReadAODMCParticles) ;
  printf("Clean std AOD       =     %d\n", fCleanOutputStdAOD) ;
  printf("Delta AOD File Name =     %s\n", fDeltaAODFileName.Data()) ;

  printf("    \n") ;
} 

//___________________________________________________
Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry, const char * currentFileName) {
  //Fill the event counter and input lists that are needed, called by the analysis maker.

  fEventNumber = iEntry;
  fCurrentFileName = TString(currentFileName);
	
  if(fOutputEvent && (fDataType != kAOD) && ((fOutputEvent->GetCaloClusters())->GetEntriesFast()!=0 ||(fOutputEvent->GetTracks())->GetEntriesFast()!=0)){
    printf("AliCaloTrackReader::AODCaloClusters or AODTracks already filled by the filter, do not use the ESD reader, use the AOD reader, STOP\n");
    abort();
  }

  //In case of analysis of events with jets, skip those with jet pt > 5 pt hard	
  if(fComparePtHardAndJetPt && GetStack()) {
		if(!ComparePtHardAndJetPt()) return kFALSE ;
  }

  //In case of mixing events with other AOD file	
  if(fDataType == kAOD && fSecondInputAODTree){
	 
	 if(fDebug > 1) 
		 printf("AliCaloTrackReader::FillInputEvent() - Get event %d from second input AOD file \n", iEntry+fSecondInputFirstEvent);
	 if(fSecondInputAODTree->GetEntriesFast() <= iEntry+fSecondInputFirstEvent) {
		 if(fSecondInputAODTree->GetEntriesFast() == iEntry+fSecondInputFirstEvent) 
			 printf("AliCaloTrackReader::FillInputEvent() - Skip events from event %d, no more events in second AOD file \n", iEntry);
		 return kFALSE;
	 }
	  printf("Reader 4 \n");

	 //Get the Event
	 Int_t nbytes = fSecondInputAODTree->GetEvent(iEntry+fSecondInputFirstEvent);
	 if ( nbytes == 0 ) {//If nothing in AOD
		 printf("AliCaloTrackReader::FillInputEvent() - Nothing in Second AOD input, STOP\n");
		 abort() ; 
	 }
	  
  }

  if(fFillCTS)   FillInputCTS();
  if(fFillEMCAL) FillInputEMCAL();
  if(fFillPHOS)  FillInputPHOS();
  if(fFillEMCALCells) FillInputEMCALCells();
  if(fFillPHOSCells)  FillInputPHOSCells();

  return kTRUE ;
}

//__________________________________________________
void AliCaloTrackReader::ResetLists() {
  //  Reset lists, called by the analysis maker 

  if(fAODCTS)   fAODCTS -> Clear();
  if(fAODEMCAL) fAODEMCAL -> Clear();
  if(fAODPHOS)  fAODPHOS -> Clear();
  if(fEMCALCells) fEMCALCells -> Clear();
  if(fPHOSCells)  fPHOSCells -> Clear();
  if(fCleanOutputStdAOD && fOutputEvent ){
	  //Only keep copied tracks and clusters if requested
	  fOutputEvent->GetTracks()      ->Clear();
	  fOutputEvent->GetCaloClusters()->Clear();
  }
}
Commit	Line	Data
1c5acb87	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, 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	/* $Id: $ */
	16
	17	//_________________________________________________________________________
	18	// Base class for reading data: MonteCarlo, ESD or AOD, of PHOS EMCAL and
	19	// Central Barrel Tracking detectors (CTS).
	20	// Not all MC particles/tracks/clusters are kept, some kinematical/fidutial restrictions are done.
591cc579	21	// Mother class of : AliCaloTrackESDReader: Fills ESD data in 3 TObjArrays (PHOS, EMCAL, CTS)
	22	// : AliCaloTrackMCReader: Fills Kinematics data in 3 TObjArrays (PHOS, EMCAL, CTS)
	23	// : AliCaloTrackAODReader: Fills AOD data in 3 TObjArrays (PHOS, EMCAL, CTS)
1c5acb87	24	//
	25	//-- Author: Gustavo Conesa (LNF-INFN)
	26	//////////////////////////////////////////////////////////////////////////////
	27
	28
	29	// --- ROOT system ---
591cc579	30	#include "TFile.h"
1c5acb87	31
	32	//---- ANALYSIS system ----
	33	#include "AliCaloTrackReader.h"
477d6cee	34	#include "AliFidutialCut.h"
477d6cee	35	#include "AliMCEvent.h"
591cc579	36	#include "AliAODMCHeader.h"
591cc579	37	#include "AliGenPythiaEventHeader.h"
8dacfd76	38	#include "AliAODEvent.h"
1c5acb87	39
	40	ClassImp(AliCaloTrackReader)
	41
	42
	43	//____________________________________________________________________________
	44	AliCaloTrackReader::AliCaloTrackReader() :
a79a2424	45	TObject(), fEventNumber(-1), fCurrentFileName(""),fDataType(0), fDebug(0),
591cc579	46	fFidutialCut(0x0), fComparePtHardAndJetPt(kFALSE), fPtHardAndJetPtFactor(7),
1c5acb87	47	fCTSPtMin(0), fEMCALPtMin(0),fPHOSPtMin(0),
591cc579	48	fAODCTS(new TObjArray()), fAODEMCAL(new TObjArray()), fAODPHOS(new TObjArray()),
1c5acb87	49	fEMCALCells(0x0), fPHOSCells(0x0),
477d6cee	50	fInputEvent(0x0), fOutputEvent(0x0),fMC(0x0),
1c5acb87	51	fFillCTS(0),fFillEMCAL(0),fFillPHOS(0),
591cc579	52	fFillEMCALCells(0),fFillPHOSCells(0),
	53	fSecondInputAODTree(0x0), fSecondInputAODEvent(0x0),
	54	fSecondInputFileName(""),fSecondInputFirstEvent(0),
	55	fAODCTSNormalInputEntries(0), fAODEMCALNormalInputEntries(0),
	56	fAODPHOSNormalInputEntries(0), fTrackStatus(0),
42dc8e7d	57	fReadStack(kFALSE), fReadAODMCParticles(kFALSE),
08a064bc	58	fCleanOutputStdAOD(kFALSE), fDeltaAODFileName("deltaAODPartCorr.root")
1c5acb87	59	{
	60	//Ctor
	61
	62	//Initialize parameters
	63	InitParameters();
	64	}
	65
	66	//____________________________________________________________________________
	67	AliCaloTrackReader::AliCaloTrackReader(const AliCaloTrackReader & g) :
a79a2424	68	TObject(g), fEventNumber(g.fEventNumber), fCurrentFileName(g.fCurrentFileName),
a79a2424	69	fDataType(g.fDataType), fDebug(g.fDebug),
6639984f	70	fFidutialCut(g.fFidutialCut),
591cc579	71	fComparePtHardAndJetPt(g.fComparePtHardAndJetPt),
591cc579	72	fPtHardAndJetPtFactor(g.fPtHardAndJetPtFactor),
1c5acb87	73	fCTSPtMin(g.fCTSPtMin), fEMCALPtMin(g.fEMCALPtMin),fPHOSPtMin(g.fPHOSPtMin),
591cc579	74	fAODCTS(new TObjArray(*g.fAODCTS)),
	75	fAODEMCAL(new TObjArray(*g.fAODEMCAL)),
	76	fAODPHOS(new TObjArray(*g.fAODPHOS)),
1c5acb87	77	fEMCALCells(new TNamed(*g.fEMCALCells)),
1c5acb87	78	fPHOSCells(new TNamed(*g.fPHOSCells)),
477d6cee	79	fInputEvent(g.fInputEvent), fOutputEvent(g.fOutputEvent), fMC(g.fMC),
1c5acb87	80	fFillCTS(g.fFillCTS),fFillEMCAL(g.fFillEMCAL),fFillPHOS(g.fFillPHOS),
591cc579	81	fFillEMCALCells(g.fFillEMCALCells),fFillPHOSCells(g.fFillPHOSCells),
	82	fSecondInputAODTree(g.fSecondInputAODTree),
	83	fSecondInputAODEvent(g.fSecondInputAODEvent),
	84	fSecondInputFileName(g.fSecondInputFileName),
	85	fSecondInputFirstEvent(g.fSecondInputFirstEvent),
	86	fAODCTSNormalInputEntries(g.fAODCTSNormalInputEntries),
	87	fAODEMCALNormalInputEntries(g.fAODEMCALNormalInputEntries),
	88	fAODPHOSNormalInputEntries(g.fAODPHOSNormalInputEntries),
	89	fTrackStatus(g.fTrackStatus),
42dc8e7d	90	fReadStack(g.fReadStack), fReadAODMCParticles(g.fReadAODMCParticles),
08a064bc	91	fCleanOutputStdAOD(g.fCleanOutputStdAOD), fDeltaAODFileName(g.fDeltaAODFileName)
1c5acb87	92	{
	93	// cpy ctor
	94
	95	}
	96
	97	//_________________________________________________________________________
	98	AliCaloTrackReader & AliCaloTrackReader::operator = (const AliCaloTrackReader & source)
	99	{
	100	// assignment operator
	101
	102	if(&source == this) return *this;
	103
6639984f	104	fDataType = source.fDataType ;
	105	fDebug = source.fDebug ;
	106	fEventNumber = source.fEventNumber ;
a79a2424	107	fCurrentFileName = source.fCurrentFileName ;
1c5acb87	108	fFidutialCut = source.fFidutialCut;
591cc579	109
29b2ceec	110	fComparePtHardAndJetPt = source.fComparePtHardAndJetPt;
591cc579	111	fPtHardAndJetPtFactor = source.fPtHardAndJetPtFactor;
591cc579	112
29b2ceec	113	fCTSPtMin = source.fCTSPtMin ;
	114	fEMCALPtMin = source.fEMCALPtMin ;
	115	fPHOSPtMin = source.fPHOSPtMin ;
1c5acb87	116
591cc579	117	fAODCTS = new TObjArray(*source.fAODCTS) ;
	118	fAODEMCAL = new TObjArray(*source.fAODEMCAL) ;
	119	fAODPHOS = new TObjArray(*source.fAODPHOS) ;
1c5acb87	120	fEMCALCells = new TNamed(*source.fEMCALCells) ;
6639984f	121	fPHOSCells = new TNamed(*source.fPHOSCells) ;
1c5acb87	122
477d6cee	123	fInputEvent = source.fInputEvent;
	124	fOutputEvent = source.fOutputEvent;
	125	fMC = source.fMC;
1c5acb87	126
6639984f	127	fFillCTS = source.fFillCTS;
	128	fFillEMCAL = source.fFillEMCAL;
	129	fFillPHOS = source.fFillPHOS;
1c5acb87	130	fFillEMCALCells = source.fFillEMCALCells;
6639984f	131	fFillPHOSCells = source.fFillPHOSCells;
1c5acb87	132
591cc579	133	fSecondInputAODTree = source.fSecondInputAODTree;
	134	fSecondInputAODEvent = source.fSecondInputAODEvent;
	135	fSecondInputFileName = source.fSecondInputFileName;
	136	fSecondInputFirstEvent = source.fSecondInputFirstEvent;
	137
	138	fAODCTSNormalInputEntries = source.fAODCTSNormalInputEntries;
	139	fAODEMCALNormalInputEntries = source.fAODEMCALNormalInputEntries;
	140	fAODPHOSNormalInputEntries = source.fAODPHOSNormalInputEntries;
	141
	142	fTrackStatus = source.fTrackStatus;
	143	fReadStack = source.fReadStack;
	144	fReadAODMCParticles = source.fReadAODMCParticles;
	145
08a064bc	146	fCleanOutputStdAOD = source.fCleanOutputStdAOD;
42dc8e7d	147	fDeltaAODFileName = source.fDeltaAODFileName;
42dc8e7d	148
1c5acb87	149	return *this;
	150
	151	}
	152
	153	//_________________________________
	154	AliCaloTrackReader::~AliCaloTrackReader() {
	155	//Dtor
	156
	157	if(fFidutialCut) delete fFidutialCut ;
29b2ceec	158
1c5acb87	159	if(fAODCTS){
	160	fAODCTS->Clear() ;
	161	delete fAODCTS ;
	162	}
	163
	164	if(fAODEMCAL){
	165	fAODEMCAL->Clear() ;
	166	delete fAODEMCAL ;
	167	}
	168
	169	if(fAODPHOS){
	170	fAODPHOS->Clear() ;
	171	delete fAODPHOS ;
	172	}
	173
	174	if(fEMCALCells){
	175	fEMCALCells->Clear() ;
	176	delete fEMCALCells ;
	177	}
	178
	179	if(fPHOSCells){
	180	fPHOSCells->Clear() ;
	181	delete fPHOSCells ;
	182	}
	183
477d6cee	184	if(fInputEvent) delete fInputEvent ;
477d6cee	185	if(fOutputEvent) delete fOutputEvent ;
591cc579	186	if(fMC) delete fMC ;
	187
	188	if(fSecondInputAODTree){
	189	fSecondInputAODTree->Clear();
	190	delete fSecondInputAODTree;
	191	}
	192
	193	if(fSecondInputAODEvent) delete fSecondInputAODEvent ;
	194
1c5acb87	195	}
1c5acb87	196
477d6cee	197
591cc579	198	//_________________________________________________________________________
	199	Bool_t AliCaloTrackReader::ComparePtHardAndJetPt(){
	200	// Check the event, if the requested ptHard is much larger than the jet pT, then there is a problem.
	201	// Only for PYTHIA.
	202	if(!fReadStack) return kTRUE; //Information not filtered to AOD
	203
	204	if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader")){
	205	TParticle * jet = new TParticle;
	206	AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
	207	Int_t nTriggerJets = pygeh->NTriggerJets();
	208	Float_t ptHard = pygeh->GetPtHard();
	209
	210	//if(fDebug > 1) printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %f\n",nTriggerJets, ptHard);
	211	Float_t tmpjet[]={0,0,0,0};
	212	for(Int_t ijet = 0; ijet< nTriggerJets; ijet++){
	213	pygeh->TriggerJet(ijet, tmpjet);
	214	jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
	215	//Compare jet pT and pt Hard
	216	//if(fDebug > 1) printf("AliMCAnalysisUtils:: %d pycell jet pT %f\n",ijet, jet->Pt());
	217	if(jet->Pt() > fPtHardAndJetPtFactor * ptHard) {
	218	printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n",
	219	nTriggerJets, ptHard, jet->Pt(), fPtHardAndJetPtFactor);
	220	return kFALSE;
	221	}
	222	}
	223	}
	224
	225	return kTRUE ;
	226
	227	}
	228
1c5acb87	229	//____________________________________________________________________________
	230	AliStack* AliCaloTrackReader::GetStack() const {
	231	//Return pointer to stack
	232	if(fMC)
	233	return fMC->Stack();
	234	else{
477d6cee	235	if(fDebug > 1) printf("AliCaloTrackReader::GetStack() - Stack is not available\n");
1c5acb87	236	return 0x0 ;
	237	}
	238	}
	239
	240	//____________________________________________________________________________
	241	AliHeader* AliCaloTrackReader::GetHeader() const {
	242	//Return pointer to header
	243	if(fMC)
	244	return fMC->Header();
	245	else{
477d6cee	246	printf("AliCaloTrackReader::Header is not available\n");
1c5acb87	247	return 0x0 ;
	248	}
	249	}
	250	//____________________________________________________________________________
	251	AliGenEventHeader* AliCaloTrackReader::GetGenEventHeader() const {
	252	//Return pointer to Generated event header
	253	if(fMC)
	254	return fMC->GenEventHeader();
	255	else{
477d6cee	256	printf("AliCaloTrackReader::GenEventHeader is not available\n");
1c5acb87	257	return 0x0 ;
	258	}
	259	}
	260
591cc579	261	//____________________________________________________________________________
	262	TClonesArray* AliCaloTrackReader::GetAODMCParticles(Int_t input) const {
	263	//Return list of particles in AOD. Do it for the corresponding input event.
	264	if(fDataType == kAOD){
	265	//Normal input AOD
	266	if(input == 0) return (TClonesArray)((AliAODEvent)fInputEvent)->FindListObject("mcparticles");
	267	//Second input AOD
	268	else if(input == 1 && fSecondInputAODEvent) return (TClonesArray*) fSecondInputAODEvent->FindListObject("mcparticles");
	269	else {
	270	printf("AliCaloTrackReader::GetAODMCParticles() - wrong AOD input index? %d, or non existing tree? \n",input);
	271	return 0x0;
	272	}
	273	}
	274	else {
	275	printf("AliCaloTrackReader::GetAODMCParticles() - Input are not AODs\n");
	276	return 0x0;
	277	}
	278	}
	279
	280	//____________________________________________________________________________
	281	AliAODMCHeader* AliCaloTrackReader::GetAODMCHeader(Int_t input) const {
	282	//Return MC header in AOD. Do it for the corresponding input event.
	283	if(fDataType == kAOD){
	284	//Normal input AOD
	285	if(input == 0) return (AliAODMCHeader)((AliAODEvent)fInputEvent)->FindListObject("mcheader");
	286	//Second input AOD
	287	else if(input == 1) return (AliAODMCHeader*) fSecondInputAODEvent->FindListObject("mcheader");
	288	else {
	289	printf("AliCaloTrackReader::GetAODMCHeader() - wrong AOD input index, %d\n",input);
	290	return 0x0;
	291	}
	292	}
	293	else {
	294	printf("AliCaloTrackReader::GetAODMCHeader() - Input are not AODs\n");
	295	return 0x0;
	296	}
	297	}
	298
	299	//_______________________________________________________________
	300	void AliCaloTrackReader::Init()
	301	{
	302	//Init reader. Method to be called in AliAnaPartCorrMaker
	303
	304	//Get the file with second input events if the filename is given
	305	//Get the tree and connect the AODEvent. Only with AODs
f1f0bd84	306
591cc579	307	if(fReadStack && fReadAODMCParticles){
	308	printf("AliCaloTrackReader::Init() - Cannot access stack and mcparticles at the same time, change them \n");
	309	fReadStack = kFALSE;
	310	fReadAODMCParticles = kFALSE;
	311	}
	312
591cc579	313	if(fSecondInputFileName!=""){
	314	if(fDataType == kAOD){
	315	TFile * input2 = new TFile(fSecondInputFileName,"read");
ddcb0d54	316	printf("AliCaloTrackReader::Init() - Second input file opened: %s, size %d \n", input2->GetName(), (Int_t) input2->GetSize());
591cc579	317	fSecondInputAODTree = (TTree*) input2->Get("aodTree");
ddcb0d54	318	if(fSecondInputAODTree) printf("AliCaloTrackReader::Init() - Second input tree opened: %s, entries %d \n",
	319	fSecondInputAODTree->GetName(), (Int_t) fSecondInputAODTree->GetEntries());
	320	else{
	321	printf("AliCaloTrackReader::Init() - Second input tree not available, STOP \n");
	322	abort();
	323	}
591cc579	324	fSecondInputAODEvent = new AliAODEvent;
	325	fSecondInputAODEvent->ReadFromTree(fSecondInputAODTree);
	326	if(fSecondInputFirstEvent >= fSecondInputAODTree->GetEntriesFast()){
	327	printf("AliCaloTrackReader::Init() - Requested first event of second input %d, is larger than number of events %d, STOP\n",
	328	fSecondInputFirstEvent, (Int_t) fSecondInputAODTree->GetEntriesFast());
	329	abort();
	330	}
	331	}
	332	else printf("AliCaloTrackReader::Init() - Second input not added, reader is not AOD\n");
	333	}
	334
	335
	336	}
1c5acb87	337	//_______________________________________________________________
	338	void AliCaloTrackReader::InitParameters()
	339	{
	340
	341	//Initialize the parameters of the analysis.
591cc579	342	fDataType = kESD ;
1c5acb87	343	fCTSPtMin = 0.2 ;
29b2ceec	344	fEMCALPtMin = 0.2 ;
29b2ceec	345	fPHOSPtMin = 0.2 ;
1c5acb87	346
902aa95c	347	//Do not filter the detectors input by default.
	348	fFillEMCAL = kFALSE;
	349	fFillPHOS = kFALSE;
	350	fFillCTS = kFALSE;
1c5acb87	351	fFillEMCALCells = kFALSE;
591cc579	352	fFillPHOSCells = kFALSE;
1c5acb87	353
591cc579	354	fFidutialCut = new AliFidutialCut();
	355	fSecondInputFileName = "" ;
	356	fSecondInputFirstEvent = 0 ;
	357	fReadStack = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
	358	fReadAODMCParticles = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
08a064bc	359	fCleanOutputStdAOD = kFALSE; // Clean the standard clusters/tracks?
42dc8e7d	360	fDeltaAODFileName = "deltaAODPartCorr.root";
1c5acb87	361	}
	362
	363
	364	//________________________________________________________________
	365	void AliCaloTrackReader::Print(const Option_t * opt) const
	366	{
	367
	368	//Print some relevant parameters set for the analysis
	369	if(! opt)
	370	return;
	371
	372	printf("*** Print: %s %s ****\n", GetName(), GetTitle() ) ;
	373	printf("Data type : %d\n", fDataType) ;
	374	printf("CTS Min pT : %2.1f GeV/c\n", fCTSPtMin) ;
	375	printf("EMCAL Min pT : %2.1f GeV/c\n", fEMCALPtMin) ;
	376	printf("PHOS Min pT : %2.1f GeV/c\n", fPHOSPtMin) ;
	377	printf("Use CTS = %d\n", fFillCTS) ;
	378	printf("Use EMCAL = %d\n", fFillEMCAL) ;
	379	printf("Use PHOS = %d\n", fFillPHOS) ;
	380	printf("Use EMCAL Cells = %d\n", fFillEMCALCells) ;
	381	printf("Use PHOS Cells = %d\n", fFillPHOSCells) ;
591cc579	382	printf("Track status = %d\n", (Int_t) fTrackStatus) ;
	383	if(fComparePtHardAndJetPt)
	384	printf("Compare jet pt and pt hard to accept event, factor = %2.2f",fPtHardAndJetPtFactor);
	385
	386	if(fSecondInputFileName!="") {
	387	printf("Second Input File Name = %s\n", fSecondInputFileName.Data()) ;
	388	printf("Second Input First Event = %d\n", fSecondInputFirstEvent) ;
	389	}
	390
42dc8e7d	391	printf("Read Kine from, stack? %d, AOD ? %d \n", fReadStack, fReadAODMCParticles) ;
08a064bc	392	printf("Clean std AOD = %d\n", fCleanOutputStdAOD) ;
42dc8e7d	393	printf("Delta AOD File Name = %s\n", fDeltaAODFileName.Data()) ;
42dc8e7d	394
1c5acb87	395	printf(" \n") ;
	396	}
	397
	398	//___________________________________________________
29b2ceec	399	Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry, const char * currentFileName) {
6639984f	400	//Fill the event counter and input lists that are needed, called by the analysis maker.
902aa95c	401
6639984f	402	fEventNumber = iEntry;
a79a2424	403	fCurrentFileName = TString(currentFileName);
a79a2424	404
902aa95c	405	if(fOutputEvent && (fDataType != kAOD) && ((fOutputEvent->GetCaloClusters())->GetEntriesFast()!=0 \|\|(fOutputEvent->GetTracks())->GetEntriesFast()!=0)){
29b2ceec	406	printf("AliCaloTrackReader::AODCaloClusters or AODTracks already filled by the filter, do not use the ESD reader, use the AOD reader, STOP\n");
477d6cee	407	abort();
477d6cee	408	}
902aa95c	409
29b2ceec	410	//In case of analysis of events with jets, skip those with jet pt > 5 pt hard
29b2ceec	411	if(fComparePtHardAndJetPt && GetStack()) {
591cc579	412	if(!ComparePtHardAndJetPt()) return kFALSE ;
29b2ceec	413	}
591cc579	414
	415	//In case of mixing events with other AOD file
	416	if(fDataType == kAOD && fSecondInputAODTree){
	417
	418	if(fDebug > 1)
	419	printf("AliCaloTrackReader::FillInputEvent() - Get event %d from second input AOD file \n", iEntry+fSecondInputFirstEvent);
	420	if(fSecondInputAODTree->GetEntriesFast() <= iEntry+fSecondInputFirstEvent) {
	421	if(fSecondInputAODTree->GetEntriesFast() == iEntry+fSecondInputFirstEvent)
	422	printf("AliCaloTrackReader::FillInputEvent() - Skip events from event %d, no more events in second AOD file \n", iEntry);
	423	return kFALSE;
	424	}
902aa95c	425	printf("Reader 4 \n");
902aa95c	426
591cc579	427	//Get the Event
	428	Int_t nbytes = fSecondInputAODTree->GetEvent(iEntry+fSecondInputFirstEvent);
	429	if ( nbytes == 0 ) {//If nothing in AOD
	430	printf("AliCaloTrackReader::FillInputEvent() - Nothing in Second AOD input, STOP\n");
	431	abort() ;
	432	}
	433
	434	}
902aa95c	435
477d6cee	436	if(fFillCTS) FillInputCTS();
1c5acb87	437	if(fFillEMCAL) FillInputEMCAL();
477d6cee	438	if(fFillPHOS) FillInputPHOS();
1c5acb87	439	if(fFillEMCALCells) FillInputEMCALCells();
477d6cee	440	if(fFillPHOSCells) FillInputPHOSCells();
1c5acb87	441
29b2ceec	442	return kTRUE ;
1c5acb87	443	}
	444
	445	//__________________________________________________
	446	void AliCaloTrackReader::ResetLists() {
	447	// Reset lists, called by the analysis maker
	448
477d6cee	449	if(fAODCTS) fAODCTS -> Clear();
1c5acb87	450	if(fAODEMCAL) fAODEMCAL -> Clear();
477d6cee	451	if(fAODPHOS) fAODPHOS -> Clear();
1c5acb87	452	if(fEMCALCells) fEMCALCells -> Clear();
477d6cee	453	if(fPHOSCells) fPHOSCells -> Clear();
902aa95c	454	if(fCleanOutputStdAOD && fOutputEvent ){
42dc8e7d	455	//Only keep copied tracks and clusters if requested
	456	fOutputEvent->GetTracks() ->Clear();
	457	fOutputEvent->GetCaloClusters()->Clear();
	458	}
1c5acb87	459	}