[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliConversionSelection.cxx

#include "AliConversionSelection.h"
#include "AliAODHeader.h"
#include "AliVVZERO.h"
#include "AliMultiplicity.h"
#include <iostream>


// Author Daniel Lohner (Daniel.Lohner@cern.ch)

using namespace std;

ClassImp(AliConversionSelection)


//________________________________________________________________________
AliConversionSelection::AliConversionSelection(AliConvEventCuts *evtCut, AliConversionPhotonCuts *convCut, AliConversionMesonCuts *mesonCut) : TObject(),
	fInputEvent(NULL),
	fMCEvent(NULL),
	fEventCut(evtCut),
	fConversionCut(convCut),
	fMesonCut(mesonCut),
	fESDTrackCuts(NULL),
	fGoodGammas(NULL),
	fPi0Candidates(NULL),
	fBGPi0s(NULL),
	fRandomizer(NULL),
	fBGHandler(NULL),
	fCurrentEventNumber(-1),
	fIsOwner(kFALSE)
{
	// Default Values
	fInvMassRange[0]=0.05;
	fInvMassRange[1]=0.3;
}

//________________________________________________________________________
AliConversionSelection::AliConversionSelection(TString evtCut, TString convCut, TString mesonCut) : TObject(),
	fInputEvent(NULL),
	fMCEvent(NULL),
	fEventCut(NULL),
	fConversionCut(NULL),
	fMesonCut(NULL),
	fESDTrackCuts(NULL),
	fGoodGammas(NULL),
	fPi0Candidates(NULL),
	fBGPi0s(NULL),
	fRandomizer(NULL),
	fBGHandler(NULL),
	fCurrentEventNumber(-1),
	fIsOwner(kTRUE)
{
	// Default Values
	fInvMassRange[0]=0.05;
	fInvMassRange[1]=0.3;

	fEventCut = new AliConvEventCuts();
	fEventCut -> InitializeCutsFromCutString(evtCut.Data());
	fConversionCut = new AliConversionPhotonCuts();
	fConversionCut -> InitializeCutsFromCutString(convCut.Data());
	fMesonCut = new AliConversionMesonCuts();
	fMesonCut -> InitializeCutsFromCutString(mesonCut.Data());

}


//________________________________________________________________________
AliConversionSelection::AliConversionSelection(const AliConversionSelection &ref) : TObject(ref),
	fInputEvent(NULL),
	fMCEvent(NULL),
	fEventCut(NULL),
	fConversionCut(NULL),
	fMesonCut(NULL),
	fESDTrackCuts(NULL),
	fGoodGammas(NULL),
	fPi0Candidates(NULL),
	fBGPi0s(NULL),
	fRandomizer(NULL),
	fBGHandler(NULL),
	fCurrentEventNumber(-1),
	fIsOwner(kTRUE)
{
	// Copy Constructor
	fEventCut=new AliConvEventCuts(*ref.fEventCut);
	fConversionCut=new AliConversionPhotonCuts(*ref.fConversionCut);
	fMesonCut=new AliConversionMesonCuts(*ref.fMesonCut);

	fInvMassRange[0]=ref.fInvMassRange[0];
	fInvMassRange[1]=ref.fInvMassRange[1];
}

//________________________________________________________________________
AliConversionSelection::~AliConversionSelection(){

	if(fBGHandler){
		delete fBGHandler;
		fBGHandler=NULL;
	}
	if(fRandomizer){
		delete fRandomizer;
		fRandomizer=NULL;
	}
	if(fPi0Candidates){
		delete fPi0Candidates;
		fPi0Candidates=NULL;
	}
	if(fBGPi0s){
		delete fBGPi0s;
		fBGPi0s=NULL;
	}
	if(fESDTrackCuts){
		delete fESDTrackCuts;
		fESDTrackCuts=NULL;
	}
	if(fIsOwner){
		if(fEventCut){
			delete fEventCut;
			fEventCut=NULL;
		}
		if(fConversionCut){
			delete fConversionCut;
			fConversionCut=NULL;
		}
		if(fMesonCut){
			delete fMesonCut;
			fMesonCut=NULL;
		}
	}
}

//________________________________________________________________________
Bool_t AliConversionSelection::ProcessEvent(TClonesArray *photons,AliVEvent *inputEvent,AliMCEvent *mcEvent){
	fInputEvent=inputEvent;
	fMCEvent=mcEvent;

	// Protection
	Int_t eventnumber=GetEventNumber(inputEvent);
	if(eventnumber==fCurrentEventNumber){
		AliWarning("Event already analyzed! Return.");
		return kFALSE;
	}
	else{
		fCurrentEventNumber=eventnumber;
	}

	// Initialize and Reset Arrays
	if(fGoodGammas == NULL){
		fGoodGammas=new TObjArray(30);
	}
	fGoodGammas->Clear();

	if(fPi0Candidates == NULL){
		fPi0Candidates = new TClonesArray("AliAODConversionMother",100);
	}
	fPi0Candidates->Delete();

	if(fBGPi0s == NULL){
		fBGPi0s = new TClonesArray("AliAODConversionMother",100);
	}
	fBGPi0s->Delete();


	if(!photons||!fInputEvent)return kFALSE;
	
	if(!fEventCut->EventIsSelected(fInputEvent,fMCEvent))return kFALSE;

	// Select photons
	for(Int_t i = 0; i < photons->GetEntriesFast(); i++) {
		AliAODConversionPhoton* gamma =dynamic_cast<AliAODConversionPhoton*>(photons->At(i));
		if(!gamma) continue;
		if(!fConversionCut->PhotonIsSelected(gamma,fInputEvent))continue;
		fGoodGammas->Add(gamma);
	}

	// Do MC Smearing
	Double_t *fUnsmearedPx=NULL;
	Double_t *fUnsmearedPy=NULL;
	Double_t *fUnsmearedPz=NULL;
	Double_t *fUnsmearedE=NULL;

	if(fMesonCut->UseMCPSmearing() && fMCEvent){
		fUnsmearedPx = new Double_t[fGoodGammas->GetEntries()]; // Store unsmeared Momenta
		fUnsmearedPy = new Double_t[fGoodGammas->GetEntries()];
		fUnsmearedPz = new Double_t[fGoodGammas->GetEntries()];
		fUnsmearedE =  new Double_t[fGoodGammas->GetEntries()];

		for(Int_t gamma=0;gamma<fGoodGammas->GetEntries();gamma++){ // Smear the AODPhotons in MC
			fUnsmearedPx[gamma] = ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->Px();
			fUnsmearedPy[gamma] = ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->Py();
			fUnsmearedPz[gamma] = ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->Pz();
			fUnsmearedE[gamma] =  ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->E();
			fMesonCut->SmearParticle(dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(gamma)));
		}
	}

	// Reconstruct Pi0 and BG
	CalculatePi0Candidates();
	CalculateBackground();
	if(fBGHandler)fBGHandler->AddEvent(fGoodGammas,fInputEvent);

	// Undo MC Smearing
	if(fMesonCut->UseMCPSmearing() && fMCEvent){
		for(Int_t gamma=0;gamma<fGoodGammas->GetEntries();gamma++){ // Smear the AODPhotons in MC
			((AliAODConversionPhoton*)fGoodGammas->At(gamma))->SetPx(fUnsmearedPx[gamma]); // Reset Unsmeared Momenta
			((AliAODConversionPhoton*)fGoodGammas->At(gamma))->SetPy(fUnsmearedPy[gamma]);
			((AliAODConversionPhoton*)fGoodGammas->At(gamma))->SetPz(fUnsmearedPz[gamma]);
			((AliAODConversionPhoton*)fGoodGammas->At(gamma))->SetE(fUnsmearedE[gamma]);
		}
		delete[] fUnsmearedPx; fUnsmearedPx = 0x0;
		delete[] fUnsmearedPy; fUnsmearedPy = 0x0;
		delete[] fUnsmearedPz; fUnsmearedPz = 0x0;
		delete[] fUnsmearedE;  fUnsmearedE  = 0x0;
	}
	return kTRUE;
}

