#include <AliHeader.h>
#include <AliStack.h>
-#include <AliLog.h>
#include <AliLog.h>
#include <AliESD.h>
#include <AliESDEvent.h>
#include <AliESDVertex.h>
+#include <AliESDRun.h>
+#include <AliVertexerTracks.h>
+#include <AliMultiplicity.h>
#include <AliGenEventHeader.h>
#include <AliGenPythiaEventHeader.h>
#include <AliGenCocktailEventHeader.h>
+#include <AliGenDPMjetEventHeader.h>
+#include <AliESDVZERO.h>
//____________________________________________________________________
ClassImp(AliPWG0Helper)
-//____________________________________________________________________
-Bool_t AliPWG0Helper::IsEventTriggered(const AliESD* aEsd, Trigger trigger)
-{
- // see function with ULong64_t argument
-
- ULong64_t triggerMask = aEsd->GetTriggerMask();
- return IsEventTriggered(triggerMask, trigger);
-}
+Int_t AliPWG0Helper::fgLastProcessType = -1;
//____________________________________________________________________
-Bool_t AliPWG0Helper::IsEventTriggered(ULong64_t triggerMask, Trigger trigger)
+Bool_t AliPWG0Helper::TestVertex(const AliESDVertex* vertex, AnalysisMode analysisMode, Bool_t debug)
{
- // check if the event was triggered
- //
- // this function needs the branch fTriggerMask
-
+ // Checks if a vertex meets the needed quality criteria
- // definitions from p-p.cfg
- ULong64_t spdFO = (1 << 14);
- ULong64_t v0left = (1 << 11);
- ULong64_t v0right = (1 << 12);
-
- switch (trigger)
+ Float_t requiredZResolution = -1;
+ if (analysisMode & kSPD || analysisMode & kTPCITS || analysisMode & kTPCSPD)
{
- case kMB1:
- {
- if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
- return kTRUE;
- break;
- }
- case kMB2:
- {
- if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
- return kTRUE;
- break;
- }
+ // disable cut on resolution
+ requiredZResolution = 1000;
}
+ else if (analysisMode & kTPC)
+ requiredZResolution = 10.;
- return kFALSE;
-}
-
-//____________________________________________________________________
-Bool_t AliPWG0Helper::IsVertexReconstructed(const AliESD* aEsd)
-{
- // see function with AliESDVertex argument
-
- const AliESDVertex* vtxESD = aEsd->GetVertex();
- return IsVertexReconstructed(vtxESD);
-}
-
-//____________________________________________________________________
-Bool_t AliPWG0Helper::IsVertexReconstructed(const AliESDVertex* vtxESD)
-{
- // checks if the vertex is reasonable
- //
- // this function needs the branches fSPDVertex*
-
- if (!vtxESD)
- return kFALSE;
-
- // the vertex should be reconstructed
- if (strcmp(vtxESD->GetName(), "default")==0)
- return kFALSE;
-
- // check Ncontributors
- if (vtxESD->GetNContributors() <= 0)
- return kFALSE;
-
- Double_t vtx_res[3];
- vtx_res[0] = vtxESD->GetXRes();
- vtx_res[1] = vtxESD->GetYRes();
- vtx_res[2] = vtxESD->GetZRes();
+ // check resolution
+ Double_t zRes = vertex->GetZRes();
- if (vtx_res[2]==0 || vtx_res[2]>0.1)
+ if (zRes > requiredZResolution) {
+ if (debug)
+ 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;
}
// also the quality criteria that are applied)
const AliESDVertex* vertex = 0;
- Float_t requiredZResolution = -1;
- if (analysisMode == kSPD || analysisMode == kTPCITS)
+ if (analysisMode & kSPD)
{
vertex = aEsd->GetPrimaryVertexSPD();
- requiredZResolution = 0.1;
if (debug)
Printf("AliPWG0Helper::GetVertex: Returning SPD vertex");
}
- else if (analysisMode == kTPC)
+ else if (analysisMode & kTPCITS || analysisMode & kTPCSPD)
{
- vertex = aEsd->GetPrimaryVertexTPC();
- requiredZResolution = 0.6;
+ 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);
// check Ncontributors
if (vertex->GetNContributors() <= 0) {
- if (debug)
- Printf("AliPWG0Helper::GetVertex: NContributors() <= 0");
+ if (debug){
+ Printf("AliPWG0Helper::GetVertex: NContributors() <= 0: %d",vertex->GetNContributors());
+ Printf("AliPWG0Helper::GetVertex: NIndices(): %d",vertex->GetNIndices());
+ vertex->Print();
+ }
return 0;
}
// check resolution
Double_t zRes = vertex->GetZRes();
-
- if (zRes == 0 || zRes > requiredZResolution) {
- if (debug)
- Printf("AliPWG0Helper::GetVertex: Resolution too poor %f (required: %f", zRes, requiredZResolution);
+ if (zRes == 0) {
+ Printf("AliPWG0Helper::GetVertex: UNEXPECTED: resolution is 0.");
return 0;
}
if (debug)
- Printf("AliPWG0Helper::GetVertex: Returning valid vertex");
+ {
+ Printf("AliPWG0Helper::GetVertex: Returning valid vertex: %s", vertex->GetTitle());
+ vertex->Print();
+ }
return vertex;
}
}
Int_t pdgCode = TMath::Abs(aParticle->GetPdgCode());
+
// skip quarks and gluon
if (pdgCode <= 10 || pdgCode == 21)
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;
}
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;
}
else if (strlen(axis) == 1)
proj2->SetXTitle(GetAxisTitle(hist2, axis[0]));
- TH1* division = dynamic_cast<TH1*> (proj->Clone(Form("%s_div_%s", proj->GetName(), proj2->GetName())));
+ TH1* division = static_cast<TH1*> (proj->Clone(Form("%s_div_%s", proj->GetName(), proj2->GetName())));
//printf("doing axis: %s, x axis has %d %d bins, min %f %f max %f %f\n", axis, division->GetNbinsX(), proj2->GetNbinsX(), division->GetXaxis()->GetBinLowEdge(1), proj2->GetXaxis()->GetBinLowEdge(1), division->GetXaxis()->GetBinUpEdge(division->GetNbinsX()), proj2->GetXaxis()->GetBinUpEdge(proj2->GetNbinsX()));
//printf("doing axis: %s, y axis has %d %d bins, min %f %f max %f %f\n", axis, division->GetNbinsY(), proj2->GetNbinsY(), division->GetYaxis()->GetBinLowEdge(1), proj2->GetYaxis()->GetBinLowEdge(1), division->GetYaxis()->GetBinUpEdge(division->GetNbinsY()), proj2->GetYaxis()->GetBinUpEdge(proj2->GetNbinsY()));
division->Divide(proj, proj2, 1, 1, "B");
return 0;
}
-//____________________________________________________________________
-Int_t AliPWG0Helper::GetPythiaEventProcessType(AliHeader* aHeader, Bool_t adebug) {
+
+AliPWG0Helper::MCProcessType AliPWG0Helper::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
+
+ AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(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<AliGenDPMjetEventHeader*>(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<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
+ if (pythiaGenHeader) {
+ return GetPythiaEventProcessType(pythiaGenHeader,adebug);
+ }
- if (!pythiaGenHeader) {
+ AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader->GenEventHeader());
+ if (dpmJetGenHeader) {
+ return GetDPMjetEventProcessType(dpmJetGenHeader,adebug);
+ }
+
- AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(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<AliGenCocktailEventHeader*>(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; i<headerList->GetEntries(); i++) {
- pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
- if (pythiaGenHeader)
- break;
+ for (Int_t i=0; i<headerList->GetEntries(); i++) {
+
+ pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(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<AliGenDPMjetEventHeader*>(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)
return label;
}
-void AliPWG0Helper::SetBranchStatusRecursive(TTree* tree, char *bname, Bool_t status, Bool_t debug)
-{
- // Function to switch on/off all data members of a top level branch
- // this is needed for branches without a trailing dot ".", for those
- // the root functionality with regular expressions does not work.
- // Usage e.g.
- // chain->SetBranchStatus("*", 0);
- // SetBranchStatusRecursive(chain,"SPDVertex",1);
- // You need to give the full name of the top level branch zou want to access
- //==========================================================================
- // Author Christian.Klein-Boesing@cern.ch
-
- if (!tree)
- return;
-
- TBranch *br = tree->GetBranch(bname);
- if(!br) {
- Printf("AliPWG0Helper::SetBranchStatusRecursive: Branch %s not found", bname);
- }
-
- TObjArray *leaves = tree->GetListOfLeaves();
- Int_t nleaves = leaves->GetEntries();
- TLeaf *leaf = 0;
- TBranch *branch = 0;
- TBranch *mother = 0;
- for (Int_t i=0;i<nleaves;i++) {
- // the matched entry e.g. SPDVertex is its own Mother
- leaf = (TLeaf*)leaves->UncheckedAt(i);
- branch = (TBranch*)leaf->GetBranch();
- mother = branch->GetMother();
- if (mother==br){
- if (debug)
- Printf(">>>> AliPWG0Helper::SetBranchStatusRecursive: Setting status %s %s to %d", mother->GetName(), leaf->GetName(), status);
- if (status) branch->ResetBit(kDoNotProcess);
- else branch->SetBit(kDoNotProcess);
- }
- }
-}
-
//____________________________________________________________________
void AliPWG0Helper::NormalizeToBinWidth(TH1* hist)
{
}
//____________________________________________________________________
-void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, Trigger trigger)
+void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, AliTriggerAnalysis::Trigger trigger, AliPWG0Helper::DiffTreatment diffTreatment)
{
//
// Prints the given configuration
//
- TString str(">>>> Running with ");
+ 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);
+}
- switch (analysisMode)
+//____________________________________________________________________
+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++)
{
- case kSPD : str += "SPD-only"; break;
- case kTPC : str += "TPC-only"; break;
- case kTPCITS : str += "Global tracking"; break;
+ 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;
- str += " and trigger ";
+ Int_t pdg1 = part1->GetPdgCode();
+ Int_t pdg2 = part2->GetPdgCode();
- switch (trigger)
+ 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)
{
- case kMB1 : str += "MB1"; break;
- case kMB2 : str += "MB2"; break;
+ if (y1 > y2)
+ arm = 0;
+ else
+ arm = 1;
}
+ else if (pdg1 == 2212)
+ arm = 0;
+ else if (pdg2 == 2212)
+ arm = 1;
- str += " <<<<";
+ Double_t M02s = 1. - 2 * part1->Energy() / cms;
+ Double_t M12s = 1. - 2 * part2->Energy() / cms;
- Printf("%s", str.Data());
+ 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;
+
+ if (!esd)
+ {
+ Printf("ERROR: AliPWG0Helper::GetEventProcessType: diffTreatment != kMCFlags and esd == 0");
+ return kInvalidProcess;
+ }
+
+ 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;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff