#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)
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
-
- // definitions from p-p.cfg
- ULong64_t spdFO = (1 << 14);
- ULong64_t v0left = (1 << 11);
- ULong64_t v0right = (1 << 12);
+ // 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 & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
- return kTRUE;
- break;
- }
- case kMB2:
- {
- if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
- return kTRUE;
- break;
- }
- case kSPDFASTOR:
- {
- if (triggerMask & spdFO)
- 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;
+ 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
//
// 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)
{
- if(bRedoTPC){
- Double_t kBz = aEsd->GetMagneticField();
- AliVertexerTracks vertexer(kBz);
- vertexer.SetTPCMode();
- AliESDVertex *vTPC = vertexer.FindPrimaryVertex(aEsd);
- aEsd->SetPrimaryVertexTPC(vTPC);
- for (Int_t i=0; i<aEsd->GetNumberOfTracks(); i++) {
- AliESDtrack *t = aEsd->GetTrack(i);
- t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
- }
- delete vTPC;
+ 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();
- requiredZResolution = 10.;
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);
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;
}
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");
}
- 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;
}
//____________________________________________________________________
-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";
- switch (analysisMode)
+ if (analysisMode & kFieldOn)
{
- 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 += " (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 += "< <<<<";
- str += " and trigger ";
+ Printf("%s", str.Data());
+}
- switch (trigger)
+//____________________________________________________________________
+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++)
{
- case kMB1 : str += "MB1"; break;
- case kMB2 : str += "MB2"; break;
- case kSPDFASTOR : str += "SPD FASTOR"; 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 += " <<<<";
+ Int_t pdg1 = part1->GetPdgCode();
+ Int_t pdg2 = part2->GetPdgCode();
- Printf("%s", str.Data());
+ 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;
+
+ 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