#include "AliESDEvent.h"
#include "AliCentrality.h"
#include "TList.h"
+#include "TFile.h"
#include "AliLog.h"
#include "AliGenCocktailEventHeader.h"
#include "AliGenDPMjetEventHeader.h"
#include "AliGenHijingEventHeader.h"
#include "AliTriggerAnalysis.h"
#include "AliV0ReaderV1.h"
+#include "AliAODMCParticle.h"
+#include "AliAODMCHeader.h"
class iostream;
fDoPhotonAsymmetryCut(kTRUE),
fMinPPhotonAsymmetryCut(100.),
fMinPhotonAsymmetry(0.),
- fIsHeavyIon(kFALSE),
+ fIsHeavyIon(0),
fDetectorCentrality(0),
fModCentralityClass(0),
fMaxVertexZ(10),
fNotRejectedEnd(NULL),
fGeneratorNames(NULL),
fCutString(NULL),
+ fUtils(NULL),
+ fEtaShift(0.0),
+ fDoEtaShift(kFALSE),
+ fDoReweightHistoMCPi0(kFALSE),
+ fDoReweightHistoMCEta(kFALSE),
+ fDoReweightHistoMCK0s(kFALSE),
+ fPathTrFReweighting(""),
+ fNameHistoReweightingPi0(""),
+ fNameHistoReweightingEta(""),
+ fNameHistoReweightingK0s(""),
hdEdxCuts(NULL),
hTPCdEdxbefore(NULL),
hTPCdEdxafter(NULL),
hCentrality(NULL),
hVertexZ(NULL),
hTriggerClass(NULL),
- hTriggerClassSelected(NULL)
+ hTriggerClassSelected(NULL),
+ hReweightMCHistPi0(NULL),
+ hReweightMCHistEta(NULL),
+ hReweightMCHistK0s(NULL),
+ fPreSelCut(kFALSE),
+ fTriggerSelectedManually(kFALSE),
+ fSpecialTriggerName("")
+
{
InitPIDResponse();
for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=0;}
fCutString=new TObjString((GetCutNumber()).Data());
fElectronLabelArray = new Int_t[fElectronArraySize];
+ fUtils = new AliAnalysisUtils();
+ //if you do not want to apply the cut on the distance between the SPD and TRK vertex:
+ //fUtils->SetCutOnZVertexSPD(kFALSE);
+
}
fNSigmaMass(ref.fNSigmaMass),
fUseEtaMinCut(ref.fUseEtaMinCut),
fUseOnFlyV0Finder(ref.fUseOnFlyV0Finder),
- fDoPhotonAsymmetryCut(fDoPhotonAsymmetryCut),
+ fDoPhotonAsymmetryCut(ref.fDoPhotonAsymmetryCut),
fMinPPhotonAsymmetryCut(ref.fMinPPhotonAsymmetryCut),
fMinPhotonAsymmetry(ref.fMinPhotonAsymmetry),
fIsHeavyIon(ref.fIsHeavyIon),
fNotRejectedEnd(NULL),
fGeneratorNames(ref.fGeneratorNames),
fCutString(NULL),
+ fUtils(NULL),
+ fEtaShift(ref.fEtaShift),
+ fDoEtaShift(ref.fDoEtaShift),
+ fDoReweightHistoMCPi0(ref.fDoReweightHistoMCPi0),
+ fDoReweightHistoMCEta(ref.fDoReweightHistoMCEta),
+ fDoReweightHistoMCK0s(ref.fDoReweightHistoMCK0s),
+ fPathTrFReweighting(ref.fPathTrFReweighting),
+ fNameHistoReweightingPi0(ref.fNameHistoReweightingPi0),
+ fNameHistoReweightingEta(ref.fNameHistoReweightingEta),
+ fNameHistoReweightingK0s(ref.fNameHistoReweightingK0s),
hdEdxCuts(NULL),
hTPCdEdxbefore(NULL),
hTPCdEdxafter(NULL),
hCentrality(NULL),
hVertexZ(NULL),
hTriggerClass(NULL),
- hTriggerClassSelected(NULL)
+ hTriggerClassSelected(NULL),
+ hReweightMCHistPi0(NULL),
+ hReweightMCHistEta(NULL),
+ hReweightMCHistK0s(NULL),
+ fPreSelCut(ref.fPreSelCut),
+ fTriggerSelectedManually(ref.fTriggerSelectedManually),
+ fSpecialTriggerName(ref.fSpecialTriggerName)
{
// Copy Constructor
for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=ref.fCuts[jj];}
fCutString=new TObjString((GetCutNumber()).Data());
fElectronLabelArray = new Int_t[fElectronArraySize];
+ fUtils = new AliAnalysisUtils();
// dont copy histograms (if you like histograms, call InitCutHistograms())
+
}
delete[] fGeneratorNames;
fGeneratorNames = NULL;
}
+ if(fUtils){
+ delete fUtils;
+ fUtils = NULL;
+ }
+
}
//________________________________________________________________________
void AliConversionCuts::InitCutHistograms(TString name, Bool_t preCut){
// Initialize Cut Histograms for QA (only initialized and filled if function is called)
+ TH1::AddDirectory(kFALSE);
if(fHistograms != NULL){
delete fHistograms;
}
if(fHistograms==NULL){
fHistograms=new TList();
+ fHistograms->SetOwner(kTRUE);
if(name=="")fHistograms->SetName(Form("ConvCuts_%s",GetCutNumber().Data()));
else fHistograms->SetName(Form("%s_%s",name.Data(),GetCutNumber().Data()));
}
+ if (hReweightMCHistPi0) fHistograms->Add(hReweightMCHistPi0);
+ if (hReweightMCHistEta) fHistograms->Add(hReweightMCHistEta);
+ if (hReweightMCHistK0s) fHistograms->Add(hReweightMCHistK0s);
+
// IsPhotonSelected
hCutIndex=new TH1F(Form("IsPhotonSelected %s",GetCutNumber().Data()),"IsPhotonSelected",10,-0.5,9.5);
hCutIndex->GetXaxis()->SetBinLabel(kPhotonIn+1,"in");
hVertexZ=new TH1F(Form("VertexZ %s",GetCutNumber().Data()),"VertexZ",1000,-50,50);
fHistograms->Add(hVertexZ);
- hTriggerClass= new TH1F(Form("OfflineTrigger %s",GetCutNumber().Data()),"OfflineTrigger",34,-0.5,33.5);
+ hTriggerClass= new TH1F(Form("OfflineTrigger %s",GetCutNumber().Data()),"OfflineTrigger",35,-0.5,34.5);
hTriggerClass->GetXaxis()->SetBinLabel( 1,"kMB");
hTriggerClass->GetXaxis()->SetBinLabel( 2,"kINT7");
hTriggerClass->GetXaxis()->SetBinLabel( 3,"kMUON");
hTriggerClass->GetXaxis()->SetBinLabel(26,"kMuonUnlikeLowPt8");
hTriggerClass->GetXaxis()->SetBinLabel(27,"kMuonUnlikeLowPt0");
hTriggerClass->GetXaxis()->SetBinLabel(28,"kUserDefined");
- hTriggerClass->GetXaxis()->SetBinLabel(29,"kFastOnly");
- hTriggerClass->GetXaxis()->SetBinLabel(30,"kAnyINT");
- hTriggerClass->GetXaxis()->SetBinLabel(31,"kAny");
- hTriggerClass->GetXaxis()->SetBinLabel(32,"V0AND");
- hTriggerClass->GetXaxis()->SetBinLabel(33,"NOT kFastOnly");
- hTriggerClass->GetXaxis()->SetBinLabel(34,"failed Physics Selection");
+ hTriggerClass->GetXaxis()->SetBinLabel(29,"kTRD");
+ hTriggerClass->GetXaxis()->SetBinLabel(30,"kFastOnly");
+ hTriggerClass->GetXaxis()->SetBinLabel(31,"kAnyINT");
+ hTriggerClass->GetXaxis()->SetBinLabel(32,"kAny");
+ hTriggerClass->GetXaxis()->SetBinLabel(33,"V0AND");
+ hTriggerClass->GetXaxis()->SetBinLabel(34,"NOT kFastOnly");
+ hTriggerClass->GetXaxis()->SetBinLabel(35,"failed Physics Selection");
fHistograms->Add(hTriggerClass);
-
- hTriggerClassSelected= new TH1F(Form("OfflineTriggerSelected %s",GetCutNumber().Data()),"OfflineTriggerSelected",33,-0.5,32.5);
+ }
+ if(!preCut){
+ hTriggerClassSelected= new TH1F(Form("OfflineTriggerSelected %s",GetCutNumber().Data()),"OfflineTriggerSelected",34,-0.5,33.5);
hTriggerClassSelected->GetXaxis()->SetBinLabel( 1,"kMB");
hTriggerClassSelected->GetXaxis()->SetBinLabel( 2,"kINT7");
hTriggerClassSelected->GetXaxis()->SetBinLabel( 3,"kMUON");
hTriggerClassSelected->GetXaxis()->SetBinLabel(26,"kMuonUnlikeLowPt8");
hTriggerClassSelected->GetXaxis()->SetBinLabel(27,"kMuonUnlikeLowPt0");
hTriggerClassSelected->GetXaxis()->SetBinLabel(28,"kUserDefined");
- hTriggerClassSelected->GetXaxis()->SetBinLabel(29,"kFastOnly");
- hTriggerClassSelected->GetXaxis()->SetBinLabel(30,"kAnyINT");
- hTriggerClassSelected->GetXaxis()->SetBinLabel(31,"kAny");
- hTriggerClassSelected->GetXaxis()->SetBinLabel(32,"V0AND");
- hTriggerClassSelected->GetXaxis()->SetBinLabel(33,"NOT kFastOnly");
+ hTriggerClassSelected->GetXaxis()->SetBinLabel(29,"kTRD");
+ hTriggerClassSelected->GetXaxis()->SetBinLabel(30,"kFastOnly");
+ hTriggerClassSelected->GetXaxis()->SetBinLabel(31,"kAnyINT");
+ hTriggerClassSelected->GetXaxis()->SetBinLabel(32,"kAny");
+ hTriggerClassSelected->GetXaxis()->SetBinLabel(33,"V0AND");
+ hTriggerClassSelected->GetXaxis()->SetBinLabel(34,"NOT kFastOnly");
fHistograms->Add(hTriggerClassSelected);
}
+ TH1::AddDirectory(kTRUE);
}
//________________________________________________________________________
cutindex++;
// Check for MC event
- if(fMCEvent){
+ if(fMCEvent && fInputEvent->IsA()==AliESDEvent::Class()){
// Check if MC event is correctly loaded
AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
if (!mcHandler){
}
// Event Trigger
- if(!IsTriggerSelected()){
+ if(!IsTriggerSelected(fInputEvent)){
if(hV0EventCuts)hV0EventCuts->Fill(cutindex);
fEventQuality = 3;
return kFALSE;
if(fInputEvent->IsA()==AliESDEvent::Class()){
AliTriggerAnalysis fTriggerAnalysis;// = new AliTriggerAnalysis;
fHasV0AND = fTriggerAnalysis.IsOfflineTriggerFired((AliESDEvent*)fInputEvent, AliTriggerAnalysis::kV0AND);
- if(fHasV0AND&&hTriggerClass)hTriggerClass->Fill(31);
+ if(fHasV0AND&&hTriggerClass)hTriggerClass->Fill(32);
}
+// cout << "event number " << ((AliESDEvent*)fInputEvent)->GetEventNumberInFile() << " entered"<< endl;
+
// Number of Contributors Cut
if(GetNumberOfContributorsVtx(fInputEvent)<=0) {
if (particle->GetPdgCode() == 22){
- if(particle->R() > fMaxR) return kFALSE;
- if(abs(particle->Eta())> fEtaCut || abs(particle->Eta())< fEtaCutMin) return kFALSE;
+
+ if( particle->Eta() > (fEtaCut + fEtaShift) || particle->Eta() < (-fEtaCut + fEtaShift) )
+ return kFALSE;
+ if(fEtaCutMin>-0.1){
+ if( particle->Eta() < (fEtaCutMin + fEtaShift) && particle->Eta() > (-fEtaCutMin + fEtaShift) )
+ return kFALSE;
+ }
if(particle->GetMother(0) >-1 && fMCStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){
return kFALSE; // no photon as mothers!
return kFALSE; // no reconstruction below the Pt cut
}
- if( abs(ePos->Eta())> fEtaCut || abs(ePos->Eta())< fEtaCutMin ||
- abs(eNeg->Eta())> fEtaCut || abs(eNeg->Eta())< fEtaCutMin ) {
+ if( ePos->Eta() > (fEtaCut + fEtaShift) || ePos->Eta() < (-fEtaCut + fEtaShift) ||
+ eNeg->Eta() > (fEtaCut + fEtaShift) || eNeg->Eta() < (-fEtaCut + fEtaShift) )
return kFALSE;
+
+ if(fEtaCutMin > -0.1){
+ if( (ePos->Eta() < (fEtaCutMin + fEtaShift) && ePos->Eta() > (-fEtaCutMin + fEtaShift)) ||
+ (eNeg->Eta() < (fEtaCutMin + fEtaShift) && eNeg->Eta() > (-fEtaCutMin + fEtaShift)) )
+ return kFALSE;
}
if(ePos->R()>fMaxR){
}
return kFALSE;
}
+///________________________________________________________________________
+Bool_t AliConversionCuts::PhotonIsSelectedAODMC(AliAODMCParticle *particle,TClonesArray *aodmcArray,Bool_t checkForConvertedGamma){
+ // MonteCarlo Photon Selection
+
+ if(!aodmcArray)return kFALSE;
+
+ if (particle->GetPdgCode() == 22){
+ if( particle->Eta() > (fEtaCut + fEtaShift) || particle->Eta() < (-fEtaCut + fEtaShift) )
+ return kFALSE;
+ if(fEtaCutMin>-0.1){
+ if( particle->Eta() < (fEtaCutMin + fEtaShift) && particle->Eta() > (-fEtaCutMin + fEtaShift) )
+ return kFALSE;
+ }
+
+ if(particle->GetMother() > -1){
+ if((static_cast<AliAODMCParticle*>(aodmcArray->At(particle->GetMother())))->GetPdgCode() == 22){
+ return kFALSE; // no photon as mothers!
+ }
+ if(!(static_cast<AliAODMCParticle*>(aodmcArray->At(particle->GetMother()))->IsPrimary())){
+ return kFALSE; // the gamma has a mother, and it is not a primary particle
+ }
+ }
+
+ if(!checkForConvertedGamma) return kTRUE; // return in case of accepted gamma
+
+ // looking for conversion gammas (electron + positron from pairbuilding (= 5) )
+ AliAODMCParticle* ePos = NULL;
+ AliAODMCParticle* eNeg = NULL;
+
+ if(particle->GetNDaughters() >= 2){
+ for(Int_t daughterIndex=particle->GetDaughter(0);daughterIndex<=particle->GetDaughter(1);daughterIndex++){
+ AliAODMCParticle *tmpDaughter = static_cast<AliAODMCParticle*>(aodmcArray->At(daughterIndex));
+ if(!tmpDaughter) continue;
+ if(((tmpDaughter->GetMCProcessCode())) == 5){ // STILL A BUG IN ALIROOT >>8 HAS TPO BE REMOVED AFTER FIX
+ if(tmpDaughter->GetPdgCode() == 11){
+ eNeg = tmpDaughter;
+ } else if(tmpDaughter->GetPdgCode() == -11){
+ ePos = tmpDaughter;
+ }
+ }
+ }
+ }
+
+ if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production
+ return kFALSE;
+ }
+
+ if(ePos->Pt()<fSinglePtCut || eNeg->Pt()<fSinglePtCut){
+ return kFALSE; // no reconstruction below the Pt cut
+ }
+
+ if( ePos->Eta() > (fEtaCut + fEtaShift) || ePos->Eta() < (-fEtaCut + fEtaShift) ||
+ eNeg->Eta() > (fEtaCut + fEtaShift) || eNeg->Eta() < (-fEtaCut + fEtaShift) )
+ return kFALSE;
+ if(fEtaCutMin > -0.1){
+ if( (ePos->Eta() < (fEtaCutMin + fEtaShift) && ePos->Eta() > (-fEtaCutMin + fEtaShift)) ||
+ (eNeg->Eta() < (fEtaCutMin + fEtaShift) && eNeg->Eta() > (-fEtaCutMin + fEtaShift)) )
+ return kFALSE;
+ }
+
+ Double_t rPos = sqrt( (ePos->Xv()*ePos->Xv()) + (ePos->Yv()*ePos->Yv()) );
+ Double_t rNeg = sqrt( (eNeg->Xv()*eNeg->Xv()) + (eNeg->Yv()*eNeg->Yv()) );
+
+ if(rPos>fMaxR){
+ return kFALSE; // cuts on distance from collision point
+ }
+ if(abs(ePos->Zv()) > fMaxZ){
+ return kFALSE; // outside material
+ }
+ if(abs(eNeg->Zv()) > fMaxZ){
+ return kFALSE; // outside material
+ }
+
+ if( rPos <= ((abs(ePos->Zv()) * fLineCutZRSlope) - fLineCutZValue)){
+ return kFALSE; // line cut to exclude regions where we do not reconstruct
+ } else if ( fEtaCutMin != -0.1 && rPos >= ((abs(ePos->Zv()) * fLineCutZRSlopeMin) - fLineCutZValueMin)){
+ return kFALSE;
+ }
+ if( rNeg <= ((abs(eNeg->Zv()) * fLineCutZRSlope) - fLineCutZValue)){
+ return kFALSE; // line cut to exclude regions where we do not reconstruct
+ } else if ( fEtaCutMin != -0.1 && rNeg >= ((abs(eNeg->Zv()) * fLineCutZRSlopeMin) - fLineCutZValueMin)){
+ return kFALSE;
+ }
+
+ return kTRUE;
+ //if(AcceptanceCut(particle,ePos,eNeg))return kTRUE;
+ }
+ return kFALSE;
+}
///________________________________________________________________________
Bool_t AliConversionCuts::PhotonCuts(AliConversionPhotonBase *photon,AliVEvent *event)
{ // Specific Photon Cuts
FillPhotonCutIndex(kPhotonIn);
+ if(event->IsA()==AliESDEvent::Class()) {
+ if(!SelectV0Finder( ( ((AliESDEvent*)event)->GetV0(photon->GetV0Index()))->GetOnFlyStatus() ) ){
+ FillPhotonCutIndex(kOnFly);
+ return kFALSE;
+ }
+ }
+ // else if(event->IsA()==AliAODEvent::Class()) {
+ // if(!SelectV0Finder( ( ((AliAODEvent*)event)->GetV0(photon->GetV0Index())) ) ){
+ // FillPhotonCutIndex(kOnFly);
+ // return kFALSE;
+ // }
+ // }
+
// Get Tracks
AliVTrack * negTrack = GetTrack(event, photon->GetTrackLabelNegative());
AliVTrack * posTrack = GetTrack(event, photon->GetTrackLabelPositive());
return kFALSE;
}
- // dEdx Cuts
- if(!dEdxCuts(negTrack) || !dEdxCuts(posTrack)) {
- FillPhotonCutIndex(kdEdxCuts);
- return kFALSE;
- }
-
// Track Cuts
if(!TracksAreSelected(negTrack, posTrack)){
FillPhotonCutIndex(kTrackCuts);
return kFALSE;
}
+ // dEdx Cuts
+ if(!dEdxCuts(negTrack) || !dEdxCuts(posTrack)) {
+ FillPhotonCutIndex(kdEdxCuts);
+ return kFALSE;
+ }
+
// Photon Cuts
if(!PhotonCuts(photon,event)){
FillPhotonCutIndex(kPhotonCuts);
cutIndex++;
- if(abs(photon->GetPhotonEta())> fEtaCut || abs(photon->GetPhotonEta())< fEtaCutMin){
+ if( photon->GetPhotonEta() > (fEtaCut + fEtaShift) || photon->GetPhotonEta() < (-fEtaCut + fEtaShift) ){
if(hAcceptanceCuts)hAcceptanceCuts->Fill(cutIndex);
return kFALSE;
}
+ if(fEtaCutMin>-0.1){
+ if( photon->GetPhotonEta() < (fEtaCutMin + fEtaShift) && photon->GetPhotonEta() > (-fEtaCutMin + fEtaShift) ){
+ if(hAcceptanceCuts)hAcceptanceCuts->Fill(cutIndex);
+ return kFALSE;
+ }
+ }
cutIndex++;
-
if(photon->GetPhotonPt()<fPtCut){
if(hAcceptanceCuts)hAcceptanceCuts->Fill(cutIndex);
return kFALSE;
cutIndex++;
// Acceptance
-
- if(abs(negTrack->Eta()) > fEtaCut || abs(negTrack->Eta()) < fEtaCutMin ||
- abs(posTrack->Eta())> fEtaCut || abs(posTrack->Eta())< fEtaCutMin) {
+ if( posTrack->Eta() > (fEtaCut + fEtaShift) || posTrack->Eta() < (-fEtaCut + fEtaShift) ||
+ negTrack->Eta() > (fEtaCut + fEtaShift) || negTrack->Eta() < (-fEtaCut + fEtaShift) ){
if(hTrackCuts)hTrackCuts->Fill(cutIndex);
return kFALSE;
}
+ if(fEtaCutMin>-0.1){
+ if( (posTrack->Eta() < (fEtaCutMin + fEtaShift) && posTrack->Eta() > (-fEtaCutMin + fEtaShift)) ||
+ (negTrack->Eta() < (fEtaCutMin + fEtaShift) && negTrack->Eta() > (-fEtaCutMin + fEtaShift)) ){
+ if(hTrackCuts)hTrackCuts->Fill(cutIndex);
+ return kFALSE;
+ }
+ }
cutIndex++;
// Single Pt Cut
return track;
} else {
- for(Int_t ii=0; ii<event->GetNumberOfTracks(); ii++) {
- AliVTrack * track = dynamic_cast<AliVTrack*>(event->GetTrack(ii));
-
- if(track) {
- if(track->GetID() == label) {
- return track;
+ AliVTrack * track = 0x0;
+ if(((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1"))->AreAODsRelabeled()){
+ track = dynamic_cast<AliVTrack*>(event->GetTrack(label));
+ return track;
+ }
+ else{
+ for(Int_t ii=0; ii<event->GetNumberOfTracks(); ii++) {
+ track = dynamic_cast<AliVTrack*>(event->GetTrack(ii));
+ if(track){
+ if(track->GetID() == label) {
+ return track;
+ }
}
}
}
}
-
//AliDebug(5,(Form("track not found %d %d",label,event->GetNumberOfTracks()));
return NULL;
}
return kFALSE;
}
- if(abs(particle->Eta())> fEtaCut || abs(particle->Eta())< fEtaCutMin){
+
+ if( particle->Eta() > (fEtaCut + fEtaShift) || particle->Eta() < (-fEtaCut + fEtaShift) ){
return kFALSE;
}
-
- if(abs(ePos->Eta())> fEtaCut || abs(ePos->Eta())< fEtaCutMin){
+ if( ePos->Eta() > (fEtaCut + fEtaShift) || ePos->Eta() < (-fEtaCut + fEtaShift) ){
return kFALSE;
}
-
- if(abs(eNeg->Eta())> fEtaCut || abs(eNeg->Eta())< fEtaCutMin){
+ if( eNeg->Eta() > (fEtaCut + fEtaShift) || eNeg->Eta() < (-fEtaCut + fEtaShift) ){
return kFALSE;
}
+ if(fEtaCutMin>-0.1){
+ if( particle->Eta() < (fEtaCutMin + fEtaShift) && particle->Eta() > (-fEtaCutMin + fEtaShift) ){
+ return kFALSE;
+ }
+ if( ePos->Eta() < (fEtaCutMin + fEtaShift) && ePos->Eta() > (-fEtaCutMin + fEtaShift) ){
+ return kFALSE;
+ }
+ if( eNeg->Eta() < (fEtaCutMin + fEtaShift) && eNeg->Eta() > (-fEtaCutMin + fEtaShift) ){
+ return kFALSE;
+ }
+ }
if( ePos->Pt()< fSinglePtCut || eNeg->Pt()< fSinglePtCut){
return kFALSE;
return kTRUE;
}
+void AliConversionCuts::LoadReweightingHistosMCFromFile() {
+
+ AliInfo("Entering loading of histograms for weighting");
+ TFile *f = TFile::Open(fPathTrFReweighting.Data());
+ if(!f){
+ AliError(Form("file for weighting %s not found",fPathTrFReweighting.Data()));
+ return;
+ }
+ if (fNameHistoReweightingPi0.CompareTo("") != 0 && fDoReweightHistoMCPi0 ){
+ hReweightMCHistPi0 = (TH1D*)f->Get(fNameHistoReweightingPi0.Data());
+ if (hReweightMCHistPi0) AliInfo(Form("%s has been loaded from %s", fNameHistoReweightingPi0.Data(),fPathTrFReweighting.Data() ));
+ else AliWarning(Form("%s not found in %s",fPathTrFReweighting.Data(), fNameHistoReweightingPi0.Data() ));
+ }
+ if (fNameHistoReweightingEta.CompareTo("") != 0 && fDoReweightHistoMCEta){
+ hReweightMCHistEta = (TH1D*)f->Get(fNameHistoReweightingEta.Data());
+ if (hReweightMCHistEta) AliInfo(Form("%s has been loaded from %s", fNameHistoReweightingEta.Data(),fPathTrFReweighting.Data() ));
+ else AliWarning(Form("%s not found in %s",fPathTrFReweighting.Data(), fNameHistoReweightingEta.Data() ));
+
+ }
+ if (fNameHistoReweightingK0s.CompareTo("") != 0 && fDoReweightHistoMCK0s){
+ hReweightMCHistK0s = (TH1D*)f->Get(fNameHistoReweightingK0s.Data());
+ if (hReweightMCHistK0s) AliInfo(Form("%s has been loaded from %s", fNameHistoReweightingK0s.Data(),fPathTrFReweighting.Data() ));
+ else AliWarning(Form("%s not found in %s",fPathTrFReweighting.Data(), fNameHistoReweightingK0s.Data() ));
+
+ }
+
+}
+
+
///________________________________________________________________________
Bool_t AliConversionCuts::InitializeCutsFromCutString(const TString analysisCutSelection ) {
// Initialize Cuts from a given Cut string
-
+ if(fDoReweightHistoMCPi0 || fDoReweightHistoMCEta || fDoReweightHistoMCK0s) LoadReweightingHistosMCFromFile();
+
AliInfo(Form("Set Photoncut Number: %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));
//PrintCuts();
- // Set StandardTriggers
-
- if(fIsHeavyIon)SelectCollisionCandidates(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral);
- else SelectCollisionCandidates(AliVEvent::kMB);
-
return kTRUE;
}
///________________________________________________________________________
fDetectorCentrality=0;
fModCentralityClass=5;
break;
+ case 8:
+ fIsHeavyIon=2;
+ fDetectorCentrality=0;
+ break;
+ case 9:
+ fIsHeavyIon=2;
+ fDetectorCentrality=1;
+ break;
default:
AliError(Form("SetHeavyIon not defined %d",isHeavyIon));
return kFALSE;
}
return kTRUE;
}
+
//___________________________________________________________________
Bool_t AliConversionCuts::SetCentralityMin(Int_t minCentrality)
{
case 3:
fSpecialTrigger=3; // V0AND plus with SDD requested
break;
-
+ // allows to run MB & 6 other different trigger classes in parallel with the same photon cut
+ case 4:
+ fSpecialTrigger=4; // different trigger class as MB
+ fTriggerSelectedManually = kTRUE;
+ break;
+ case 5:
+ fSpecialTrigger=4; // different trigger class as MB
+ fTriggerSelectedManually = kTRUE;
+ break;
+ case 6:
+ fSpecialTrigger=4; // different trigger class as MB
+ fTriggerSelectedManually = kTRUE;
+ break;
+ case 7:
+ fSpecialTrigger=4; // different trigger class as MB
+ fTriggerSelectedManually = kTRUE;
+ break;
+ case 8:
+ fSpecialTrigger=4; // different trigger class as MB
+ fTriggerSelectedManually = kTRUE;
+ break;
+ case 9:
+ fSpecialTrigger=4; // different trigger class as MB
+ fTriggerSelectedManually = kTRUE;
+ break;
default:
AliError("Warning: Special Trigger Not known");
return kFALSE;
{ // Set Cut
switch (v0FinderType){
case 0: // on fly V0 finder
+ cout << "have chosen onfly V0" << endl;
fUseOnFlyV0Finder=kTRUE;
break;
case 1: // offline V0 finder
+ cout << "have chosen offline V0" << endl;
fUseOnFlyV0Finder=kFALSE;
break;
default:
fEtaCutMin = -0.1;
fLineCutZRSlopeMin = 0.;
break;
- case 1: // 1.2
- fEtaCut = 1.2;
+ case 1: // 1.2 // changed from 1.2 to 0.6 on 2013.06.10
+ fEtaCut = 0.6;
fLineCutZRSlope = tan(2*atan(exp(-fEtaCut)));
fEtaCutMin = -0.1;
fLineCutZRSlopeMin = 0.;
fEtaCutMin = -0.1;
fLineCutZRSlopeMin = 0.;
break;
- case 5: // 0.9 - 1.4
- fEtaCut = 1.4;
+ case 5: // 0.5
+ fEtaCut = 0.5;
fLineCutZRSlope = tan(2*atan(exp(-fEtaCut)));
- fEtaCutMin = 0.9;
- fLineCutZRSlopeMin = tan(2*atan(exp(-fEtaCutMin)));
+ fEtaCutMin = -0.1;
+ fLineCutZRSlopeMin = 0.;
break;
case 6: // 5.
fEtaCut = 5.;
fEtaCutMin = -0.1;
fLineCutZRSlopeMin = 0.;
break;
- case 7: // 0.1 - 0.8
- fEtaCut = 0.8;
+ case 7:
+ fEtaCut = 0.3;
fLineCutZRSlope = tan(2*atan(exp(-fEtaCut)));
- fEtaCutMin = 0.1;
- fLineCutZRSlopeMin = tan(2*atan(exp(-fEtaCutMin)));
+ fEtaCutMin = -0.1;
+ fLineCutZRSlopeMin = 0.;
break;
- case 8: // 0.1 - 0.8
- fEtaCut = 0.9;
+ // case 8: // 0.1 - 0.8
+ // fEtaCut = 0.9;
+ // fLineCutZRSlope = tan(2*atan(exp(-fEtaCut)));
+ // fEtaCutMin = 0.1;
+ // fLineCutZRSlopeMin = tan(2*atan(exp(-fEtaCutMin)));
+ // break;
+ case 8: // 0.4
+ fEtaCut = 0.4;
fLineCutZRSlope = tan(2*atan(exp(-fEtaCut)));
- fEtaCutMin = 0.1;
- fLineCutZRSlopeMin = tan(2*atan(exp(-fEtaCutMin)));
+ fEtaCutMin = -0.1;
+ fLineCutZRSlopeMin = 0.;
break;
case 9: // 10
fEtaCut = 10;
case 0: // 0
fMinClsTPC= 0.;
break;
- case 1: // 70
- fMinClsTPC= 70.;
+ case 1: // 60
+ fMinClsTPC= 60.;
break;
case 2: // 80
fMinClsTPC= 80.;
case 3: // 100
fMinClsTPC= 100.;
break;
- case 4: // 60% of findable clusters
- fMinClsTPCToF= 0.6;
- fUseCorrectedTPCClsInfo=0;
+ case 4: // 95% of findable clusters
+ fMinClsTPCToF= 0.95;
+ fUseCorrectedTPCClsInfo=1;
break;
case 5: // 0% of findable clusters
fMinClsTPCToF= 0.0;
AliCentrality *fESDCentrality=(AliCentrality*)esdEvent->GetCentrality();
if(fDetectorCentrality==0){
- return fESDCentrality->GetCentralityPercentile("V0M"); // default
+ if (fIsHeavyIon==2){
+ return fESDCentrality->GetCentralityPercentile("V0A"); // default for pPb
+ } else{
+ return fESDCentrality->GetCentralityPercentile("V0M"); // default
+ }
}
if(fDetectorCentrality==1){
return fESDCentrality->GetCentralityPercentile("CL1");
if(!fIsHeavyIon)return kTRUE;
if(fCentralityMin == 0 && fCentralityMax == 0) return kTRUE;//0-100%
- if(fCentralityMin >= fCentralityMax) return kTRUE;//0-100%
-
+ if(fCentralityMin >= fCentralityMax ){
+ if(fCentralityMax==0 && fModCentralityClass ==0) fCentralityMax=10; //CentralityRange = fCentralityMin-100%
+ else return kTRUE;//0-100%
+ }
Double_t centrality=GetCentrality(event);
if(centrality<0)return kFALSE;
// Use strict V0 amplitude cut for MC centrality
Float_t nv0amplitude = event->GetVZEROData()->GetMTotV0A()+event->GetVZEROData()->GetMTotV0C();
Float_t V0Amplitude10[10] = {9999999.0,13670,9345,6209,3944,2352,1272,611,255, 83};
- // 0 10 20 30 40 50 60 70 80 90%
+ // 0 10 20 30 40 50 60 70 80 90%
Float_t V0Amplitude5a[10] = {9999999.0,16612,13670,11290,9345,7650,6209,4984,3944,3074};
- // 0 5 10 15 20 25 30 35 40 45%
+ // 0 5 10 15 20 25 30 35 40 45%
Float_t V0Amplitude5b[10] = {3074,2352,1725,1272,899,611,402,255,152,83};
// 45 50 55 60 65 70 75 80 85 90%
-
+
if (fModCentralityClass == 3){
if(fMCEvent){
if(nv0amplitude > V0Amplitude10[fCentralityMax] && nv0amplitude <= V0Amplitude10[fCentralityMin])
Double_t fVertexZ=event->GetPrimaryVertex()->GetZ();
if(abs(fVertexZ)>fMaxVertexZ)return kFALSE;
+
+ if (fIsHeavyIon == 2){
+ if(fUtils->IsFirstEventInChunk(event)) return kFALSE;
+ if(!fUtils->IsVertexSelected2013pA(event)) return kFALSE;
+
+ }
+
return kTRUE;
}
///________________________________________________________________________
if(fESDEvent){
if (fESDEvent->GetPrimaryVertex() != NULL){
if(fESDEvent->GetPrimaryVertex()->GetNContributors()>0) {
+// cout << "accepted global" << fESDEvent->GetEventNumberInFile() << " with NCont: " << fESDEvent->GetPrimaryVertex()->GetNContributors() << endl;
return fESDEvent->GetPrimaryVertex()->GetNContributors();
}
}
if(fESDEvent->GetPrimaryVertexSPD() !=NULL){
if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
+// cout << "accepted SPD" << fESDEvent->GetEventNumberInFile() << " with NCont: " << fESDEvent->GetPrimaryVertexSPD()->GetNContributors() << endl;
return fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
} else {
AliWarning(Form("Number of contributors from bad vertex type:: %s",fESDEvent->GetPrimaryVertex()->GetName()));
+// cout << "rejected " << fESDEvent->GetEventNumberInFile() << endl;
return 0;
}
}
}
}
}
-
+ // cout << "rejected " << fESDEvent->GetEventNumberInFile() << endl;
return 0;
}
///________________________________________________________________________
-Bool_t AliConversionCuts::IsTriggerSelected()
+Bool_t AliConversionCuts::IsTriggerSelected(AliVEvent *fInputEvent)
{
AliInputEventHandler *fInputHandler=(AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
UInt_t isSelected = AliVEvent::kAny;
- if( fInputHandler && fInputHandler->GetEventSelection()) {
+ if (fInputHandler==NULL) return kFALSE;
+ if( fInputHandler->GetEventSelection() || fInputEvent->IsA()==AliAODEvent::Class()) {
+ if (!fTriggerSelectedManually){
+ if (fPreSelCut) fOfflineTriggerMask = AliVEvent::kAny;
+ else {
+ if (fIsHeavyIon == 1) fOfflineTriggerMask = AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral;
+ else if (fIsHeavyIon == 2) fOfflineTriggerMask = AliVEvent::kINT7;
+ else fOfflineTriggerMask = AliVEvent::kMB;
+ }
+ }
// Get the actual offline trigger mask for the event and AND it with the
// requested mask. If no mask requested select by default the event.
+// if (fPreSelCut) cout << "Trigger selected from outside: "<< fTriggerSelectedManually <<"\t Offline Trigger mask for Precut: " << fOfflineTriggerMask << endl;
+// else cout << "Trigger selected from outside: "<< fTriggerSelectedManually <<"\t Offline Trigger mask: " << fOfflineTriggerMask << endl;
+
if (fOfflineTriggerMask)
isSelected = fOfflineTriggerMask & fInputHandler->IsEventSelected();
}
// Fill Histogram
if(hTriggerClass){
- if (fIsSDDFired) hTriggerClass->Fill(32);
+ if (fIsSDDFired) hTriggerClass->Fill(33);
if (fInputHandler->IsEventSelected() & AliVEvent::kMB)hTriggerClass->Fill(0);
if (fInputHandler->IsEventSelected() & AliVEvent::kINT7)hTriggerClass->Fill(1);
if (fInputHandler->IsEventSelected() & AliVEvent::kMUON)hTriggerClass->Fill(2);
if (fInputHandler->IsEventSelected() & AliVEvent::kMuonUnlikeLowPt8)hTriggerClass->Fill(25);
if (fInputHandler->IsEventSelected() & AliVEvent::kMuonUnlikeLowPt0)hTriggerClass->Fill(26);
if (fInputHandler->IsEventSelected() & AliVEvent::kUserDefined)hTriggerClass->Fill(27);
- if (fInputHandler->IsEventSelected() & AliVEvent::kFastOnly)hTriggerClass->Fill(28);
- if (fInputHandler->IsEventSelected() & AliVEvent::kAnyINT)hTriggerClass->Fill(29);
- if (fInputHandler->IsEventSelected() & AliVEvent::kAny)hTriggerClass->Fill(30);
- if (!fInputHandler->IsEventSelected()) hTriggerClass->Fill(33);
+ if (fInputHandler->IsEventSelected() & AliVEvent::kTRD)hTriggerClass->Fill(28);
+ if (fInputHandler->IsEventSelected() & AliVEvent::kFastOnly)hTriggerClass->Fill(29);
+ if (fInputHandler->IsEventSelected() & AliVEvent::kAnyINT)hTriggerClass->Fill(30);
+ if (fInputHandler->IsEventSelected() & AliVEvent::kAny)hTriggerClass->Fill(31);
+ if (!fInputHandler->IsEventSelected()) hTriggerClass->Fill(34);
}
if(hTriggerClassSelected && isSelected){
- if (!fIsSDDFired) hTriggerClassSelected->Fill(32);
+ if (!fIsSDDFired) hTriggerClassSelected->Fill(33);
if (fInputHandler->IsEventSelected() & AliVEvent::kMB)hTriggerClassSelected->Fill(0);
if (fInputHandler->IsEventSelected() & AliVEvent::kINT7)hTriggerClassSelected->Fill(1);
if (fInputHandler->IsEventSelected() & AliVEvent::kMUON)hTriggerClassSelected->Fill(2);
if (fInputHandler->IsEventSelected() & AliVEvent::kMuonUnlikeLowPt8)hTriggerClassSelected->Fill(25);
if (fInputHandler->IsEventSelected() & AliVEvent::kMuonUnlikeLowPt0)hTriggerClassSelected->Fill(26);
if (fInputHandler->IsEventSelected() & AliVEvent::kUserDefined)hTriggerClassSelected->Fill(27);
- if (fInputHandler->IsEventSelected() & AliVEvent::kFastOnly)hTriggerClassSelected->Fill(28);
- if (fInputHandler->IsEventSelected() & AliVEvent::kAnyINT)hTriggerClassSelected->Fill(29);
- if (fInputHandler->IsEventSelected() & AliVEvent::kAny)hTriggerClassSelected->Fill(30);
+ if (fInputHandler->IsEventSelected() & AliVEvent::kTRD)hTriggerClassSelected->Fill(28);
+ if (fInputHandler->IsEventSelected() & AliVEvent::kFastOnly)hTriggerClassSelected->Fill(29);
+ if (fInputHandler->IsEventSelected() & AliVEvent::kAnyINT)hTriggerClassSelected->Fill(30);
+ if (fInputHandler->IsEventSelected() & AliVEvent::kAny)hTriggerClassSelected->Fill(31);
}
if(!isSelected)return kFALSE;
return kTRUE;
}
///________________________________________________________________________
-void AliConversionCuts::GetNotRejectedParticles(Int_t rejection, TList *HeaderList, AliMCEvent *MCEvent){
+void AliConversionCuts::GetNotRejectedParticles(Int_t rejection, TList *HeaderList, AliVEvent *MCEvent){
if(fNotRejectedStart){
delete[] fNotRejectedStart;
}
if(rejection == 0) return; // No Rejection
- AliGenCocktailEventHeader *cHeader = dynamic_cast<AliGenCocktailEventHeader*>(MCEvent->GenEventHeader());
- if(cHeader){
- TList *genHeaders = cHeader->GetHeaders();
+
+ AliGenCocktailEventHeader *cHeader = 0x0;
+ AliAODMCHeader *cHeaderAOD = 0x0;
+ Bool_t headerFound = kFALSE;
+
+ if(MCEvent->IsA()==AliMCEvent::Class()){
+ cHeader = dynamic_cast<AliGenCocktailEventHeader*>(dynamic_cast<AliMCEvent*>(MCEvent)->GenEventHeader());
+ if(cHeader) headerFound = kTRUE;
+ }
+ if(MCEvent->IsA()==AliAODEvent::Class()){
+ cHeaderAOD = dynamic_cast<AliAODMCHeader*>(MCEvent->FindListObject(AliAODMCHeader::StdBranchName()));
+ if(cHeaderAOD) headerFound = kTRUE;
+ }
+
+ if(headerFound){
+ TList *genHeaders = 0x0;
+ if(cHeader) genHeaders = cHeader->GetHeaders();
+ if(cHeaderAOD){
+ genHeaders = cHeaderAOD->GetCocktailHeaders();
+ if(genHeaders->GetEntries()==1){
+ SetRejectExtraSignalsCut(0);
+ return;
+ }
+ }
AliGenEventHeader* gh = 0;
fnHeaders = 0;
if(rejection == 1 || rejection == 3) fnHeaders = 1; // MinBiasHeader
fnHeaders = 1;
fNotRejectedStart[0] = 0;
- fNotRejectedEnd[0] = MCEvent->Stack()->GetNprimary()-1;
-// if(rejection == 2){
+ fNotRejectedEnd[0] = static_cast<AliMCEvent*>(MCEvent)->Stack()->GetNprimary()-1;
+ // if(rejection == 2){
fGeneratorNames = new TString[1];
fGeneratorNames[0] = "NoCocktailGeneratorFound";
// }
-
- AliGenPythiaEventHeader *mcHeaderPythia = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent->GenEventHeader());
+
+ AliGenPythiaEventHeader *mcHeaderPythia = dynamic_cast<AliGenPythiaEventHeader*>(static_cast<AliMCEvent*>(MCEvent)->GenEventHeader());
if (mcHeaderPythia) fGeneratorNames[0] = "NoCocktailGeneratorFound_Pythia";
- AliGenDPMjetEventHeader *mcHeaderPhojet = dynamic_cast<AliGenDPMjetEventHeader*>(MCEvent->GenEventHeader());
+ AliGenDPMjetEventHeader *mcHeaderPhojet = dynamic_cast<AliGenDPMjetEventHeader*>(static_cast<AliMCEvent*>(MCEvent)->GenEventHeader());
if (mcHeaderPhojet) fGeneratorNames[0] = "NoCocktailGeneratorFound_Phojet";
- AliGenHijingEventHeader *mcHeaderHijing = dynamic_cast<AliGenHijingEventHeader*>(MCEvent->GenEventHeader());
+ AliGenHijingEventHeader *mcHeaderHijing = dynamic_cast<AliGenHijingEventHeader*>(static_cast<AliMCEvent*>(MCEvent)->GenEventHeader());
if (mcHeaderHijing) fGeneratorNames[0] = "NoCocktailGeneratorFound_Hijing";
SetRejectExtraSignalsCut(0);
}
//_________________________________________________________________________
-Int_t AliConversionCuts::IsParticleFromBGEvent(Int_t index, AliStack *MCStack){
+Int_t AliConversionCuts::IsParticleFromBGEvent(Int_t index, AliStack *MCStack, AliVEvent *InputEvent){
// Not Accepted == kFALSE == 0
// Accepted == kTRUE == 1
// FirstHeader == kTRUE == 3
-
if(index < 0) return 0; // No Particle
Int_t accepted = 0;
- if( index >= MCStack->GetNprimary()){ // Secondary Particle
- if( ((TParticle*)MCStack->Particle(index))->GetMother(0) < 0) return 1; // Secondary Particle without Mother??
- return IsParticleFromBGEvent(((TParticle*)MCStack->Particle(index))->GetMother(0),MCStack);
+ if(!InputEvent || InputEvent->IsA()==AliESDEvent::Class()){
+ if( index >= MCStack->GetNprimary()){ // Secondary Particle
+ if( ((TParticle*)MCStack->Particle(index))->GetMother(0) < 0) return 1; // Secondary Particle without Mother??
+ return IsParticleFromBGEvent(((TParticle*)MCStack->Particle(index))->GetMother(0),MCStack,InputEvent);
+ }
+ for(Int_t i = 0;i<fnHeaders;i++){
+ if(index >= fNotRejectedStart[i] && index <= fNotRejectedEnd[i]){
+ accepted = 1;
+ if(i == 0) accepted = 2; // MB Header
+ }
+ }
}
- for(Int_t i = 0;i<fnHeaders;i++){
- if(index >= fNotRejectedStart[i] && index <= fNotRejectedEnd[i]){
- accepted = 1;
- if(i == 0) accepted = 2; // MB Header
+ else if(InputEvent->IsA()==AliAODEvent::Class()){
+ TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(InputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
+ AliAODMCParticle *aodMCParticle = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(index));
+ if(!aodMCParticle->IsPrimary()){
+ if( aodMCParticle->GetMother() < 0) return 1;// Secondary Particle without Mother??
+ return IsParticleFromBGEvent(aodMCParticle->GetMother(),MCStack,InputEvent);
+ }
+ index = abs(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(index))->GetLabel());
+ for(Int_t i = 0;i<fnHeaders;i++){
+ if(index >= fNotRejectedStart[i] && index <= fNotRejectedEnd[i]){
+ accepted = 1;
+ if(i == 0) accepted = 2; // MB Header
+ }
}
}
//_________________________________________________________________________
Int_t AliConversionCuts::IsEventAcceptedByConversionCut(AliConversionCuts *ReaderCuts, AliVEvent *InputEvent, AliMCEvent *MCEvent, Bool_t isHeavyIon){
+ if ( !IsTriggerSelected(InputEvent) )
+ return 3;
+
if(isHeavyIon && !(IsCentralitySelected(InputEvent,MCEvent)))
return 1; // Check Centrality --> Not Accepted => eventQuality = 1
if(!isHeavyIon && GetIsFromPileup()){
if(InputEvent->IsPileupFromSPD(3,0.8,3.,2.,5.)){
+
return 6; // Check Pileup --> Not Accepted => eventQuality = 6
}
}
-
+
Bool_t hasV0And = ReaderCuts->HasV0AND();
Bool_t isSDDFired = ReaderCuts->IsSDDFired();
if( (IsSpecialTrigger() == 2 || IsSpecialTrigger() == 3) && !isSDDFired && !MCEvent)
if( (IsSpecialTrigger() == 1 || IsSpecialTrigger() == 3) && !hasV0And)
return 8; // V0AND requested but no fired
-
+
return 0;
}
-
-Float_t AliConversionCuts::GetWeightForMeson(TString period, Int_t index, AliStack *MCStack){
- if (!(period.CompareTo("LHC12f1a") == 0 || period.CompareTo("LHC12f1b") == 0 || period.CompareTo("LHC12i3") == 0 )) return 1.;
+//_________________________________________________________________________
+Float_t AliConversionCuts::GetWeightForMeson(TString period, Int_t index, AliStack *MCStack, AliVEvent *InputEvent){
+ if (!(period.CompareTo("LHC12f1a") == 0 || period.CompareTo("LHC12f1b") == 0 || period.CompareTo("LHC12i3") == 0 || period.CompareTo("LHC11a10a") == 0 || period.CompareTo("LHC11a10b") == 0 || period.CompareTo("LHC11a10b_bis") == 0 || period.CompareTo("LHC11a10a_bis") == 0 || period.CompareTo("LHC11a10b_plus") == 0 || period.CompareTo("LHC13d2") == 0)) return 1.;
Int_t kCaseGen = 0;
for (Int_t i = 0; i < fnHeaders; i++){
if (index >= fNotRejectedStart[i] && index < fNotRejectedEnd[i]+1){
-// cout << fGeneratorNames[i].Data() << endl;
if (fGeneratorNames[i].CompareTo("Pythia") == 0){
kCaseGen = 1;
} else if (fGeneratorNames[i].CompareTo("DPMJET") == 0){
kCaseGen = 2;
- } else if (fGeneratorNames[i].CompareTo("HIJING") == 0){
+ } else if (fGeneratorNames[i].CompareTo("HIJING") == 0 ||
+ fGeneratorNames[i].CompareTo("Hijing") == 0 ||
+ fGeneratorNames[i].Contains("hijing")){
kCaseGen = 3;
} else if (fGeneratorNames[i].CompareTo("BOX") == 0){
kCaseGen = 4;
kCaseGen = 2;
} else if (fGeneratorNames[i].CompareTo("NoCocktailGeneratorFound_Hijing") == 0){
kCaseGen = 3;
- }
-// cout << "resulting kCaseGen :" << kCaseGen << endl;
+ }
}
}
+ if (kCaseGen == 0) return 1;
- Double_t mesonPt = ((TParticle*)MCStack->Particle(index))->Pt();
-
-// Double_t mesonY = 10.;
-// if(((TParticle*)MCStack->Particle(index))->Energy() - ((TParticle*)MCStack->Particle(index))->Pz() == 0 || ((TParticle*)MCStack->Particle(index))->Energy() + ((TParticle*)MCStack->Particle(index))->Pz() == 0){
-// mesonY=10.;
-// } else{
-// mesonY = 0.5*(TMath::Log((((TParticle*)MCStack->Particle(index))->Energy()+((TParticle*)MCStack->Particle(index))->Pz()) / (((TParticle*)MCStack->Particle(index))->Energy()-((TParticle*)MCStack->Particle(index))->Pz())));
-// }
- Double_t mesonMass = ((TParticle*)MCStack->Particle(index))->GetCalcMass();
- Float_t functionResult = 1.;
- if (kCaseGen == 1){
+
+ Double_t mesonPt = 0;
+ Double_t mesonMass = 0;
+ Int_t PDGCode = 0;
+ if(!InputEvent || InputEvent->IsA()==AliESDEvent::Class()){
+ mesonPt = ((TParticle*)MCStack->Particle(index))->Pt();
+ mesonMass = ((TParticle*)MCStack->Particle(index))->GetCalcMass();
+ PDGCode = ((TParticle*)MCStack->Particle(index))->GetPdgCode();
+ }
+ else if(InputEvent->IsA()==AliAODEvent::Class()){
+ TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(InputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
+ AliAODMCParticle *aodMCParticle = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(index));
+ mesonPt = aodMCParticle->Pt();
+ mesonMass = aodMCParticle->GetCalcMass();
+ PDGCode = aodMCParticle->GetPdgCode();
+ }
+
+ Float_t functionResultMC = 1.;
+ if (kCaseGen == 1){ // Pythia 6
Float_t dNdyMC = 2.1462;
Float_t nMC = 7.06055;
Float_t tMC = 0.12533;
- if ( ((TParticle*)MCStack->Particle(index))->GetPdgCode() == 111){
+ if ( PDGCode == 111){
dNdyMC = 2.1462;
nMC = 7.06055;
tMC = 0.12533;
- } else if ( ((TParticle*)MCStack->Particle(index))->GetPdgCode() == 221){
+ } else if ( PDGCode == 221){
dNdyMC = 0.2357;
nMC = 5.9105;
tMC = 0.1525;
}
- functionResult = dNdyMC / ( 2 * TMath::Pi())*(nMC-1.)*(nMC-2.) / (nMC*tMC*(nMC*tMC+mesonMass*(nMC-2.))) * pow(1.+(sqrt(mesonPt*mesonPt+mesonMass*mesonMass)-mesonMass)/(nMC*tMC), -nMC);
- } else if (kCaseGen == 2){
+ functionResultMC = dNdyMC / ( 2 * TMath::Pi())*(nMC-1.)*(nMC-2.) / (nMC*tMC*(nMC*tMC+mesonMass*(nMC-2.))) * TMath::Power(1.+(TMath::Sqrt(mesonPt*mesonPt+mesonMass*mesonMass)-mesonMass)/(nMC*tMC), -nMC);
+ } else if (kCaseGen == 2){ // Phojet
Float_t dNdyMC = 2.35978;
Float_t nMC = 6.81795;
Float_t tMC = 0.11492;
- if ( ((TParticle*)MCStack->Particle(index))->GetPdgCode() == 111){
+ if ( PDGCode == 111){
dNdyMC = 2.35978;
nMC = 6.81795;
tMC = 0.11492;
- } else if ( ((TParticle*)MCStack->Particle(index))->GetPdgCode() == 221){
+ } else if ( PDGCode == 221){
dNdyMC = 0.3690;
nMC = 5.55809;
tMC = 0.13387;
}
- functionResult = dNdyMC / ( 2 * TMath::Pi())*(nMC-1.)*(nMC-2.) / (nMC*tMC*(nMC*tMC+mesonMass*(nMC-2.))) * pow(1.+(sqrt(mesonPt*mesonPt+mesonMass*mesonMass)-mesonMass)/(nMC*tMC), -nMC);
- } else if (kCaseGen == 4){
-// functionResult = 1./sqrt(1.-mesonMass*mesonMass/((mesonMass*mesonMass+mesonPt*mesonPt)*cosh(mesonY)*cosh(mesonY)));
+ functionResultMC = dNdyMC / ( 2 * TMath::Pi())*(nMC-1.)*(nMC-2.) / (nMC*tMC*(nMC*tMC+mesonMass*(nMC-2.))) * TMath::Power(1.+(TMath::Sqrt(mesonPt*mesonPt+mesonMass*mesonMass)-mesonMass)/(nMC*tMC), -nMC);
+ } else if (kCaseGen == 4){ // BOX generators pp
+// functionResultMC = 1./sqrt(1.-mesonMass*mesonMass/((mesonMass*mesonMass+mesonPt*mesonPt)*cosh(mesonY)*cosh(mesonY)));
Float_t a = 0.23437;
Float_t b = 5.6661;
Float_t c = -1430.5863;
Float_t d = -0.6966624;
Float_t e = 252.3742;
- if ( ((TParticle*)MCStack->Particle(index))->GetPdgCode() == 111){
+ if ( PDGCode == 111){
a = 0.23437;
b = 5.6661;
c = -1430.5863;
d = -0.6966624;
e = 252.3742;
- } else if ( ((TParticle*)MCStack->Particle(index))->GetPdgCode() == 221){
+ } else if ( PDGCode == 221){
a = 0.10399;
b = 4.35311;
c = -12.17723;
d = -0.01172;
e =1.85140;
}
- functionResult = a*TMath::Power(mesonPt,-1.*(b+c/(TMath::Power(mesonPt,d)+e)))*1./mesonPt *1./1.6 *1./(2.* TMath::Pi());
+ functionResultMC = a*TMath::Power(mesonPt,-1.*(b+c/(TMath::Power(mesonPt,d)+e)))*1./mesonPt *1./1.6 *1./(2.* TMath::Pi());
+ } else if (kCaseGen == 3 ){ // HIJING
+ if ( PDGCode == 111 && fDoReweightHistoMCPi0 && hReweightMCHistPi0!= 0x0){
+ functionResultMC = hReweightMCHistPi0->Interpolate(mesonPt);
+ }
+ if ( PDGCode == 310 && fDoReweightHistoMCK0s && hReweightMCHistK0s!= 0x0){
+ functionResultMC = hReweightMCHistK0s->Interpolate(mesonPt);
+ }
}
-
- Float_t dNdyData = 2.2328;
- Float_t nData = 7.1473;
- Float_t tData = 0.1346;
- if ( ((TParticle*)MCStack->Particle(index))->GetPdgCode() == 111){
+
+ Float_t functionResultData = 1;
+ if (kCaseGen == 1 || kCaseGen == 2 || kCaseGen == 4 ){
+ Float_t dNdyData = 2.2328;
+ Float_t nData = 7.1473;
+ Float_t tData = 0.1346;
+ if ( PDGCode == 111){
dNdyData = 2.2328;
nData = 7.1473;
tData = 0.1346;
- } else if ( ((TParticle*)MCStack->Particle(index))->GetPdgCode() == 221){
+ } else if ( PDGCode == 221){
dNdyData = 0.38992; //be careful this fit is not optimal, eta in data still has problems
nData = 5.72778;
tData = 0.13835;
}
- Float_t tsallisData = dNdyData / ( 2 * TMath::Pi())*(nData-1.)*(nData-2.) / (nData*tData*(nData*tData+mesonMass*(nData-2.))) * pow(1.+(sqrt(mesonPt*mesonPt+mesonMass*mesonMass)-mesonMass)/(nData*tData), -nData);
-// cout << "tsallisData/functionResult: " << tsallisData/functionResult << endl;
- return tsallisData/functionResult;
-// return
+ functionResultData = dNdyData / ( 2 * TMath::Pi())*(nData-1.)*(nData-2.) / (nData*tData*(nData*tData+mesonMass*(nData-2.))) * TMath::Power(1.+(TMath::Sqrt(mesonPt*mesonPt+mesonMass*mesonMass)-mesonMass)/(nData*tData), -nData);
+ } else {
+ Float_t a = 0.;
+ Float_t b = 0.;
+ Float_t c = 0.;
+ Float_t d = 0.;
+ Float_t e = 0.;
+ if ( PDGCode == 111 ){
+ if (fModCentralityClass == 1 && fCentralityMin == 0 && fCentralityMax == 1 ){ // 0-5 % PbPb
+ a = 25.8747458223;
+ b = 5.8761820045;
+ c = -33.9928191673;
+ d = 3.0731850142;
+ e = 13.2500447620;
+ } else if (fModCentralityClass == 1 && fCentralityMin == 1 && fCentralityMax == 2){ // 5-10% PbPb
+ a = 21.7518148922;
+ b = 5.8441200081;
+ c = -17.1497051691;
+ d = 2.3799090842;
+ e = 5.4346404718;
+ } else if (fModCentralityClass == 0 && fCentralityMin == 0 && fCentralityMax == 1){ // 0-10% PbPb
+ a = 22.9852133622;
+ b = 5.8602063916;
+ c = -17.0992478654;
+ d = 2.4426218039;
+ e = 5.1194526345;
+ } else if (fModCentralityClass == 0 && fCentralityMin == 1 && fCentralityMax == 2){ // 10-20% PbPb
+ a = 19.3237333776;
+ b = 5.8145906958;
+ c = -13.8316665424;
+ d = 2.3737630637;
+ e = 4.7690300693;
+ } else if (fModCentralityClass == 0 && fCentralityMin == 2 && fCentralityMax == 4){ // 20-40% PbPb
+ a = 11.2656032751;
+ b = 5.8003194354;
+ c = -13.3936105929;
+ d = 2.3371452334;
+ e = 4.4726244958;
+ } else if (fModCentralityClass == 0 && fCentralityMin == 4 && fCentralityMax == 6){ // 40-60% PbPb
+ a = 4.1578154081;
+ b = 5.6450005163;
+ c = -8.4309375240;
+ d = 1.8918308704;
+ e = 2.9429194709;
+ } else if (fModCentralityClass == 0 && fCentralityMin == 6 && fCentralityMax == 8){ // 60-80% PbPb
+ a = 1.0635443810;
+ b = 5.1337469970;
+ c = -8.5906997238;
+ d = 2.9794995997;
+ e = 3.9294980048;
+ } else if (fModCentralityClass == 0 && fCentralityMin == 0 && fCentralityMax == 2){ // 0-20% PbPb
+ a = 21.7018745556;
+ b = 5.9019352094;
+ c = -14.2295510326;
+ d = 2.2104490688;
+ e = 4.2969671500;
+ } else if (fModCentralityClass == 0 && fCentralityMin == 0 && fCentralityMax == 4){ // 0-40% PbPb
+ a = 16.8227412106;
+ b = 5.8660502207;
+ c = -12.0978551215;
+ d = 2.1695068981;
+ e = 3.5349621182;
+ } else if (fModCentralityClass == 0 && fCentralityMin == 0 && fCentralityMax == 8){ // 0-80% PbPb
+ a = 9.4675681080;
+ b = 5.8114944205;
+ c = -10.4901523616;
+ d = 2.0607982712;
+ e = 2.9262259130;
+ } else if (fModCentralityClass == 0 && fCentralityMin == 4 && fCentralityMax == 8){ // 60-80% PbPb
+ a = 2.5985551785;
+ b = 5.4118895738;
+ c = -8.2510958428;
+ d = 2.2551249190;
+ e = 3.0700919491;
+ }
+
+ functionResultData = a*TMath::Power(mesonPt,-1*(b+c/(TMath::Power(mesonPt,d)+e)));
+ }
+
+ }
+
+ Double_t weight = 1;
+ if (PDGCode == 111 || PDGCode == 221){
+ if (functionResultData != 0. && functionResultMC != 0. && isfinite(functionResultData) && isfinite(functionResultMC)){
+ weight = functionResultData/functionResultMC;
+ if ( !(kCaseGen == 3 && fDoReweightHistoMCPi0 && hReweightMCHistPi0!= 0x0 && PDGCode == 111)){
+ weight = 1.;
+ }
+ if (!isfinite(functionResultData)) weight = 1.;
+ if (!isfinite(weight)) weight = 1.;
+ }
+ } else if (PDGCode == 310 && functionResultMC != 0 && isfinite(functionResultMC)){
+ weight = functionResultMC;
+ }
+
+// if (fModCentralityClass == 0 && fCentralityMin == 4 && fCentralityMax == 6 && PDGCode == 111){
+// cout << period.Data() << "\t" << kCaseGen << "\t" <<fModCentralityClass<< "\t" <<fCentralityMin<< "\t" <<fCentralityMax << "\t" << mesonPt << "\t" <<mesonMass<< "\t"<<functionResultData << "\t"<< functionResultMC << "\t" << weight <<endl;
+// }
+ return weight;
}
-
///________________________________________________________________________
AliConversionCuts* AliConversionCuts::GetStandardCuts2010PbPb(){
//Create and return standard 2010 PbPb cuts
cout<<"Warning: Initialization of Standardcuts2010pp failed"<<endl;}
return cuts;
}
-
+///________________________________________________________________________
+void AliConversionCuts::GetCorrectEtaShiftFromPeriod(TString periodName){
+
+ if(periodName.CompareTo("LHC12g") == 0 || //pilot run 2012
+ periodName.CompareTo("LHC13b") == 0 || //mainly minimum bias
+ periodName.CompareTo("LHC13c") == 0 || //mainly minimum bias
+ periodName.CompareTo("LHC13d") == 0 || //mainly triggered
+ periodName.CompareTo("LHC13e") == 0 || //mainly triggered
+ periodName.CompareTo("LHC13c3") == 0 || //MC Starlight, anchor LHC13d+e
+ periodName.CompareTo("LHC13c2") == 0 || //MC Starlight, coherent J/Psi, UPC muon anchor LHC13d+e
+ periodName.CompareTo("LHC13b4") == 0 || //MC Pythia 6 (Jet-Jet), anchor LHC13b
+ periodName.CompareTo("LHC13b2_fix_1") == 0 || //MC DPMJET, anchr LHC13b+c
+ periodName.CompareTo("LHC13b3") == 0 || //MC HIJING, weighted to number of events per run, anchor LHC13b
+ periodName.CompareTo("LHC13b2") == 0 || // MC DPMJET, wrong energy, anchor LHC13b
+ periodName.CompareTo("LHC13b2_plus") == 0 || // MC DPMJET, weighted to number event per run, anchor LHC13b
+ periodName.CompareTo("LHC13c1_bis") == 0 || // MC AMPT fast generation, pT hardbin, anchor ?
+ periodName.CompareTo("LHC13c1") == 0 || // MC AMPT fast generation, anchor ?
+ periodName.CompareTo("LHC13b1") == 0 || // MC DPMJET, fragments, with fixed label 0, anchor LHC12g
+ periodName.CompareTo("LHC12g4b_fix") == 0 || // MC DPMJET, with fixed label 0, anchor LHC12g
+ periodName.CompareTo("LHC12g1_fix") == 0 || // MC ?, with fixed label 0, anchor LHC12g
+ periodName.CompareTo("LHC12g4c") == 0 || // MC DPMJET, shifted vertex runs, anchor LHC12g
+ periodName.CompareTo("LHC12h6") == 0 || // MC muon cocktail, anchor LHC12g
+ periodName.CompareTo("LHC12g4b") == 0 || // MC DPMJET 3rd iteration, anchor LHC12g
+ periodName.CompareTo("LHC12g4a") == 0 || // MC DPMJET improved, anchor LHC12g
+ periodName.CompareTo("LHC12g4") == 0 || // MC DPMJET, anchor LHC12g
+ periodName.CompareTo("LHC12g5") == 0 || // MC PHOJET, anchor LHC12g
+ periodName.CompareTo("LHC12g2") == 0 || // MC Starlight background, anchor LHC12g
+ periodName.CompareTo("LHC12g1") == 0 ) // MC ?, anchor LHC12g
+ {
+ printf(" Gamma Conversion Cuts %s :: pPb Run doing Eta Shift of %f \n\n",(GetCutNumber()).Data(),-0.465);
+ SetEtaShift(-0.465);
+ }
+ else if(periodName.CompareTo("LHC13f") == 0 ||
+ periodName.CompareTo("LHC13c6b") == 0 ||// MC Jpsi -> mumu, anchor LHC13f
+ periodName.CompareTo("LHC13c5") == 0 || //MC Starlight, gamma gamma UPC muon, anchor LHC13f
+ periodName.CompareTo("LHC13c4") == 0 )//MC Starlight, coherent JPsi, UPC muon, anchor LHC13f
+ {
+ printf(" Gamma Conversion Cuts %s :: Pbp Run doing Eta Shift of %f \n\n",(GetCutNumber()).Data(),0.465);
+ SetEtaShift(+0.465);
+ }
+ else printf(" Gamma Conversion Cuts %s :: Automatic Eta Shift requested but Period is not known -> No Shift \n\n",(GetCutNumber()).Data());
+}