//________________________________________________________________________
AliAODConversionMother* AliConversionSelection::GetPi0(Int_t index){

	if(index>=0&&index<GetNumberOfPi0s()){
		return dynamic_cast<AliAODConversionMother*>(fPi0Candidates->At(index));
	}
	return NULL;
}

//________________________________________________________________________
AliAODConversionMother* AliConversionSelection::GetBG(Int_t index){

	if(index>=0&&index<GetNumberOfBGs()){
		return dynamic_cast<AliAODConversionMother*>(fBGPi0s->At(index));
	}
	return NULL;
}

//________________________________________________________________________
AliAODConversionPhoton* AliConversionSelection::GetPhoton(Int_t index){

	if(index>=0&&index<GetNumberOfPhotons()){
		return dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(index));
	}
	return NULL;
}

//________________________________________________________________________
void AliConversionSelection::CalculatePi0Candidates(){
	// Conversion Gammas
	if(fGoodGammas->GetEntriesFast()>1){
		for(Int_t firstGammaIndex=0;firstGammaIndex<fGoodGammas->GetEntriesFast()-1;firstGammaIndex++){
			AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(firstGammaIndex));
			if (gamma0==NULL) continue;
			// Combine Photons

			for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fGoodGammas->GetEntriesFast();secondGammaIndex++){

				AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(secondGammaIndex));
				if (gamma1==NULL) continue;
				//Check for same Electron ID
				if(gamma0->GetTrackLabelPositive()==gamma1->GetTrackLabelPositive()||gamma0->GetTrackLabelNegative()==gamma1->GetTrackLabelNegative()
				||gamma0->GetTrackLabelNegative()==gamma1->GetTrackLabelPositive()||gamma0->GetTrackLabelPositive()==gamma1->GetTrackLabelNegative())continue;

				AliAODConversionMother pi0cand(gamma0,gamma1);
				pi0cand.SetLabels(firstGammaIndex,secondGammaIndex);

				// Set MC Label

				if(fMCEvent){

					TParticle *mcgam0=gamma0->GetMCParticle(fMCEvent->Stack());
					TParticle *mcgam1=gamma1->GetMCParticle(fMCEvent->Stack());

					if(mcgam0&&mcgam1){
					// Have same Mother?

					if(mcgam0->GetMother(0)==mcgam1->GetMother(0)){

						pi0cand.SetMCLabel(mcgam0->GetMother(0));
					}
					}
				}

				if((fMesonCut->MesonIsSelected(&pi0cand))){
					if(MesonInMassWindow(&pi0cand)){

					// Add Pi0 to Stack
					new((*fPi0Candidates)[fPi0Candidates->GetEntriesFast()]) AliAODConversionMother(pi0cand);
					}
				}
			}
		}
	}
}

