[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(AliConversionCuts *convCut, AliConversionMesonCuts *mesonCut) : TObject(),
    fInputEvent(NULL),
    fMCEvent(NULL),
    fConversionCut(convCut),
    fMesonCut(mesonCut),
    fESDTrackCuts(NULL),
    fGoodGammas(NULL),
    fPi0Candidates(NULL),
    fBGPi0s(NULL),
    fRandomizer(NULL),
    fBGHandler(NULL),
    fInvMassRange(NULL),
    fUnsmearedPx(NULL),
    fUnsmearedPy(NULL),
    fUnsmearedPz(NULL),
    fUnsmearedE(NULL),
    fCurrentEventNumber(-1),
    fIsOwner(kFALSE)
{
    // Default Values
    fInvMassRange=new Double_t[2];
    fInvMassRange[0]=0.05;
    fInvMassRange[1]=0.3;
}

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

    fConversionCut = new AliConversionCuts();
    fConversionCut -> InitializeCutsFromCutString(convCut.Data());
    fMesonCut = new AliConversionMesonCuts();
    fMesonCut -> InitializeCutsFromCutString(mesonCut.Data());

}

//________________________________________________________________________
AliConversionSelection::~AliConversionSelection(){

    if(fBGHandler){
	delete fBGHandler;
	fBGHandler=NULL;
    }
    if(fInvMassRange){
	delete fInvMassRange;
	fInvMassRange=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(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(!fConversionCut->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
    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));

	    // Combine Photons

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

		AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(secondGammaIndex));

		//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));

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