1 #include "AliConversionSelection.h"
2 #include "AliAODHeader.h"
4 #include "AliMultiplicity.h"
8 // Author Daniel Lohner (Daniel.Lohner@cern.ch)
12 ClassImp(AliConversionSelection)
15 //________________________________________________________________________
16 AliConversionSelection::AliConversionSelection(AliConversionCuts *convCut, AliConversionMesonCuts *mesonCut) : TObject(),
19 fConversionCut(convCut),
27 fCurrentEventNumber(-1),
31 fInvMassRange[0]=0.05;
35 //________________________________________________________________________
36 AliConversionSelection::AliConversionSelection(TString convCut, TString mesonCut) : TObject(),
47 fCurrentEventNumber(-1),
51 fInvMassRange[0]=0.05;
54 fConversionCut = new AliConversionCuts();
55 fConversionCut -> InitializeCutsFromCutString(convCut.Data());
56 fMesonCut = new AliConversionMesonCuts();
57 fMesonCut -> InitializeCutsFromCutString(mesonCut.Data());
62 //________________________________________________________________________
63 AliConversionSelection::AliConversionSelection(const AliConversionSelection &ref) : TObject(ref),
74 fCurrentEventNumber(-1),
79 fConversionCut=new AliConversionCuts(*ref.fConversionCut);
80 fMesonCut=new AliConversionMesonCuts(*ref.fMesonCut);
82 fInvMassRange[0]=ref.fInvMassRange[0];
83 fInvMassRange[1]=ref.fInvMassRange[1];
86 //________________________________________________________________________
87 AliConversionSelection::~AliConversionSelection(){
98 delete fPi0Candidates;
106 delete fESDTrackCuts;
111 delete fConversionCut;
121 //________________________________________________________________________
122 Bool_t AliConversionSelection::ProcessEvent(TClonesArray *photons,AliVEvent *inputEvent,AliMCEvent *mcEvent){
123 fInputEvent=inputEvent;
127 Int_t eventnumber=GetEventNumber(inputEvent);
128 if(eventnumber==fCurrentEventNumber){
129 AliWarning("Event already analyzed! Return.");
133 fCurrentEventNumber=eventnumber;
136 // Initialize and Reset Arrays
137 if(fGoodGammas == NULL){
138 fGoodGammas=new TObjArray(30);
140 fGoodGammas->Clear();
142 if(fPi0Candidates == NULL){
143 fPi0Candidates = new TClonesArray("AliAODConversionMother",100);
145 fPi0Candidates->Delete();
148 fBGPi0s = new TClonesArray("AliAODConversionMother",100);
153 if(!photons||!fInputEvent)return kFALSE;
155 if(!fConversionCut->EventIsSelected(fInputEvent,fMCEvent))return kFALSE;
158 for(Int_t i = 0; i < photons->GetEntriesFast(); i++) {
159 AliAODConversionPhoton* gamma =dynamic_cast<AliAODConversionPhoton*>(photons->At(i));
161 if(!fConversionCut->PhotonIsSelected(gamma,fInputEvent))continue;
162 fGoodGammas->Add(gamma);
166 Double_t *fUnsmearedPx=NULL;
167 Double_t *fUnsmearedPy=NULL;
168 Double_t *fUnsmearedPz=NULL;
169 Double_t *fUnsmearedE=NULL;
171 if(fMesonCut->UseMCPSmearing() && fMCEvent){
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()];
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();
182 fMesonCut->SmearParticle(dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(gamma)));
186 // Reconstruct Pi0 and BG
187 CalculatePi0Candidates();
188 CalculateBackground();
189 if(fBGHandler)fBGHandler->AddEvent(fGoodGammas,fInputEvent);
192 if(fMesonCut->UseMCPSmearing() && fMCEvent){
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]);
199 delete[] fUnsmearedPx; fUnsmearedPx = 0x0;
200 delete[] fUnsmearedPy; fUnsmearedPy = 0x0;
201 delete[] fUnsmearedPz; fUnsmearedPz = 0x0;
202 delete[] fUnsmearedE; fUnsmearedE = 0x0;
208 //________________________________________________________________________
209 AliAODConversionMother* AliConversionSelection::GetPi0(Int_t index){
211 if(index>=0&&index<GetNumberOfPi0s()){
212 return dynamic_cast<AliAODConversionMother*>(fPi0Candidates->At(index));
217 //________________________________________________________________________
218 AliAODConversionMother* AliConversionSelection::GetBG(Int_t index){
220 if(index>=0&&index<GetNumberOfBGs()){
221 return dynamic_cast<AliAODConversionMother*>(fBGPi0s->At(index));
226 //________________________________________________________________________
227 AliAODConversionPhoton* AliConversionSelection::GetPhoton(Int_t index){
229 if(index>=0&&index<GetNumberOfPhotons()){
230 return dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(index));
235 //________________________________________________________________________
236 void AliConversionSelection::CalculatePi0Candidates(){
239 if(fGoodGammas->GetEntriesFast()>1){
241 for(Int_t firstGammaIndex=0;firstGammaIndex<fGoodGammas->GetEntriesFast()-1;firstGammaIndex++){
243 AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(firstGammaIndex));
244 if (gamma0==NULL) continue;
247 for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fGoodGammas->GetEntriesFast();secondGammaIndex++){
249 AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(secondGammaIndex));
250 if (gamma1==NULL) continue;
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;
255 AliAODConversionMother pi0cand(gamma0,gamma1);
256 pi0cand.SetLabels(firstGammaIndex,secondGammaIndex);
262 TParticle *mcgam0=gamma0->GetMCParticle(fMCEvent->Stack());
263 TParticle *mcgam1=gamma1->GetMCParticle(fMCEvent->Stack());
268 if(mcgam0->GetMother(0)==mcgam1->GetMother(0)){
270 pi0cand.SetMCLabel(mcgam0->GetMother(0));
275 if((fMesonCut->MesonIsSelected(&pi0cand))){
276 if(MesonInMassWindow(&pi0cand)){
279 new((*fPi0Candidates)[fPi0Candidates->GetEntriesFast()]) AliAODConversionMother(pi0cand);
287 //________________________________________________________________________
288 Bool_t AliConversionSelection::MesonInMassWindow(AliAODConversionMother *pi0cand)
290 if (pi0cand->M() > fInvMassRange[0] && pi0cand->M() < fInvMassRange[1] ){
296 //________________________________________________________________________
297 void AliConversionSelection::RotateParticle(AliAODConversionPhoton *gamma,Int_t nDegreesPMBackground){
300 fRandomizer=new TRandom3();
301 fRandomizer->SetSeed(0);
303 Double_t nRadiansPM = nDegreesPMBackground*TMath::Pi()/180;
304 Double_t rotationValue = fRandomizer->Rndm()*2*nRadiansPM + TMath::Pi()-nRadiansPM;
305 gamma->RotateZ(rotationValue);
308 //________________________________________________________________________
310 void AliConversionSelection::CalculateBackground(){
313 if(fMesonCut->UseRotationMethod()){
315 // Correct for the number of rotations
316 // BG is for rotation the same, except for factor NRotations
317 Double_t weight=1./Double_t(fMesonCut->GetNumberOfBGEvents());
319 for(Int_t firstGammaIndex=0;firstGammaIndex<fGoodGammas->GetEntriesFast();firstGammaIndex++){
321 AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(firstGammaIndex));
322 if (gamma0 ==NULL) continue;
323 for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fGoodGammas->GetEntriesFast();secondGammaIndex++){
324 AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(secondGammaIndex));
325 if (gamma1==NULL) continue;
326 if(!fConversionCut->PhotonIsSelected(gamma1,fInputEvent))continue;
327 for(Int_t nRandom=0;nRandom<fMesonCut->GetNumberOfBGEvents();nRandom++){
329 RotateParticle(gamma1,fMesonCut->NDegreesRotation());
331 AliAODConversionMother BGcandidate(gamma0,gamma1);
333 if(fMesonCut->MesonIsSelected(&BGcandidate,kFALSE)){
334 if(MesonInMassWindow(&BGcandidate)){
336 new((*fBGPi0s)[fBGPi0s->GetEntriesFast()]) AliAODConversionMother(BGcandidate);
338 dynamic_cast<AliAODConversionMother*>(fBGPi0s->At(fBGPi0s->GetEntriesFast()-1))->SetWeight(weight);
349 if(fBGHandler==NULL){
350 fBGHandler=new AliConversionAODBGHandlerRP(fConversionCut->IsHeavyIon(),fMesonCut->UseTrackMultiplicity());
353 for(Int_t nEventsInBG=0;nEventsInBG <fBGHandler->GetNBGEvents(fGoodGammas,fInputEvent);nEventsInBG++){
355 AliGammaConversionPhotonVector *previousEventGammas = fBGHandler->GetBGGoodGammas(fGoodGammas,fInputEvent,nEventsInBG);
357 if(previousEventGammas){
359 // test weighted background
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());
367 for(Int_t iCurrent=0;iCurrent<fGoodGammas->GetEntriesFast();iCurrent++){
369 AliAODConversionPhoton *gamma0 = (AliAODConversionPhoton*)(fGoodGammas->At(iCurrent));
371 for(UInt_t iPrevious=0;iPrevious<previousEventGammas->size();iPrevious++){
373 AliAODConversionPhoton *gamma1 = (AliAODConversionPhoton*)(previousEventGammas->at(iPrevious));
375 AliAODConversionMother BGcandidate(gamma0,gamma1);
377 if(fMesonCut->MesonIsSelected(&BGcandidate,kFALSE)){
378 if(MesonInMassWindow(&BGcandidate)){
380 new((*fBGPi0s)[fBGPi0s->GetEntriesFast()]) AliAODConversionMother(BGcandidate);
381 dynamic_cast<AliAODConversionMother*>(fBGPi0s->At(fBGPi0s->GetEntriesFast()-1))->SetWeight(weight);
390 //________________________________________________________________________
391 Double_t AliConversionSelection::GetMultiplicity(AliVEvent *inputEvent){
393 switch(fConversionCut->GetMultiplicityMethod())
396 return Double_t(GetNumberOfPhotons());
398 return Double_t(GetNumberOfChargedTracks(inputEvent));
400 return GetVZEROMult(inputEvent);
402 return GetSPDMult(inputEvent);
404 return 1; // if mult is used as a weight, this number can be used to switch off weighting
410 //________________________________________________________________________
411 Int_t AliConversionSelection::GetNumberOfChargedTracks(AliVEvent *inputEvent){
415 AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(inputEvent);
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);
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++;
430 for(Int_t ii=0; ii<inputEvent->GetNumberOfTracks(); ii++) {
431 AliVTrack * track = dynamic_cast<AliVTrack*>(inputEvent->GetTrack(ii));
432 if (track==NULL) continue;
433 if(TMath::Abs(track->Eta())>fConversionCut->GetEtaCut())continue;
441 //________________________________________________________________________
442 Double_t AliConversionSelection::GetVZEROMult(AliVEvent *inputEvent){
444 AliVVZERO *vzero=inputEvent->GetVZEROData();
445 Double_t multV0A=vzero->GetMTotV0A();
446 Double_t multV0C=vzero->GetMTotV0C();
447 Double_t mult=multV0A+multV0C;
452 //________________________________________________________________________
453 Double_t AliConversionSelection::GetSPDMult(AliVEvent *inputEvent){
455 AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(inputEvent);
457 const AliMultiplicity *esdmult=esdEvent->GetMultiplicity();
458 return esdmult->GetNumberOfITSClusters(1);
460 // AOD implementation
461 AliAODHeader *header=(AliAODHeader*)inputEvent->GetHeader();
462 return header->GetNumberOfITSClusters(1);
467 //________________________________________________________________________
468 Int_t AliConversionSelection::GetEventNumber(AliVEvent *inputEvent){
470 AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(inputEvent);
472 return esdEvent->GetEventNumberInFile();
475 AliAODHeader *header=(AliAODHeader*)inputEvent->GetHeader();
476 return header->GetEventNumberESDFile();
481 //________________________________________________________________________
482 TString AliConversionSelection::GetCutString(){
483 TString a= Form("%s%s",fConversionCut->GetCutNumber().Data(),fMesonCut->GetCutNumber().Data());