//________________________________________________________________________
Bool_t AliConversionSelection::MesonInMassWindow(AliAODConversionMother *pi0cand)
{
	if (pi0cand->M() > fInvMassRange[0] && pi0cand->M() < fInvMassRange[1] ){
		return kTRUE;
	}
	return kFALSE;
}

//________________________________________________________________________
void AliConversionSelection::RotateParticle(AliAODConversionPhoton *gamma,Int_t nDegreesPMBackground){

	if(!fRandomizer){
		fRandomizer=new TRandom3();
		fRandomizer->SetSeed(0);
	}
	Double_t nRadiansPM = nDegreesPMBackground*TMath::Pi()/180;
	Double_t rotationValue = fRandomizer->Rndm()*2*nRadiansPM + TMath::Pi()-nRadiansPM;
	gamma->RotateZ(rotationValue);
}

//________________________________________________________________________

void AliConversionSelection::CalculateBackground(){

	//Rotation Method
	if(fMesonCut->UseRotationMethod()){
		// Correct for the number of rotations
		// BG is for rotation the same, except for factor NRotations
		Double_t weight=1./Double_t(fMesonCut->GetNumberOfBGEvents());
		for(Int_t firstGammaIndex=0;firstGammaIndex<fGoodGammas->GetEntriesFast();firstGammaIndex++){
			AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(firstGammaIndex));
			if (gamma0 ==NULL) continue;
			for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fGoodGammas->GetEntriesFast();secondGammaIndex++){
				AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(secondGammaIndex));
				if (gamma1==NULL) continue;
				if(!fConversionCut->PhotonIsSelected(gamma1,fInputEvent))continue;
				for(Int_t nRandom=0;nRandom<fMesonCut->GetNumberOfBGEvents();nRandom++){
					RotateParticle(gamma1,fMesonCut->NDegreesRotation());
					AliAODConversionMother BGcandidate(gamma0,gamma1);
					if(fMesonCut->MesonIsSelected(&BGcandidate,kFALSE)){
						if(MesonInMassWindow(&BGcandidate)){
							new((*fBGPi0s)[fBGPi0s->GetEntriesFast()]) AliAODConversionMother(BGcandidate);
							dynamic_cast<AliAODConversionMother*>(fBGPi0s->At(fBGPi0s->GetEntriesFast()-1))->SetWeight(weight);
						}
					}
				}
			}
		}
	} else {
		// Do Event Mixing
		if(fBGHandler==NULL){
			fBGHandler=new AliConversionAODBGHandlerRP(fEventCut->IsHeavyIon(),fMesonCut->UseTrackMultiplicity());
		}

		for(Int_t nEventsInBG=0;nEventsInBG <fBGHandler->GetNBGEvents(fGoodGammas,fInputEvent);nEventsInBG++){
			AliGammaConversionPhotonVector *previousEventGammas = fBGHandler->GetBGGoodGammas(fGoodGammas,fInputEvent,nEventsInBG);
			if(previousEventGammas){
				// test weighted background
				Double_t weight=1.0;
				// Correct for the number of eventmixing:
				// N gammas -> (N-1) + (N-2) +(N-3) ...+ (N-(N-1))  using sum formula sum(i)=N*(N-1)/2  -> N*(N-1)/2
				// real combinations (since you cannot combine a photon with its own)
				// but BG leads to N_{a}*N_{b} combinations
				weight*=0.5*(Double_t(fGoodGammas->GetEntriesFast()-1))/Double_t(previousEventGammas->size());
				for(Int_t iCurrent=0;iCurrent<fGoodGammas->GetEntriesFast();iCurrent++){
					AliAODConversionPhoton *gamma0 = (AliAODConversionPhoton*)(fGoodGammas->At(iCurrent));
					for(UInt_t iPrevious=0;iPrevious<previousEventGammas->size();iPrevious++){
						AliAODConversionPhoton *gamma1 = (AliAODConversionPhoton*)(previousEventGammas->at(iPrevious));
						AliAODConversionMother BGcandidate(gamma0,gamma1);
						if(fMesonCut->MesonIsSelected(&BGcandidate,kFALSE)){
							if(MesonInMassWindow(&BGcandidate)){
								new((*fBGPi0s)[fBGPi0s->GetEntriesFast()]) AliAODConversionMother(BGcandidate);
								dynamic_cast<AliAODConversionMother*>(fBGPi0s->At(fBGPi0s->GetEntriesFast()-1))->SetWeight(weight);
							}
						}
					}
				}
			}
		}
    }
}
//________________________________________________________________________
Double_t AliConversionSelection::GetMultiplicity(AliVEvent *inputEvent){

    switch(fEventCut->GetMultiplicityMethod()){
		case 0:
			return Double_t(GetNumberOfPhotons());
		case 1:
			return Double_t(GetNumberOfChargedTracks(inputEvent));
		case 2:
			return GetVZEROMult(inputEvent);
		case 3:
			return GetSPDMult(inputEvent);
		case 9:
			return 1; // if mult is used as a weight, this number can be used to switch off weighting
		default:
		return 0;
    }
}

