X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWG0%2FAliPWG0Helper.cxx;h=7ca5f738a75b1787b76cb32c2b7236df12485cef;hb=b3659bade5509cd8023cd9ca6d0e1480552fd7ad;hp=99c24a633ca554cd2cb247e7b1105b4d9961740e;hpb=116083b1e4e27625a7a0144360efd52d32986522;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWG0/AliPWG0Helper.cxx b/PWG0/AliPWG0Helper.cxx index 99c24a633ca..7ca5f738a75 100644 --- a/PWG0/AliPWG0Helper.cxx +++ b/PWG0/AliPWG0Helper.cxx @@ -5,98 +5,137 @@ #include #include #include +#include #include #include +#include +#include +#include #include #include -#include #include #include +#include #include +#include +#include +#include #include #include #include +#include +#include //____________________________________________________________________ ClassImp(AliPWG0Helper) +Int_t AliPWG0Helper::fgLastProcessType = -1; + //____________________________________________________________________ -Bool_t AliPWG0Helper::IsEventTriggered(const AliESD* aEsd, Trigger trigger) +Bool_t AliPWG0Helper::TestVertex(const AliESDVertex* vertex, AnalysisMode analysisMode, Bool_t debug) { - // see function with ULong64_t argument + // Checks if a vertex meets the needed quality criteria - ULong64_t triggerMask = aEsd->GetTriggerMask(); - return IsEventTriggered(triggerMask, trigger); -} + Float_t requiredZResolution = -1; + if (analysisMode & kSPD || analysisMode & kTPCITS || analysisMode & kTPCSPD) + { + // disable cut on resolution + requiredZResolution = 1000; + } + else if (analysisMode & kTPC) + requiredZResolution = 10.; -//____________________________________________________________________ -Bool_t AliPWG0Helper::IsEventTriggered(ULong64_t triggerMask, Trigger trigger) -{ - // check if the event was triggered - // - // this function needs the branch fTriggerMask - // - // MB should be - // ITS_SPD_GFO_L0 : 32 - // VZERO_OR_LEFT : 1 - // VZERO_OR_RIGHT : 2 + // check resolution + Double_t zRes = vertex->GetZRes(); - switch (trigger) + if (zRes > requiredZResolution) { + if (debug) + Printf("AliPWG0Helper::TestVertex: Resolution too poor %f (required: %f", zRes, requiredZResolution); + return kFALSE; + } + + if (vertex->IsFromVertexerZ()) { - case kMB1: + if (vertex->GetDispersion() > 0.02) { - if (triggerMask&32 || ((triggerMask&1) || (triggerMask&2))) - return kTRUE; - break; - } - case kMB2: - { - if (triggerMask&32 && ((triggerMask&1) || (triggerMask&2))) - return kTRUE; - break; + if (debug) + Printf("AliPWG0Helper::TestVertex: Delta Phi too large in Vertexer Z: %f (required: %f", vertex->GetDispersion(), 0.02); + return kFALSE; } } - return kFALSE; -} - -//____________________________________________________________________ -Bool_t AliPWG0Helper::IsVertexReconstructed(const AliESD* aEsd) -{ - // see function with AliESDVertex argument - - const AliESDVertex* vtxESD = aEsd->GetVertex(); - return IsVertexReconstructed(vtxESD); + return kTRUE; } //____________________________________________________________________ -Bool_t AliPWG0Helper::IsVertexReconstructed(const AliESDVertex* vtxESD) +const AliESDVertex* AliPWG0Helper::GetVertex(const AliESDEvent* aEsd, AnalysisMode analysisMode, Bool_t debug) { - // checks if the vertex is reasonable + // Get the vertex from the ESD and returns it if the vertex is valid // - // this function needs the branches fSPDVertex* + // Second argument decides which vertex is used (this selects + // also the quality criteria that are applied) - if (!vtxESD) - return kFALSE; + const AliESDVertex* vertex = 0; + if (analysisMode & kSPD) + { + vertex = aEsd->GetPrimaryVertexSPD(); + if (debug) + Printf("AliPWG0Helper::GetVertex: Returning SPD vertex"); + } + else if (analysisMode & kTPCITS || analysisMode & kTPCSPD) + { + vertex = aEsd->GetPrimaryVertexTracks(); + if (debug) + Printf("AliPWG0Helper::GetVertex: Returning vertex from tracks"); + if (!vertex || vertex->GetNContributors() <= 0) + { + if (debug) + 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 TPC-only tracks"); + } + else + Printf("AliPWG0Helper::GetVertex: ERROR: Invalid second argument %d", analysisMode); - // the vertex should be reconstructed - if (strcmp(vtxESD->GetName(), "default")==0) - return kFALSE; + if (!vertex) { + if (debug) + Printf("AliPWG0Helper::GetVertex: No vertex found in ESD"); + return 0; + } - Double_t vtx_res[3]; - vtx_res[0] = vtxESD->GetXRes(); - vtx_res[1] = vtxESD->GetYRes(); - vtx_res[2] = vtxESD->GetZRes(); + // check Ncontributors + if (vertex->GetNContributors() <= 0) { + if (debug){ + Printf("AliPWG0Helper::GetVertex: NContributors() <= 0: %d",vertex->GetNContributors()); + Printf("AliPWG0Helper::GetVertex: NIndices(): %d",vertex->GetNIndices()); + vertex->Print(); + } + return 0; + } - if (vtx_res[2]==0 || vtx_res[2]>0.1) - return kFALSE; + // check resolution + Double_t zRes = vertex->GetZRes(); + if (zRes == 0) { + Printf("AliPWG0Helper::GetVertex: UNEXPECTED: resolution is 0."); + return 0; + } - // check Ncontributors, if <0 it means error *gna* + if (debug) + { + Printf("AliPWG0Helper::GetVertex: Returning valid vertex: %s", vertex->GetTitle()); + vertex->Print(); + } - return kTRUE; + return vertex; } //____________________________________________________________________ @@ -108,6 +147,8 @@ Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimari // // This function or a equivalent should be available in some common place of AliRoot // + // WARNING: Call this function only for particles that are among the particles from the event generator! + // --> stack->Particle(id) with id < stack->GetNprimary() // if the particle has a daughter primary, we do not want to count it if (aParticle->GetFirstDaughter() != -1 && aParticle->GetFirstDaughter() < aTotalPrimaries) @@ -118,6 +159,7 @@ Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimari } Int_t pdgCode = TMath::Abs(aParticle->GetPdgCode()); + // skip quarks and gluon if (pdgCode <= 10 || pdgCode == 21) @@ -127,10 +169,17 @@ Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimari return kFALSE; } + Int_t status = aParticle->GetStatusCode(); + // skip non final state particles.. + if(status!=1){ + if (adebug) + printf("Dropping particle because it is not a final state particle.\n"); + return kFALSE; + } + if (strcmp(aParticle->GetName(),"XXX") == 0) { - if (adebug) - printf("WARNING: There is a particle named XXX.\n"); + Printf("WARNING: There is a particle named XXX (pdg code %d).", pdgCode); return kFALSE; } @@ -138,8 +187,7 @@ Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimari if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0) { - if (adebug) - printf("WARNING: There is a particle with an unknown particle class (pdg code %d).\n", pdgCode); + Printf("WARNING: There is a particle with an unknown particle class (pdg code %d).", pdgCode); return kFALSE; } @@ -263,46 +311,125 @@ const char* AliPWG0Helper::GetAxisTitle(TH3* hist, const char axis) return 0; } -//____________________________________________________________________ -Int_t AliPWG0Helper::GetPythiaEventProcessType(AliHeader* aHeader, Bool_t adebug) { + +AliPWG0Helper::MCProcessType AliPWG0Helper::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) { + + AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast(aHeader); + + if (!pythiaGenHeader) { + printf("AliPWG0Helper::GetProcessType : Unknown gen Header type). \n"); + return kInvalidProcess; + } + + + Int_t pythiaType = pythiaGenHeader->ProcessType(); + fgLastProcessType = pythiaType; + MCProcessType globalType = kInvalidProcess; + + + if (adebug) { + printf("AliPWG0Helper::GetProcessType : Pythia process type found: %d \n",pythiaType); + } + + + if(pythiaType==92||pythiaType==93){ + globalType = kSD; + } + else if(pythiaType==94){ + globalType = kDD; + } + //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB?? + else { + globalType = kND; + } + return globalType; +} + + +AliPWG0Helper::MCProcessType AliPWG0Helper::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) { // // get the process type of the event. // // can only read pythia headers, either directly or from cocktalil header + AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast(aHeader); + + if (!dpmJetGenHeader) { + printf("AliPWG0Helper::GetDPMjetProcessType : Unknown header type (not DPMjet or). \n"); + return kInvalidProcess; + } + + Int_t dpmJetType = dpmJetGenHeader->ProcessType(); + fgLastProcessType = dpmJetType; + MCProcessType globalType = kInvalidProcess; + + + if (adebug) { + printf("AliPWG0Helper::GetDPMJetProcessType : DPMJet process type found: %d \n",dpmJetType); + } + + + if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction + globalType = kND; + } + else if (dpmJetType==5 || dpmJetType==6) { + globalType = kSD; + } + else if (dpmJetType==7) { + globalType = kDD; + } + return globalType; +} + + +AliPWG0Helper::MCProcessType AliPWG0Helper::GetEventProcessType(AliHeader* aHeader, Bool_t adebug) { + // + // get the process type of the event. + // + + + // Check for simple headers first + AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast(aHeader->GenEventHeader()); + if (pythiaGenHeader) { + return GetPythiaEventProcessType(pythiaGenHeader,adebug); + } - if (!pythiaGenHeader) { + AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast(aHeader->GenEventHeader()); + if (dpmJetGenHeader) { + return GetDPMjetEventProcessType(dpmJetGenHeader,adebug); + } + - AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast(aHeader->GenEventHeader()); - if (!genCocktailHeader) { - printf("AliPWG0Helper::GetProcessType : Unknown header type (not Pythia or Cocktail). \n"); - return -1; - } + // check for cocktail - TList* headerList = genCocktailHeader->GetHeaders(); - if (!headerList) { - return -1; - } + AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast(aHeader->GenEventHeader()); + if (!genCocktailHeader) { + printf("AliPWG0Helper::GetProcessType : Unknown header type (not Pythia or Cocktail). \n"); + return kInvalidProcess; + } + + TList* headerList = genCocktailHeader->GetHeaders(); + if (!headerList) { + return kInvalidProcess; + } + + for (Int_t i=0; iGetEntries(); i++) { - for (Int_t i=0; iGetEntries(); i++) { - pythiaGenHeader = dynamic_cast(headerList->At(i)); - if (pythiaGenHeader) - break; + pythiaGenHeader = dynamic_cast(headerList->At(i)); + if (pythiaGenHeader) { + return GetPythiaEventProcessType(pythiaGenHeader,adebug); } - if (!pythiaGenHeader) { - printf("AliPWG0Helper::GetProcessType : Could not find Pythia header. \n"); - return -1; + dpmJetGenHeader = dynamic_cast(headerList->At(i)); + if (dpmJetGenHeader) { + return GetDPMjetEventProcessType(dpmJetGenHeader,adebug); } } + return kInvalidProcess; +} - if (adebug) { - printf("AliPWG0Helper::GetProcessType : Pythia process type found: %d \n",pythiaGenHeader->ProcessType()); - } - return pythiaGenHeader->ProcessType(); -} //____________________________________________________________________ TParticle* AliPWG0Helper::FindPrimaryMother(AliStack* stack, Int_t label) @@ -354,3 +481,253 @@ Int_t AliPWG0Helper::FindPrimaryMotherLabel(AliStack* stack, Int_t label) return label; } + +//____________________________________________________________________ +void AliPWG0Helper::NormalizeToBinWidth(TH1* hist) +{ + // + // normalizes a 1-d histogram to its bin width + // + + for (Int_t i=1; i<=hist->GetNbinsX(); ++i) + { + hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i)); + hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i)); + } +} + +//____________________________________________________________________ +void AliPWG0Helper::NormalizeToBinWidth(TH2* hist) +{ + // + // normalizes a 2-d histogram to its bin width (x width * y width) + // + + for (Int_t i=1; i<=hist->GetNbinsX(); ++i) + for (Int_t j=1; j<=hist->GetNbinsY(); ++j) + { + Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j); + hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor); + hist->SetBinError(i, j, hist->GetBinError(i, j) / factor); + } +} + +//____________________________________________________________________ +void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, AliTriggerAnalysis::Trigger trigger, AliPWG0Helper::DiffTreatment diffTreatment) +{ + // + // Prints the given configuration + // + + TString str(">>>> Running with >"); + + if (analysisMode & kSPD) + str += "SPD-only"; + + if (analysisMode & kSPDOnlyL0) + str += " (only L0 clusters)"; + + if (analysisMode & kTPC) + str += "TPC-only"; + + if (analysisMode & kTPCITS) + str += "Global tracking"; + + if (analysisMode & kTPCSPD) + str += "Tracks and tracklets"; + + if (analysisMode & kFieldOn) + { + str += " (with field)"; + } + else + str += " (WITHOUT field)"; + + str += "< and 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; +}