/************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * Authors: Svein Lindal, Daniel Lohner * * Version 1.0 * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ //////////////////////////////////////////////// //--------------------------------------------- // Class handling all kinds of selection cuts for // Gamma Conversion analysis //--------------------------------------------- //////////////////////////////////////////////// #include "AliConversionMesonCuts.h" #include "AliKFVertex.h" #include "AliAODTrack.h" #include "AliESDtrack.h" #include "AliAnalysisManager.h" #include "AliInputEventHandler.h" #include "AliMCEventHandler.h" #include "AliAODHandler.h" #include "AliPIDResponse.h" #include "TH1.h" #include "TH2.h" #include "AliStack.h" #include "AliAODConversionMother.h" #include "TObjString.h" #include "AliAODEvent.h" #include "AliESDEvent.h" #include "AliCentrality.h" #include "TList.h" class iostream; using namespace std; ClassImp(AliConversionMesonCuts) const char* AliConversionMesonCuts::fgkCutNames[AliConversionMesonCuts::kNCuts] = { "MesonKind", //0 "BackgroundScheme", //1 "NumberOfBGEvents", //2 "DegreesForRotationMethod", //3 "RapidityMesonCut", //4 "RCut", //5 "AlphaMesonCut", //6 "Chi2MesonCut", //7 "SharedElectronCuts", //8 "RejectToCloseV0s", //9 "UseMCPSmearing", //10 }; //________________________________________________________________________ AliConversionMesonCuts::AliConversionMesonCuts(const char *name,const char *title) : AliAnalysisCuts(name,title), fHistograms(NULL), fMesonKind(0), fMaxR(200), fChi2CutMeson(1000), fAlphaMinCutMeson(0), fAlphaCutMeson(1), fRapidityCutMeson(1), fUseRotationMethodInBG(kFALSE), fDoBG(kTRUE), fdoBGProbability(kFALSE), fUseTrackMultiplicityForBG(kFALSE), fnDegreeRotationPMForBG(0), fNumberOfBGEvents(0), fOpeningAngle(0.005), fDoToCloseV0sCut(kFALSE), fminV0Dist(200.), fDoSharedElecCut(kFALSE), fUseMCPSmearing(kFALSE), fPBremSmearing(0), fPSigSmearing(0), fPSigSmearingCte(0), fBrem(NULL), fRandom(0), fElectronLabelArraySize(500), fElectronLabelArray(NULL), fBackgroundHandler(0), fCutString(NULL), hMesonCuts(NULL), hMesonBGCuts(NULL) { for(Int_t jj=0;jjSetParameter(0,fPBremSmearing); fBrem->SetNpx(100000); } } //________________________________________________________________________ AliConversionMesonCuts::AliConversionMesonCuts(const AliConversionMesonCuts &ref) : AliAnalysisCuts(ref), fHistograms(NULL), fMesonKind(ref.fMesonKind), fMaxR(ref.fMaxR), fChi2CutMeson(ref.fChi2CutMeson), fAlphaMinCutMeson(ref.fAlphaMinCutMeson), fAlphaCutMeson(ref.fAlphaCutMeson), fRapidityCutMeson(ref.fRapidityCutMeson), fUseRotationMethodInBG(ref.fUseRotationMethodInBG), fDoBG(ref.fDoBG), fdoBGProbability(ref.fdoBGProbability), fUseTrackMultiplicityForBG(ref.fUseTrackMultiplicityForBG), fnDegreeRotationPMForBG(ref.fnDegreeRotationPMForBG), fNumberOfBGEvents(ref. fNumberOfBGEvents), fOpeningAngle(ref.fOpeningAngle), fDoToCloseV0sCut(ref.fDoToCloseV0sCut), fminV0Dist(ref.fminV0Dist), fDoSharedElecCut(ref.fDoSharedElecCut), fUseMCPSmearing(ref.fUseMCPSmearing), fPBremSmearing(ref.fPBremSmearing), fPSigSmearing(ref.fPSigSmearing), fPSigSmearingCte(ref.fPSigSmearingCte), fBrem(NULL), fRandom(ref.fRandom), fElectronLabelArraySize(ref.fElectronLabelArraySize), fElectronLabelArray(NULL), fBackgroundHandler(ref.fBackgroundHandler), fCutString(NULL), hMesonCuts(NULL), hMesonBGCuts(NULL) { // Copy Constructor for(Int_t jj=0;jjClone("fBrem"); // Histograms are not copied, if you need them, call InitCutHistograms } //________________________________________________________________________ AliConversionMesonCuts::~AliConversionMesonCuts() { // Destructor //Deleting fHistograms leads to seg fault it it's added to output collection of a task // if(fHistograms) // delete fHistograms; // fHistograms = NULL; if(fCutString != NULL){ delete fCutString; fCutString = NULL; } if(fElectronLabelArray){ delete fElectronLabelArray; fElectronLabelArray = NULL; } } //________________________________________________________________________ void AliConversionMesonCuts::InitCutHistograms(TString name){ // Initialize Cut Histograms for QA (only initialized and filled if function is called) TH1::AddDirectory(kFALSE); if(fHistograms != NULL){ delete fHistograms; fHistograms=NULL; } if(fHistograms==NULL){ fHistograms=new TList(); fHistograms->SetOwner(kTRUE); if(name=="")fHistograms->SetName(Form("ConvMesonCuts_%s",GetCutNumber().Data())); else fHistograms->SetName(Form("%s_%s",name.Data(),GetCutNumber().Data())); } // Meson Cuts hMesonCuts=new TH1F(Form("MesonCuts %s",GetCutNumber().Data()),"MesonCuts",10,-0.5,9.5); hMesonCuts->GetXaxis()->SetBinLabel(1,"in"); hMesonCuts->GetXaxis()->SetBinLabel(2,"undef rapidity"); hMesonCuts->GetXaxis()->SetBinLabel(3,"rapidity cut"); hMesonCuts->GetXaxis()->SetBinLabel(4,"opening angle"); hMesonCuts->GetXaxis()->SetBinLabel(5,"alpha max"); hMesonCuts->GetXaxis()->SetBinLabel(6,"alpha min"); hMesonCuts->GetXaxis()->SetBinLabel(7,"out"); fHistograms->Add(hMesonCuts); hMesonBGCuts=new TH1F(Form("MesonBGCuts %s",GetCutNumber().Data()),"MesonBGCuts",10,-0.5,9.5); hMesonBGCuts->GetXaxis()->SetBinLabel(1,"in"); hMesonBGCuts->GetXaxis()->SetBinLabel(2,"undef rapidity"); hMesonBGCuts->GetXaxis()->SetBinLabel(3,"rapidity cut"); hMesonBGCuts->GetXaxis()->SetBinLabel(4,"opening angle"); hMesonBGCuts->GetXaxis()->SetBinLabel(5,"alpha max"); hMesonBGCuts->GetXaxis()->SetBinLabel(6,"alpha min"); hMesonBGCuts->GetXaxis()->SetBinLabel(7,"out"); fHistograms->Add(hMesonBGCuts); TH1::AddDirectory(kTRUE); } //________________________________________________________________________ Bool_t AliConversionMesonCuts::MesonIsSelectedMC(TParticle *fMCMother,AliStack *fMCStack, Double_t fRapidityShift){ // Returns true for all pions within acceptance cuts for decay into 2 photons // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts if(!fMCStack)return kFALSE; if(fMCMother->GetPdgCode()==111 || fMCMother->GetPdgCode()==221){ if(fMCMother->R()>fMaxR) return kFALSE; // cuts on distance from collision point Double_t rapidity = 10.; if(fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0){ rapidity=8.-fRapidityShift; } else{ rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift; } // Rapidity Cut if(abs(rapidity)>fRapidityCutMeson)return kFALSE; // Select only -> 2y decay channel if(fMCMother->GetNDaughters()!=2)return kFALSE; for(Int_t i=0;i<2;i++){ TParticle *MDaughter=fMCStack->Particle(fMCMother->GetDaughter(i)); // Is Daughter a Photon? if(MDaughter->GetPdgCode()!=22)return kFALSE; // Is Photon in Acceptance? // if(bMCDaughtersInAcceptance){ // if(!PhotonIsSelectedMC(MDaughter,fMCStack)){return kFALSE;} // } } return kTRUE; } return kFALSE; } //________________________________________________________________________ Bool_t AliConversionMesonCuts::MesonIsSelectedMCDalitz(TParticle *fMCMother,AliStack *fMCStack, Double_t fRapidityShift){ // Returns true for all pions within acceptance cuts for decay into 2 photons // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts if(!fMCStack)return kFALSE; if(fMCMother->GetPdgCode()==111 || fMCMother->GetPdgCode()==221){ if(fMCMother->R()>fMaxR) return kFALSE; // cuts on distance from collision point Double_t rapidity = 10.; if(fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0){ rapidity=8.-fRapidityShift; } else{ rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift; } // Rapidity Cut if(abs(rapidity)>fRapidityCutMeson)return kFALSE; // Select only -> Dalitz decay channel if(fMCMother->GetNDaughters()!=3)return kFALSE; Int_t daughterPDGs[3] = {0,0,0}; Int_t index = 0; // iParticle->GetFirstDaughter(); idxPi0 <= iParticle->GetLastDaughter() for(Int_t i=fMCMother->GetFirstDaughter(); i<= fMCMother->GetLastDaughter();i++){ TParticle *MDaughter=fMCStack->Particle(i); // Is Daughter a Photon or an electron? daughterPDGs[index]=MDaughter->GetPdgCode(); index++; } for (Int_t j=0;j<2;j++){ for (Int_t i=0;i<2;i++){ if (daughterPDGs[i] > daughterPDGs[i+1]){ Int_t interMed = daughterPDGs[i] ; daughterPDGs[i] = daughterPDGs[i+1]; daughterPDGs[i+1] = interMed; } } } if (daughterPDGs[0] == -11 && daughterPDGs[1] == 11 && daughterPDGs[2] == 22) return kTRUE; } return kFALSE; } //________________________________________________________________________ Bool_t AliConversionMesonCuts::MesonIsSelectedMCChiC(TParticle *fMCMother,AliStack *fMCStack,Int_t & labelelectronChiC, Int_t & labelpositronChiC, Int_t & labelgammaChiC, Double_t fRapidityShift){ // Returns true for all ChiC within acceptance cuts for decay into JPsi + gamma -> e+ + e- + gamma // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts if(!fMCStack)return kFALSE; // if(fMCMother->GetPdgCode()==20443 ){ // return kFALSE; // } if(fMCMother->GetPdgCode()==10441 || fMCMother->GetPdgCode()==10443 || fMCMother->GetPdgCode()==445 ){ if(fMCMother->R()>fMaxR) return kFALSE; // cuts on distance from collision point Double_t rapidity = 10.; if(fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0){ rapidity=8.-fRapidityShift; } else{ rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift; } // Rapidity Cut if(abs(rapidity)>fRapidityCutMeson)return kFALSE; // Select only -> ChiC radiative (JPsi+gamma) decay channel if(fMCMother->GetNDaughters()!=2)return kFALSE; TParticle *jpsi = 0x0; TParticle *gamma = 0x0; TParticle *positron = 0x0; TParticle *electron = 0x0; Int_t labeljpsiChiC = -1; for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){ TParticle* temp = (TParticle*)fMCStack->Particle( index ); switch( temp->GetPdgCode() ) { case 443: jpsi = temp; labeljpsiChiC = index; break; case 22: gamma = temp; labelgammaChiC = index; break; } } if ( !jpsi || ! gamma) return kFALSE; if(jpsi->GetNDaughters()!=2)return kFALSE; for(Int_t index= jpsi->GetFirstDaughter();index<= jpsi->GetLastDaughter();index++){ TParticle* temp = (TParticle*)fMCStack->Particle( index ); switch( temp->GetPdgCode() ) { case -11: electron = temp; labelelectronChiC = index; break; case 11: positron = temp; labelpositronChiC = index; break; } } if( !electron || !positron) return kFALSE; if( positron && electron && gamma) return kTRUE; } return kFALSE; } ///________________________________________________________________________ Bool_t AliConversionMesonCuts::MesonIsSelected(AliAODConversionMother *pi0,Bool_t IsSignal, Double_t fRapidityShift) { // Selection of reconstructed Meson candidates // Use flag IsSignal in order to fill Fill different // histograms for Signal and Background TH1 *hist=0x0; if(IsSignal){hist=hMesonCuts;} else{hist=hMesonBGCuts;} Int_t cutIndex=0; if(hist)hist->Fill(cutIndex); cutIndex++; // Undefined Rapidity -> Floating Point exception if((pi0->E()+pi0->Pz())/(pi0->E()-pi0->Pz())<=0){ if(hist)hist->Fill(cutIndex); cutIndex++; return kFALSE; } else{ // PseudoRapidity Cut --> But we cut on Rapidity !!! cutIndex++; if(abs(pi0->Rapidity()-fRapidityShift)>fRapidityCutMeson){ if(hist)hist->Fill(cutIndex); return kFALSE; } } cutIndex++; // Opening Angle Cut //fOpeningAngle=2*TMath::ATan(0.134/pi0->P());// physical minimum opening angle if( pi0->GetOpeningAngle() < fOpeningAngle){ if(hist)hist->Fill(cutIndex); return kFALSE; } cutIndex++; // Alpha Max Cut if(pi0->GetAlpha()>fAlphaCutMeson){ if(hist)hist->Fill(cutIndex); return kFALSE; } cutIndex++; // Alpha Min Cut if(pi0->GetAlpha()Fill(cutIndex); return kFALSE; } cutIndex++; if(hist)hist->Fill(cutIndex); return kTRUE; } ///________________________________________________________________________ ///________________________________________________________________________ Bool_t AliConversionMesonCuts::UpdateCutString() { ///Update the cut string (if it has been created yet) if(fCutString && fCutString->GetString().Length() == kNCuts) { fCutString->SetString(GetCutNumber()); } else { return kFALSE; } return kTRUE; } ///________________________________________________________________________ Bool_t AliConversionMesonCuts::InitializeCutsFromCutString(const TString analysisCutSelection ) { // Initialize Cuts from a given Cut string AliInfo(Form("Set Meson Cutnumber: %s",analysisCutSelection.Data())); if(analysisCutSelection.Length()!=kNCuts) { AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts)); return kFALSE; } if(!analysisCutSelection.IsDigit()){ AliError("Cut selection contains characters"); return kFALSE; } const char *cutSelection = analysisCutSelection.Data(); #define ASSIGNARRAY(i) fCuts[i] = cutSelection[i] - '0' for(Int_t ii=0;iiGetTrackLabelPositive(); Int_t negLabel = photon->GetTrackLabelNegative(); fElectronLabelArray[nV0*2] = posLabel; fElectronLabelArray[(nV0*2)+1] = negLabel; } ///________________________________________________________________________ Bool_t AliConversionMesonCuts::RejectSharedElectronV0s(AliAODConversionPhoton* photon, Int_t nV0, Int_t nV0s){ Int_t posLabel = photon->GetTrackLabelPositive(); Int_t negLabel = photon->GetTrackLabelNegative(); for(Int_t i = 0; iGetConversionX(); Double_t posY = photon->GetConversionY(); Double_t posZ = photon->GetConversionZ(); for(Int_t i = 0;iGetEntries();i++){ if(nV0 == i) continue; AliAODConversionPhoton *photonComp = (AliAODConversionPhoton*) photons->At(i); Double_t posCompX = photonComp->GetConversionX(); Double_t posCompY = photonComp->GetConversionY(); Double_t posCompZ = photonComp->GetConversionZ(); Double_t dist = pow((posX - posCompX),2)+pow((posY - posCompY),2)+pow((posZ - posCompZ),2); if(dist < fminV0Dist*fminV0Dist){ if(photon->GetChi2perNDF() < photonComp->GetChi2perNDF()) return kTRUE; else { return kFALSE;} } } return kTRUE; } ///________________________________________________________________________ void AliConversionMesonCuts::SmearParticle(AliAODConversionPhoton* photon) { if (photon==NULL) return; Double_t facPBrem = 1.; Double_t facPSig = 0.; Double_t phi=0.; Double_t theta=0.; Double_t P=0.; P=photon->P(); phi=photon->Phi(); if( photon->P()!=0){ theta=acos( photon->Pz()/ photon->P()); } if( fPSigSmearing != 0. || fPSigSmearingCte!=0. ){ facPSig = TMath::Sqrt(fPSigSmearingCte*fPSigSmearingCte+fPSigSmearing*fPSigSmearing*P*P)*fRandom.Gaus(0.,1.); } if( fPBremSmearing != 1.){ if(fBrem!=NULL){ facPBrem = fBrem->GetRandom(); } } photon->SetPx(facPBrem* (1+facPSig)* P*sin(theta)*cos(phi)) ; photon->SetPy(facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)) ; photon->SetPz(facPBrem* (1+facPSig)* P*cos(theta)) ; photon->SetE(photon->P()); }