//________________________________________________________________________
Int_t AliConversionSelection::GetNumberOfChargedTracks(AliVEvent *inputEvent){

	Int_t ntracks = 0;

	AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(inputEvent);
	if(esdEvent) {
		if(!fESDTrackCuts){
			fESDTrackCuts= AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
			fESDTrackCuts->SetMaxDCAToVertexZ(2);
			Double_t etamax=fConversionCut->GetEtaCut();
			fESDTrackCuts->SetEtaRange(-etamax, etamax);
			fESDTrackCuts->SetPtRange(0.15);
		}
		for(Int_t iTracks = 0; iTracks < inputEvent->GetNumberOfTracks(); iTracks++){
			AliESDtrack* currentTrack = esdEvent->GetTrack(iTracks);
			if(!currentTrack) continue;
			if(fESDTrackCuts->AcceptTrack(currentTrack))ntracks++;
		}
	} else {
		for(Int_t ii=0; ii<inputEvent->GetNumberOfTracks(); ii++) {
			AliVTrack * track = dynamic_cast<AliVTrack*>(inputEvent->GetTrack(ii));
			if (track==NULL) continue;
			if(TMath::Abs(track->Eta())>fConversionCut->GetEtaCut())continue;
			if(track)ntracks++;
		}
	}
	return ntracks;
}

//________________________________________________________________________
Double_t AliConversionSelection::GetVZEROMult(AliVEvent *inputEvent){

	AliVVZERO *vzero=inputEvent->GetVZEROData();
	Double_t multV0A=vzero->GetMTotV0A();
	Double_t multV0C=vzero->GetMTotV0C();
	Double_t mult=multV0A+multV0C;

	return mult;
}

//________________________________________________________________________
Double_t AliConversionSelection::GetSPDMult(AliVEvent *inputEvent){

	AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(inputEvent);
	if(esdEvent) {
		const AliMultiplicity *esdmult=esdEvent->GetMultiplicity();
		return esdmult->GetNumberOfITSClusters(1);
	} else {
		// AOD implementation
		AliAODHeader *header=(AliAODHeader*)inputEvent->GetHeader();
		return header->GetNumberOfITSClusters(1);
	}
	return 0;
}

//________________________________________________________________________
Int_t AliConversionSelection::GetEventNumber(AliVEvent *inputEvent){

	AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(inputEvent);
	if(esdEvent) {
		return esdEvent->GetEventNumberInFile();
	} else{
		AliAODHeader *header=(AliAODHeader*)inputEvent->GetHeader();
		return header->GetEventNumberESDFile();
	}
	return 0;
}

//________________________________________________________________________
TString AliConversionSelection::GetCutString(){
	TString a= Form("%s_%s_%s",fEventCut->GetCutNumber().Data(), fConversionCut->GetCutNumber().Data(),fMesonCut->GetCutNumber().Data());
	return a;
}
Commit	Line	Data
92efd725	1	#include "AliConversionSelection.h"
	2	#include "AliAODHeader.h"
	3	#include "AliVVZERO.h"
	4	#include "AliMultiplicity.h"
	5	#include <iostream>
	6
	7
	8	// Author Daniel Lohner (Daniel.Lohner@cern.ch)
	9
	10	using namespace std;
	11
	12	ClassImp(AliConversionSelection)
	13
	14
	15	//________________________________________________________________________
344100c4	16	AliConversionSelection::AliConversionSelection(AliConvEventCuts evtCut, AliConversionPhotonCuts convCut, AliConversionMesonCuts *mesonCut) : TObject(),
	17	fInputEvent(NULL),
	18	fMCEvent(NULL),
	19	fEventCut(evtCut),
	20	fConversionCut(convCut),
	21	fMesonCut(mesonCut),
	22	fESDTrackCuts(NULL),
	23	fGoodGammas(NULL),
	24	fPi0Candidates(NULL),
	25	fBGPi0s(NULL),
	26	fRandomizer(NULL),
	27	fBGHandler(NULL),
	28	fCurrentEventNumber(-1),
	29	fIsOwner(kFALSE)
92efd725	30	{
344100c4	31	// Default Values
	32	fInvMassRange[0]=0.05;
	33	fInvMassRange[1]=0.3;
92efd725	34	}
e5b6e8a6	35
e5b6e8a6	36	//________________________________________________________________________
344100c4	37	AliConversionSelection::AliConversionSelection(TString evtCut, TString convCut, TString mesonCut) : TObject(),
	38	fInputEvent(NULL),
	39	fMCEvent(NULL),
	40	fEventCut(NULL),
	41	fConversionCut(NULL),
	42	fMesonCut(NULL),
	43	fESDTrackCuts(NULL),
	44	fGoodGammas(NULL),
	45	fPi0Candidates(NULL),
	46	fBGPi0s(NULL),
	47	fRandomizer(NULL),
	48	fBGHandler(NULL),
	49	fCurrentEventNumber(-1),
	50	fIsOwner(kTRUE)
