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

//________________________________________________________________________
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(fConversionCut->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();
	    fConversionCut->SmearParticle(dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(gamma)));
	}
    }

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

    // Undo MC Smearing
    if(fConversionCut->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((fConversionCut->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(fConversionCut->UseRotationMethod()){

	// Correct for the number of rotations
	// BG is for rotation the same, except for factor NRotations
	Double_t weight=1./Double_t(fConversionCut->NumberOfRotationEvents());

	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<fConversionCut->NumberOfRotationEvents();nRandom++){

		    RotateParticle(gamma1,fConversionCut->NDegreesRotation());

		    AliAODConversionMother BGcandidate(gamma0,gamma1);

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