[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),
    fCurrentEventNumber(-1),
    fIsOwner(kFALSE)
{
    // Default Values
    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),
    fCurrentEventNumber(-1),
    fIsOwner(kTRUE)
{
    // Default Values
    fInvMassRange[0]=0.05;
    fInvMassRange[1]=0.3;

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

}


//________________________________________________________________________
AliConversionSelection::AliConversionSelection(const AliConversionSelection &ref) : TObject(ref),
    fInputEvent(NULL),
    fMCEvent(NULL),
    fConversionCut(NULL),
    fMesonCut(NULL),
    fESDTrackCuts(NULL),
    fGoodGammas(NULL),
    fPi0Candidates(NULL),
    fBGPi0s(NULL),
    fRandomizer(NULL),
    fBGHandler(NULL),
    fCurrentEventNumber(-1),
    fIsOwner(kTRUE)
{
    // Copy Constructor

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