e5b6e8a6	51	{
344100c4	52	// Default Values
	53	fInvMassRange[0]=0.05;
	54	fInvMassRange[1]=0.3;
e5b6e8a6	55
344100c4	56	fEventCut = new AliConvEventCuts();
	57	fEventCut -> InitializeCutsFromCutString(evtCut.Data());
	58	fConversionCut = new AliConversionPhotonCuts();
	59	fConversionCut -> InitializeCutsFromCutString(convCut.Data());
	60	fMesonCut = new AliConversionMesonCuts();
	61	fMesonCut -> InitializeCutsFromCutString(mesonCut.Data());
e5b6e8a6	62
	63	}
	64
1d9e6011	65
	66	//________________________________________________________________________
	67	AliConversionSelection::AliConversionSelection(const AliConversionSelection &ref) : TObject(ref),
344100c4	68	fInputEvent(NULL),
	69	fMCEvent(NULL),
	70	fEventCut(NULL),
	71	fConversionCut(NULL),
	72	fMesonCut(NULL),
	73	fESDTrackCuts(NULL),
	74	fGoodGammas(NULL),
	75	fPi0Candidates(NULL),
	76	fBGPi0s(NULL),
	77	fRandomizer(NULL),
	78	fBGHandler(NULL),
	79	fCurrentEventNumber(-1),
	80	fIsOwner(kTRUE)
1d9e6011	81	{
344100c4	82	// Copy Constructor
	83	fEventCut=new AliConvEventCuts(*ref.fEventCut);
	84	fConversionCut=new AliConversionPhotonCuts(*ref.fConversionCut);
	85	fMesonCut=new AliConversionMesonCuts(*ref.fMesonCut);
1d9e6011	86
344100c4	87	fInvMassRange[0]=ref.fInvMassRange[0];
344100c4	88	fInvMassRange[1]=ref.fInvMassRange[1];
1d9e6011	89	}
1d9e6011	90
92efd725	91	//________________________________________________________________________
	92	AliConversionSelection::~AliConversionSelection(){
	93
344100c4	94	if(fBGHandler){
	95	delete fBGHandler;
	96	fBGHandler=NULL;
e5b6e8a6	97	}
344100c4	98	if(fRandomizer){
	99	delete fRandomizer;
	100	fRandomizer=NULL;
	101	}
	102	if(fPi0Candidates){
	103	delete fPi0Candidates;
	104	fPi0Candidates=NULL;
	105	}
	106	if(fBGPi0s){
	107	delete fBGPi0s;
	108	fBGPi0s=NULL;
	109	}
	110	if(fESDTrackCuts){
	111	delete fESDTrackCuts;
	112	fESDTrackCuts=NULL;
	113	}
	114	if(fIsOwner){
	115	if(fEventCut){
	116	delete fEventCut;
	117	fEventCut=NULL;
	118	}
	119	if(fConversionCut){
	120	delete fConversionCut;
	121	fConversionCut=NULL;
	122	}
	123	if(fMesonCut){
	124	delete fMesonCut;
	125	fMesonCut=NULL;
	126	}
e5b6e8a6	127	}
92efd725	128	}
	129
	130	//________________________________________________________________________
	131	Bool_t AliConversionSelection::ProcessEvent(TClonesArray photons,AliVEvent inputEvent,AliMCEvent *mcEvent){
344100c4	132	fInputEvent=inputEvent;
	133	fMCEvent=mcEvent;
	134
	135	// Protection
	136	Int_t eventnumber=GetEventNumber(inputEvent);
	137	if(eventnumber==fCurrentEventNumber){
	138	AliWarning("Event already analyzed! Return.");
	139	return kFALSE;
	140	}
	141	else{
	142	fCurrentEventNumber=eventnumber;
	143	}
92efd725	144
344100c4	145	// Initialize and Reset Arrays
	146	if(fGoodGammas == NULL){
	147	fGoodGammas=new TObjArray(30);
	148	}
	149	fGoodGammas->Clear();
92efd725	150
344100c4	151	if(fPi0Candidates == NULL){
	152	fPi0Candidates = new TClonesArray("AliAODConversionMother",100);
	153	}
	154	fPi0Candidates->Delete();
92efd725	155
344100c4	156	if(fBGPi0s == NULL){
	157	fBGPi0s = new TClonesArray("AliAODConversionMother",100);
	158	}
	159	fBGPi0s->Delete();
92efd725	160
92efd725	161
344100c4	162	if(!photons\|\|!fInputEvent)return kFALSE;
	163
	164	if(!fEventCut->EventIsSelected(fInputEvent,fMCEvent))return kFALSE;
92efd725	165
344100c4	166	// Select photons
	167	for(Int_t i = 0; i < photons->GetEntriesFast(); i++) {
	168	AliAODConversionPhoton* gamma =dynamic_cast<AliAODConversionPhoton*>(photons->At(i));
	169	if(!gamma) continue;
	170	if(!fConversionCut->PhotonIsSelected(gamma,fInputEvent))continue;
	171	fGoodGammas->Add(gamma);
92efd725	172	}
92efd725	173
344100c4	174	// Do MC Smearing
	175	Double_t *fUnsmearedPx=NULL;
	176	Double_t *fUnsmearedPy=NULL;
	177	Double_t *fUnsmearedPz=NULL;
	178	Double_t *fUnsmearedE=NULL;
	179
	180	if(fMesonCut->UseMCPSmearing() && fMCEvent){
	181	fUnsmearedPx = new Double_t[fGoodGammas->GetEntries()]; // Store unsmeared Momenta
	182	fUnsmearedPy = new Double_t[fGoodGammas->GetEntries()];
	183	fUnsmearedPz = new Double_t[fGoodGammas->GetEntries()];
	184	fUnsmearedE = new Double_t[fGoodGammas->GetEntries()];
	185
	186	for(Int_t gamma=0;gamma<fGoodGammas->GetEntries();gamma++){ // Smear the AODPhotons in MC
	187	fUnsmearedPx[gamma] = ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->Px();
	188	fUnsmearedPy[gamma] = ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->Py();
	189	fUnsmearedPz[gamma] = ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->Pz();
	190	fUnsmearedE[gamma] = ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->E();
	191	fMesonCut->SmearParticle(dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(gamma)));
	192	}
92efd725	193	}
92efd725	194
344100c4	195	// Reconstruct Pi0 and BG
	196	CalculatePi0Candidates();
	197	CalculateBackground();
	198	if(fBGHandler)fBGHandler->AddEvent(fGoodGammas,fInputEvent);
	199
	200	// Undo MC Smearing
	201	if(fMesonCut->UseMCPSmearing() && fMCEvent){
	202	for(Int_t gamma=0;gamma<fGoodGammas->GetEntries();gamma++){ // Smear the AODPhotons in MC
	203	((AliAODConversionPhoton*)fGoodGammas->At(gamma))->SetPx(fUnsmearedPx[gamma]); // Reset Unsmeared Momenta
	204	((AliAODConversionPhoton*)fGoodGammas->At(gamma))->SetPy(fUnsmearedPy[gamma]);
	205	((AliAODConversionPhoton*)fGoodGammas->At(gamma))->SetPz(fUnsmearedPz[gamma]);
	206	((AliAODConversionPhoton*)fGoodGammas->At(gamma))->SetE(fUnsmearedE[gamma]);
	207	}
	208	delete[] fUnsmearedPx; fUnsmearedPx = 0x0;
	209	delete[] fUnsmearedPy; fUnsmearedPy = 0x0;
	210	delete[] fUnsmearedPz; fUnsmearedPz = 0x0;
	211	delete[] fUnsmearedE; fUnsmearedE = 0x0;
	212	}
	213	return kTRUE;
