X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWG0%2FAliPWG0Helper.cxx;h=4ce88d4b4ca35a26887c2513942880d2cb34a70f;hb=540898293f7ba3c3896ecbd1a2c8906a8540aa21;hp=bd5a9ad6a7d8b0e5d8b9c2341de0d00bed3a22e6;hpb=7a11141c871fc9d4b21c9c63e7caa1c304f76afb;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWG0/AliPWG0Helper.cxx b/PWG0/AliPWG0Helper.cxx index bd5a9ad6a7d..4ce88d4b4ca 100644 --- a/PWG0/AliPWG0Helper.cxx +++ b/PWG0/AliPWG0Helper.cxx @@ -19,6 +19,7 @@ #include #include #include +#include #include #include @@ -28,87 +29,10 @@ #include #include -#include - //____________________________________________________________________ ClassImp(AliPWG0Helper) Int_t AliPWG0Helper::fgLastProcessType = -1; -AliOfflineTrigger* AliPWG0Helper::fgOfflineTrigger = 0; - -//____________________________________________________________________ -AliOfflineTrigger* AliPWG0Helper::GetOfflineTrigger() -{ - // returns the current AliOfflineTrigger object - // creates one if needed - - if (!fgOfflineTrigger) - fgOfflineTrigger = new AliOfflineTrigger; - - return fgOfflineTrigger; -} - -//____________________________________________________________________ -Bool_t AliPWG0Helper::IsEventTriggered(const AliESDEvent* aEsd, Trigger trigger) -{ - // checks if an event has been triggered - - if (trigger & kOfflineFlag) - return GetOfflineTrigger()->IsEventTriggered(aEsd, trigger); - - return IsEventTriggered(aEsd->GetTriggerMask(), trigger); -} - -//____________________________________________________________________ -Bool_t AliPWG0Helper::IsEventTriggered(ULong64_t triggerMask, Trigger trigger) -{ - // check if the event was triggered - // - // this function needs the branch fTriggerMask - - // definitions from p-p.cfg - ULong64_t spdFO = (1 << 14); - ULong64_t v0left = (1 << 10); - ULong64_t v0right = (1 << 11); - - switch (trigger) - { - case kAcceptAll: - { - return kTRUE; - break; - } - case kMB1: - { - if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right))) - return kTRUE; - break; - } - case kMB2: - { - if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right))) - return kTRUE; - break; - } - case kMB3: - { - if (triggerMask & spdFO && (triggerMask & v0left) && (triggerMask & v0right)) - return kTRUE; - break; - } - case kSPDGFO: - { - if (triggerMask & spdFO) - return kTRUE; - break; - } - default: - Printf("IsEventTriggered: ERROR: Trigger type %d not implemented in this method", (Int_t) trigger); - break; - } - - return kFALSE; -} //____________________________________________________________________ Bool_t AliPWG0Helper::TestVertex(const AliESDVertex* vertex, AnalysisMode analysisMode, Bool_t debug) @@ -118,7 +42,8 @@ Bool_t AliPWG0Helper::TestVertex(const AliESDVertex* vertex, AnalysisMode analys Float_t requiredZResolution = -1; if (analysisMode & kSPD || analysisMode & kTPCITS) { - requiredZResolution = 0.1; + // disable cut on resolution + requiredZResolution = 1000; } else if (analysisMode & kTPC) requiredZResolution = 10.; @@ -131,12 +56,22 @@ Bool_t AliPWG0Helper::TestVertex(const AliESDVertex* vertex, AnalysisMode analys Printf("AliPWG0Helper::TestVertex: Resolution too poor %f (required: %f", zRes, requiredZResolution); return kFALSE; } + + if (vertex->IsFromVertexerZ()) + { + if (vertex->GetDispersion() > 0.02) + { + if (debug) + Printf("AliPWG0Helper::TestVertex: Delta Phi too large in Vertexer Z: %f (required: %f", vertex->GetDispersion(), 0.02); + return kFALSE; + } + } return kTRUE; } //____________________________________________________________________ -const AliESDVertex* AliPWG0Helper::GetVertex(AliESDEvent* aEsd, AnalysisMode analysisMode, Bool_t debug, Bool_t bRedoTPC) +const AliESDVertex* AliPWG0Helper::GetVertex(const AliESDEvent* aEsd, AnalysisMode analysisMode, Bool_t debug) { // Get the vertex from the ESD and returns it if the vertex is valid // @@ -144,32 +79,29 @@ const AliESDVertex* AliPWG0Helper::GetVertex(AliESDEvent* aEsd, AnalysisMode ana // also the quality criteria that are applied) const AliESDVertex* vertex = 0; - if (analysisMode & kSPD || analysisMode & kTPCITS) + if (analysisMode & kSPD) { vertex = aEsd->GetPrimaryVertexSPD(); if (debug) Printf("AliPWG0Helper::GetVertex: Returning SPD vertex"); } - else if (analysisMode & kTPC) + else if (analysisMode & kTPCITS) { - if(bRedoTPC){ + vertex = aEsd->GetPrimaryVertexTracks(); + if (debug) + Printf("AliPWG0Helper::GetVertex: Returning vertex from tracks"); + if (vertex && vertex->GetNContributors() <= 0) + { if (debug) - Printf("AliPWG0Helper::GetVertex: Redoing vertex"); - Double_t kBz = aEsd->GetMagneticField(); - AliVertexerTracks vertexer(kBz); - vertexer.SetTPCMode(); - AliESDVertex *vTPC = vertexer.FindPrimaryVertex(aEsd); - aEsd->SetPrimaryVertexTPC(vTPC); - for (Int_t i=0; iGetNumberOfTracks(); i++) { - AliESDtrack *t = aEsd->GetTrack(i); - t->RelateToVertexTPC(vTPC, kBz, kVeryBig); - } - delete vTPC; + Printf("AliPWG0Helper::GetVertex: Vertex from tracks has no contributors. Falling back to SPD vertex."); + vertex = aEsd->GetPrimaryVertexSPD(); } - + } + else if (analysisMode & kTPC) + { vertex = aEsd->GetPrimaryVertexTPC(); if (debug) - Printf("AliPWG0Helper::GetVertex: Returning vertex from tracks"); + Printf("AliPWG0Helper::GetVertex: Returning vertex from TPC-only tracks"); } else Printf("AliPWG0Helper::GetVertex: ERROR: Invalid second argument %d", analysisMode); @@ -437,13 +369,13 @@ AliPWG0Helper::MCProcessType AliPWG0Helper::GetDPMjetEventProcessType(AliGenEven } - if(dpmJetType == 1){ // this is explicitly inelastic + if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction globalType = kND; } - else if(dpmJetType==5||dpmJetType==6){ + else if (dpmJetType==5 || dpmJetType==6) { globalType = kSD; } - else if(dpmJetType==7||dpmJetType==4){// DD and double pomeron + else if (dpmJetType==7) { globalType = kDD; } return globalType; @@ -581,41 +513,7 @@ void AliPWG0Helper::NormalizeToBinWidth(TH2* hist) } //____________________________________________________________________ -const char* AliPWG0Helper::GetTriggerName(Trigger trigger) -{ - // returns the name of the requested trigger - // the returned string will only be valid until the next call to this function [not thread-safe] - - static TString str; - - UInt_t triggerNoFlags = (UInt_t) trigger % (UInt_t) kStartOfFlags; - - switch (triggerNoFlags) - { - case kAcceptAll : str = "ACCEPT ALL (bypass!)"; break; - case kMB1 : str = "MB1"; break; - case kMB2 : str = "MB2"; break; - case kMB3 : str = "MB3"; break; - case kSPDGFO : str = "SPD GFO"; break; - case kV0A : str = "V0 A"; break; - case kV0C : str = "V0 C"; break; - case kZDC : str = "ZDC"; break; - case kZDCA : str = "ZDC A"; break; - case kZDCC : str = "ZDC C"; break; - case kFMDA : str = "FMD A"; break; - case kFMDC : str = "FMD C"; break; - case kFPANY : str = "SPD GFO | V0 | ZDC | FMD"; break; - default: str = ""; break; - } - - if (trigger & kOfflineFlag) - str += " OFFLINE"; - - return str; -} - -//____________________________________________________________________ -void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, Trigger trigger) +void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, AliTriggerAnalysis::Trigger trigger, AliPWG0Helper::DiffTreatment diffTreatment) { // // Prints the given configuration @@ -626,6 +524,9 @@ void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, Trigger trigger) if (analysisMode & kSPD) str += "SPD-only"; + if (analysisMode & kSPDOnlyL0) + str += " (only L0 clusters)"; + if (analysisMode & kTPC) str += "TPC-only"; @@ -640,9 +541,190 @@ void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, Trigger trigger) str += " (WITHOUT field)"; str += "< and trigger >"; - str += GetTriggerName(trigger); + str += AliTriggerAnalysis::GetTriggerName(trigger); + str += "< and diffractive treatment >"; + + switch (diffTreatment) + { + case kMCFlags: + str += "MC flags"; + break; + + case kUA5Cuts: + str += "UA5 cuts"; + break; + + case kE710Cuts: + str += "E710 cuts"; + break; + + case kALICEHadronLevel: + str += "ALICE hadron level"; + break; + } + str += "< <<<<"; Printf("%s", str.Data()); } +//____________________________________________________________________ +Double_t AliPWG0Helper::Rapidity(Double_t pt, Double_t pz, Double_t m) +{ + // + // calculates rapidity keeping the sign in case E == pz + // + + Double_t energy = TMath::Sqrt(pt*pt+pz*pz+m*m); + if (energy != TMath::Abs(pz)) + return 0.5*TMath::Log((energy+pz)/(energy-pz)); + + Printf("W- mt=0"); + return TMath::Sign(1.e30,pz); +} + +//____________________________________________________________________ +Bool_t AliPWG0Helper::IsHadronLevelSingleDiffractive(AliStack* stack, Float_t cms, Float_t xiMin, Float_t xiMax) +{ + // + // return if process is single diffractive on hadron level + // + // xiMax and xiMin cut on M^2/s + // + // Based on code from Martin Poghoysan + // + + TParticle* part1 = 0; + TParticle* part2 = 0; + + Double_t smallestY = 1e10; + Double_t largestY = -1e10; + + for (Int_t iParticle = 0; iParticle < stack->GetNprimary(); iParticle++) + { + TParticle* part = stack->Particle(iParticle); + if (!part) + continue; + + Int_t pdg = TMath::Abs(part->GetPdgCode()); + + Int_t child1 = part->GetFirstDaughter(); + Int_t ist = part->GetStatusCode(); + + Int_t mfl = Int_t (pdg / TMath::Power(10, Int_t(TMath::Log10(pdg)))); + if (child1 > -1 || ist != 1) + mfl = 0; // select final state charm and beauty + if (!(stack->IsPhysicalPrimary(iParticle) || pdg == 111 || pdg == 3212 || pdg==3124 || mfl >= 4)) + continue; + Int_t imother = part->GetFirstMother(); + if (imother>0) + { + TParticle *partM = stack->Particle(imother); + Int_t pdgM=TMath::Abs(partM->GetPdgCode()); + if (pdgM==111 || pdgM==3124 || pdgM==3212) + continue; + } + + Double_t y = 0; + + // fix for problem with getting mass of particle 3124 + if (pdg != 3124) + y = Rapidity(part->Pt(), part->Pz(), part->GetMass()); + else + y = Rapidity(part->Pt(), part->Pz(), 1.5195); + + if (y < smallestY) + { + smallestY = y; + part1 = part; + } + + if (y > largestY) + { + largestY = y; + part2 = part; + } + } + + if (part1 == 0 || part2 == 0) + return kFALSE; + + Int_t pdg1 = part1->GetPdgCode(); + Int_t pdg2 = part2->GetPdgCode(); + + Double_t pt1 = part1->Pt(); + Double_t pt2 = part2->Pt(); + Double_t pz1 = part1->Pz(); + Double_t pz2 = part2->Pz(); + + Double_t y1 = TMath::Abs(Rapidity(pt1, pz1, 0.938)); + Double_t y2 = TMath::Abs(Rapidity(pt2, pz2, 0.938)); + + Int_t arm = -99999; + if (pdg1 == 2212 && pdg2 == 2212) + { + if (y1 > y2) + arm = 0; + else + arm = 1; + } + else if (pdg1 == 2212) + arm = 0; + else if (pdg2 == 2212) + arm = 1; + + Double_t M02s = 1. - 2 * part1->Energy() / cms; + Double_t M12s = 1. - 2 * part2->Energy() / cms; + + if (arm == 0 && M02s > xiMin && M02s < xiMax) + return kTRUE; + else if (arm == 1 && M12s > xiMin && M12s < xiMax) + return kTRUE; + + return kFALSE; +} + +//____________________________________________________________________ +AliPWG0Helper::MCProcessType AliPWG0Helper::GetEventProcessType(AliESDEvent* esd, AliHeader* header, AliStack* stack, AliPWG0Helper::DiffTreatment diffTreatment) +{ + // + // get process type + // + // diffTreatment: + // kMCFlags: use MC flags + // kUA5Cuts: SD events are those that have the MC SD flag and fulfill M^2/s < 0.05; DD from MC flag; Remainder is ND + // kE710Cuts: SD events are those that have the MC SD flag and fulfill 2 < M^2 < 0.05s; DD from MC flag; Remainder is ND + // kALICEHadronLevel: SD events are those that fulfill M^2/s < 0.05; DD from MC flag; Remainder is ND + // + + MCProcessType mcProcessType = GetEventProcessType(header); + + if (diffTreatment == kMCFlags) + return mcProcessType; + + Float_t cms = esd->GetESDRun()->GetBeamEnergy(); + if (esd->GetESDRun()->IsBeamEnergyIsSqrtSHalfGeV()) + cms *= 2; + //Printf("cms = %f", cms); + + if (diffTreatment == kUA5Cuts && mcProcessType == kSD) + { + if (IsHadronLevelSingleDiffractive(stack, cms, 0, 0.05)) + return kSD; + } + else if (diffTreatment == kE710Cuts && mcProcessType == kSD) + { + if (IsHadronLevelSingleDiffractive(stack, cms, 2. / cms / cms, 0.05)) + return kSD; + } + else if (diffTreatment == kALICEHadronLevel) + { + if (IsHadronLevelSingleDiffractive(stack, cms, 0, 0.05)) + return kSD; + } + + if (mcProcessType == kSD) + return kND; + + return mcProcessType; +}