X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWG0%2FAliPWG0Helper.cxx;h=ef7f85910b28b6a94e526a6e0591a8298df56017;hb=91dc3cdd71328fd460f2b5c3a028d9ea7f4cc33f;hp=99c24a633ca554cd2cb247e7b1105b4d9961740e;hpb=116083b1e4e27625a7a0144360efd52d32986522;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWG0/AliPWG0Helper.cxx b/PWG0/AliPWG0Helper.cxx index 99c24a633ca..ef7f85910b2 100644 --- a/PWG0/AliPWG0Helper.cxx +++ b/PWG0/AliPWG0Helper.cxx @@ -5,8 +5,12 @@ #include #include #include +#include #include #include +#include +#include +#include #include #include @@ -14,22 +18,94 @@ #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::IsEventTriggered(const AliESDEvent* aEsd, Trigger trigger) { - // see function with ULong64_t argument + // checks if an event has been triggered + // this function implements the "offline" methods that use the ESD, other trigger requests are passed to the function prototype with ULong_t + + Int_t firedChips = 0; + Bool_t v0A = kFALSE; + Bool_t v0C = kFALSE; - ULong64_t triggerMask = aEsd->GetTriggerMask(); - return IsEventTriggered(triggerMask, trigger); + // offline triggers have to be dealt with here, because we need the esd pointer + if (trigger == kOfflineFASTOR || trigger == kOfflineMB1 || trigger == kOfflineMB2 || trigger == kOfflineMB3) + { + const AliMultiplicity* mult = aEsd->GetMultiplicity(); + if (!mult) + { + Printf("AliPWG0Helper::IsEventTriggered: ERROR: AliMultiplicity not available"); + return kFALSE; + } + firedChips = mult->GetNumberOfFiredChips(0) + mult->GetNumberOfFiredChips(1); + } + if (trigger == kOfflineMB1 || trigger == kOfflineMB2 || trigger == kOfflineMB3) + { + AliESDVZERO* v0Data = aEsd->GetVZEROData(); + if (!v0Data) + { + Printf("AliPWG0Helper::IsEventTriggered: ERROR: AliESDVZERO not available"); + return kFALSE; + } + for (Int_t i=0; i<32; i++) + { + if (v0Data->BBTriggerV0A(i)) + v0A = kTRUE; + if (v0Data->BBTriggerV0C(i)) + v0C = kTRUE; + } + } + + switch (trigger) + { + case kOfflineFASTOR: + { + if (firedChips > 0) + return kTRUE; + break; + } + case kOfflineMB1: + { + if ((firedChips > 0) || v0A || v0C) + return kTRUE; + break; + } + case kOfflineMB2: + { + if ((firedChips > 0) && (v0A || v0C)) + return kTRUE; + break; + } + case kOfflineMB3: + { + if ((firedChips > 0) && v0A && v0C) + return kTRUE; + break; + } + default: + { + return IsEventTriggered(aEsd->GetTriggerMask(), trigger); + break; + } + } + + return kFALSE; } //____________________________________________________________________ @@ -38,65 +114,140 @@ 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 + + // definitions from p-p.cfg + ULong64_t spdFO = (1 << 14); + ULong64_t v0left = (1 << 10); + ULong64_t v0right = (1 << 11); switch (trigger) { case kMB1: { - if (triggerMask&32 || ((triggerMask&1) || (triggerMask&2))) + if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right))) return kTRUE; break; } case kMB2: { - if (triggerMask&32 && ((triggerMask&1) || (triggerMask&2))) + if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right))) + return kTRUE; + break; + } + case kMB3: + { + if (triggerMask & spdFO && (triggerMask & v0left) && (triggerMask & v0right)) + return kTRUE; + break; + } + case kSPDFASTOR: + { + 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::IsVertexReconstructed(const AliESD* aEsd) +Bool_t AliPWG0Helper::TestVertex(const AliESDVertex* vertex, AnalysisMode analysisMode, Bool_t debug) { - // see function with AliESDVertex argument + // Checks if a vertex meets the needed quality criteria + + Float_t requiredZResolution = -1; + if (analysisMode == kSPD || analysisMode == kTPCITS) + { + requiredZResolution = 0.1; + } + else if (analysisMode == kTPC) + requiredZResolution = 10.; + + // check resolution + Double_t zRes = vertex->GetZRes(); - const AliESDVertex* vtxESD = aEsd->GetVertex(); - return IsVertexReconstructed(vtxESD); + if (zRes > requiredZResolution) { + if (debug) + Printf("AliPWG0Helper::TestVertex: Resolution too poor %f (required: %f", zRes, requiredZResolution); + return kFALSE; + } + + return kTRUE; } //____________________________________________________________________ -Bool_t AliPWG0Helper::IsVertexReconstructed(const AliESDVertex* vtxESD) +const AliESDVertex* AliPWG0Helper::GetVertex(AliESDEvent* aEsd, AnalysisMode analysisMode, Bool_t debug, Bool_t bRedoTPC) { - // 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 || analysisMode == kTPCITS) + { + vertex = aEsd->GetPrimaryVertexSPD(); + if (debug) + Printf("AliPWG0Helper::GetVertex: Returning SPD vertex"); + } + else if (analysisMode == kTPC) + { + if(bRedoTPC){ + 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; + } - // the vertex should be reconstructed - if (strcmp(vtxESD->GetName(), "default")==0) - return kFALSE; + vertex = aEsd->GetPrimaryVertexTPC(); + if (debug) + Printf("AliPWG0Helper::GetVertex: Returning vertex from tracks"); + } + else + Printf("AliPWG0Helper::GetVertex: ERROR: Invalid second argument %d", analysisMode); - Double_t vtx_res[3]; - vtx_res[0] = vtxESD->GetXRes(); - vtx_res[1] = vtxESD->GetYRes(); - vtx_res[2] = vtxESD->GetZRes(); + if (!vertex) { + if (debug) + Printf("AliPWG0Helper::GetVertex: No vertex found in ESD"); + return 0; + } - if (vtx_res[2]==0 || vtx_res[2]>0.1) - return kFALSE; + // 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; + } - // check Ncontributors, if <0 it means error *gna* + // check resolution + Double_t zRes = vertex->GetZRes(); + if (zRes == 0) { + Printf("AliPWG0Helper::GetVertex: UNEXPECTED: resolution is 0."); + return 0; + } - return kTRUE; + if (debug) + { + Printf("AliPWG0Helper::GetVertex: Returning valid vertex: %s", vertex->GetTitle()); + vertex->Print(); + } + + return vertex; } //____________________________________________________________________ @@ -108,6 +259,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 +271,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 +281,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 +299,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 +423,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){ // this is explicitly inelastic + globalType = kND; + } + else if(dpmJetType==5||dpmJetType==6){ + globalType = kSD; + } + else if(dpmJetType==7||dpmJetType==4){// DD and double pomeron + 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++) { - pythiaGenHeader = dynamic_cast(headerList->At(i)); - if (pythiaGenHeader) - break; + for (Int_t i=0; iGetEntries(); i++) { + + 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 +593,70 @@ 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, Trigger trigger) +{ + // + // Prints the given configuration + // + + TString str(">>>> Running with "); + + switch (analysisMode) + { + case kInvalid: str += "invalid setting"; break; + case kSPD : str += "SPD-only"; break; + case kTPC : str += "TPC-only"; break; + case kTPCITS : str += "Global tracking"; break; + } + + str += " and trigger "; + + switch (trigger) + { + case kMB1 : str += "MB1"; break; + case kMB2 : str += "MB2"; break; + case kMB3 : str += "MB3"; break; + case kSPDFASTOR : str += "SPD FASTOR"; break; + case kOfflineMB1 : str += "Offline MB1"; break; + case kOfflineMB2 : str += "Offline MB2"; break; + case kOfflineMB3 : str += "Offline MB3"; break; + case kOfflineFASTOR : str += "Offline SPD FASTOR"; break; + } + + str += " <<<<"; + + Printf("%s", str.Data()); +} +