92efd725	214	}
	215
	216	//________________________________________________________________________
	217	AliAODConversionMother* AliConversionSelection::GetPi0(Int_t index){
	218
344100c4	219	if(index>=0&&index<GetNumberOfPi0s()){
	220	return dynamic_cast<AliAODConversionMother*>(fPi0Candidates->At(index));
	221	}
	222	return NULL;
92efd725	223	}
	224
	225	//________________________________________________________________________
	226	AliAODConversionMother* AliConversionSelection::GetBG(Int_t index){
	227
344100c4	228	if(index>=0&&index<GetNumberOfBGs()){
	229	return dynamic_cast<AliAODConversionMother*>(fBGPi0s->At(index));
	230	}
	231	return NULL;
92efd725	232	}
	233
	234	//________________________________________________________________________
	235	AliAODConversionPhoton* AliConversionSelection::GetPhoton(Int_t index){
	236
344100c4	237	if(index>=0&&index<GetNumberOfPhotons()){
	238	return dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(index));
	239	}
	240	return NULL;
92efd725	241	}
	242
	243	//________________________________________________________________________
	244	void AliConversionSelection::CalculatePi0Candidates(){
344100c4	245	// Conversion Gammas
	246	if(fGoodGammas->GetEntriesFast()>1){
	247	for(Int_t firstGammaIndex=0;firstGammaIndex<fGoodGammas->GetEntriesFast()-1;firstGammaIndex++){
	248	AliAODConversionPhoton gamma0=dynamic_cast<AliAODConversionPhoton>(fGoodGammas->At(firstGammaIndex));
	249	if (gamma0==NULL) continue;
	250	// Combine Photons
92efd725	251
344100c4	252	for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fGoodGammas->GetEntriesFast();secondGammaIndex++){
92efd725	253
344100c4	254	AliAODConversionPhoton gamma1=dynamic_cast<AliAODConversionPhoton>(fGoodGammas->At(secondGammaIndex));
	255	if (gamma1==NULL) continue;
	256	//Check for same Electron ID
	257	if(gamma0->GetTrackLabelPositive()==gamma1->GetTrackLabelPositive()\|\|gamma0->GetTrackLabelNegative()==gamma1->GetTrackLabelNegative()
	258	\|\|gamma0->GetTrackLabelNegative()==gamma1->GetTrackLabelPositive()\|\|gamma0->GetTrackLabelPositive()==gamma1->GetTrackLabelNegative())continue;
92efd725	259
344100c4	260	AliAODConversionMother pi0cand(gamma0,gamma1);
344100c4	261	pi0cand.SetLabels(firstGammaIndex,secondGammaIndex);
92efd725	262
344100c4	263	// Set MC Label
92efd725	264
344100c4	265	if(fMCEvent){
92efd725	266
344100c4	267	TParticle *mcgam0=gamma0->GetMCParticle(fMCEvent->Stack());
344100c4	268	TParticle *mcgam1=gamma1->GetMCParticle(fMCEvent->Stack());
92efd725	269
344100c4	270	if(mcgam0&&mcgam1){
344100c4	271	// Have same Mother?
92efd725	272
344100c4	273	if(mcgam0->GetMother(0)==mcgam1->GetMother(0)){
92efd725	274
344100c4	275	pi0cand.SetMCLabel(mcgam0->GetMother(0));
	276	}
	277	}
	278	}
92efd725	279
344100c4	280	if((fMesonCut->MesonIsSelected(&pi0cand))){
344100c4	281	if(MesonInMassWindow(&pi0cand)){
92efd725	282
344100c4	283	// Add Pi0 to Stack
	284	new((*fPi0Candidates)[fPi0Candidates->GetEntriesFast()]) AliAODConversionMother(pi0cand);
	285	}
	286	}
92efd725	287	}
92efd725	288	}
92efd725	289	}
92efd725	290	}
	291
	292	//________________________________________________________________________
	293	Bool_t AliConversionSelection::MesonInMassWindow(AliAODConversionMother *pi0cand)
	294	{
344100c4	295	if (pi0cand->M() > fInvMassRange[0] && pi0cand->M() < fInvMassRange[1] ){
	296	return kTRUE;
	297	}
	298	return kFALSE;
92efd725	299	}
	300
	301	//________________________________________________________________________
	302	void AliConversionSelection::RotateParticle(AliAODConversionPhoton *gamma,Int_t nDegreesPMBackground){
	303
344100c4	304	if(!fRandomizer){
	305	fRandomizer=new TRandom3();
	306	fRandomizer->SetSeed(0);
	307	}
	308	Double_t nRadiansPM = nDegreesPMBackground*TMath::Pi()/180;
	309	Double_t rotationValue = fRandomizer->Rndm()2nRadiansPM + TMath::Pi()-nRadiansPM;
	310	gamma->RotateZ(rotationValue);
92efd725	311	}
	312
	313	//________________________________________________________________________
	314
	315	void AliConversionSelection::CalculateBackground(){
	316
344100c4	317	//Rotation Method
	318	if(fMesonCut->UseRotationMethod()){
	319	// Correct for the number of rotations
	320	// BG is for rotation the same, except for factor NRotations
	321	Double_t weight=1./Double_t(fMesonCut->GetNumberOfBGEvents());
	322	for(Int_t firstGammaIndex=0;firstGammaIndex<fGoodGammas->GetEntriesFast();firstGammaIndex++){
	323	AliAODConversionPhoton gamma0=dynamic_cast<AliAODConversionPhoton>(fGoodGammas->At(firstGammaIndex));
	324	if (gamma0 ==NULL) continue;
	325	for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fGoodGammas->GetEntriesFast();secondGammaIndex++){
	326	AliAODConversionPhoton gamma1=dynamic_cast<AliAODConversionPhoton>(fGoodGammas->At(secondGammaIndex));
	327	if (gamma1==NULL) continue;
	328	if(!fConversionCut->PhotonIsSelected(gamma1,fInputEvent))continue;
	329	for(Int_t nRandom=0;nRandom<fMesonCut->GetNumberOfBGEvents();nRandom++){
	330	RotateParticle(gamma1,fMesonCut->NDegreesRotation());
	331	AliAODConversionMother BGcandidate(gamma0,gamma1);
	332	if(fMesonCut->MesonIsSelected(&BGcandidate,kFALSE)){
	333	if(MesonInMassWindow(&BGcandidate)){
	334	new((*fBGPi0s)[fBGPi0s->GetEntriesFast()]) AliAODConversionMother(BGcandidate);
	335	dynamic_cast<AliAODConversionMother*>(fBGPi0s->At(fBGPi0s->GetEntriesFast()-1))->SetWeight(weight);
	336	}
	337	}
	338	}
92efd725	339	}
92efd725	340	}
344100c4	341	} else {
	342	// Do Event Mixing
	343	if(fBGHandler==NULL){
	344	fBGHandler=new AliConversionAODBGHandlerRP(fEventCut->IsHeavyIon(),fMesonCut->UseTrackMultiplicity());
	345	}
92efd725	346
344100c4	347	for(Int_t nEventsInBG=0;nEventsInBG <fBGHandler->GetNBGEvents(fGoodGammas,fInputEvent);nEventsInBG++){
	348	AliGammaConversionPhotonVector *previousEventGammas = fBGHandler->GetBGGoodGammas(fGoodGammas,fInputEvent,nEventsInBG);
	349	if(previousEventGammas){
	350	// test weighted background
	351	Double_t weight=1.0;
	352	// Correct for the number of eventmixing:
	353	// N gammas -> (N-1) + (N-2) +(N-3) ...+ (N-(N-1)) using sum formula sum(i)=N(N-1)/2 -> N(N-1)/2
	354	// real combinations (since you cannot combine a photon with its own)
	355	// but BG leads to N_{a}*N_{b} combinations
	356	weight=0.5(Double_t(fGoodGammas->GetEntriesFast()-1))/Double_t(previousEventGammas->size());
	357	for(Int_t iCurrent=0;iCurrent<fGoodGammas->GetEntriesFast();iCurrent++){
	358	AliAODConversionPhoton gamma0 = (AliAODConversionPhoton)(fGoodGammas->At(iCurrent));
	359	for(UInt_t iPrevious=0;iPrevious<previousEventGammas->size();iPrevious++){
	360	AliAODConversionPhoton gamma1 = (AliAODConversionPhoton)(previousEventGammas->at(iPrevious));
	361	AliAODConversionMother BGcandidate(gamma0,gamma1);
	362	if(fMesonCut->MesonIsSelected(&BGcandidate,kFALSE)){
	363	if(MesonInMassWindow(&BGcandidate)){
	364	new((*fBGPi0s)[fBGPi0s->GetEntriesFast()]) AliAODConversionMother(BGcandidate);
	365	dynamic_cast<AliAODConversionMother*>(fBGPi0s->At(fBGPi0s->GetEntriesFast()-1))->SetWeight(weight);
	366	}
	367	}
	368	}
	369	}
92efd725	370	}
92efd725	371	}
92efd725	372	}
	373	}
	374	//________________________________________________________________________
	375	Double_t AliConversionSelection::GetMultiplicity(AliVEvent *inputEvent){
	376
344100c4	377	switch(fEventCut->GetMultiplicityMethod()){
	378	case 0:
	379	return Double_t(GetNumberOfPhotons());
	380	case 1:
	381	return Double_t(GetNumberOfChargedTracks(inputEvent));
	382	case 2:
	383	return GetVZEROMult(inputEvent);
	384	case 3:
	385	return GetSPDMult(inputEvent);
	386	case 9:
	387	return 1; // if mult is used as a weight, this number can be used to switch off weighting
	388	default:
	389	return 0;
92efd725	390	}
	391	}
	392
	393	//________________________________________________________________________
	394	Int_t AliConversionSelection::GetNumberOfChargedTracks(AliVEvent *inputEvent){
	395
344100c4	396	Int_t ntracks = 0;
92efd725	397
344100c4	398	AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(inputEvent);
	399	if(esdEvent) {
	400	if(!fESDTrackCuts){
	401	fESDTrackCuts= AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
	402	fESDTrackCuts->SetMaxDCAToVertexZ(2);
	403	Double_t etamax=fConversionCut->GetEtaCut();
	404	fESDTrackCuts->SetEtaRange(-etamax, etamax);
	405	fESDTrackCuts->SetPtRange(0.15);
	406	}
	407	for(Int_t iTracks = 0; iTracks < inputEvent->GetNumberOfTracks(); iTracks++){
	408	AliESDtrack* currentTrack = esdEvent->GetTrack(iTracks);
	409	if(!currentTrack) continue;
	410	if(fESDTrackCuts->AcceptTrack(currentTrack))ntracks++;
	411	}
	412	} else {
	413	for(Int_t ii=0; ii<inputEvent->GetNumberOfTracks(); ii++) {
	414	AliVTrack * track = dynamic_cast<AliVTrack*>(inputEvent->GetTrack(ii));
	415	if (track==NULL) continue;
	416	if(TMath::Abs(track->Eta())>fConversionCut->GetEtaCut())continue;
	417	if(track)ntracks++;
	418	}
	419	}
	420	return ntracks;
92efd725	421	}
	422
	423	//________________________________________________________________________
	424	Double_t AliConversionSelection::GetVZEROMult(AliVEvent *inputEvent){
	425
344100c4	426	AliVVZERO *vzero=inputEvent->GetVZEROData();
	427	Double_t multV0A=vzero->GetMTotV0A();
	428	Double_t multV0C=vzero->GetMTotV0C();
	429	Double_t mult=multV0A+multV0C;
92efd725	430
344100c4	431	return mult;
92efd725	432	}
	433
	434	//________________________________________________________________________
	435	Double_t AliConversionSelection::GetSPDMult(AliVEvent *inputEvent){
	436
344100c4	437	AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(inputEvent);
	438	if(esdEvent) {
	439	const AliMultiplicity *esdmult=esdEvent->GetMultiplicity();
	440	return esdmult->GetNumberOfITSClusters(1);
	441	} else {
	442	// AOD implementation
	443	AliAODHeader header=(AliAODHeader)inputEvent->GetHeader();
	444	return header->GetNumberOfITSClusters(1);
	445	}
	446	return 0;
92efd725	447	}
	448
	449	//________________________________________________________________________
	450	Int_t AliConversionSelection::GetEventNumber(AliVEvent *inputEvent){
	451
344100c4	452	AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(inputEvent);
	453	if(esdEvent) {
	454	return esdEvent->GetEventNumberInFile();
	455	} else{
	456	AliAODHeader header=(AliAODHeader)inputEvent->GetHeader();
	457	return header->GetEventNumberESDFile();
	458	}
	459	return 0;
92efd725	460	}
e5b6e8a6	461
	462	//________________________________________________________________________
	463	TString AliConversionSelection::GetCutString(){
344100c4	464	TString a= Form("%s_%s_%s",fEventCut->GetCutNumber().Data(), fConversionCut->GetCutNumber().Data(),fMesonCut->GetCutNumber().Data());
344100c4	465	return a;
e5b6e8a6	466	}
e5b6e8a6	467