// mult axis for event histogram
Int_t nBinsN = 22;
Float_t binLimitsN[] = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 12.5, 14.5, 16.5, 18.5, 20.5, 25.5, 30.5, 40.5, 50.5, 100.5, 300.5};
+ //Int_t nBinsN = 8;
+ //Float_t binLimitsN[] = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 300.5};
//Int_t nBinsN = 52;
//Float_t binLimitsN[] = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 100.5, 300.5};
//Float_t binLimitsVtx[] = {-20,-15,-10,-6,-3,0,3,6,10,15,20};
//Float_t binLimitsVtx[] = {-20,-19,-18,-17,-16,-15,-14,-13,-12,-11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20};
- Float_t binLimitsVtx[] = {-30,-25,-20,-15,-10,-8,-6,-4,-2,0,2,4,6,8,10,15,20,25,30};
+ //Float_t binLimitsVtx[] = {-30,-25,-20,-15,-10,-8,-6,-4,-2,0,2,4,6,8,10,15,20,25,30};
+ Float_t binLimitsVtx[] = {-30,-25,-20,-15,-10,-7,-5.5,-4,-3,-2,-1,0,1,2,3,4,5.5,7,10,15,20,25,30};
/*Float_t binLimitsEta[] = {-3.0,-2.8,-2.6,-2.4,-2.2,
-2.0,-1.8,-1.6,-1.4,-1.2,
-1.0,-0.8,-0.6,-0.4,-0.2, 0.0,
0,0.25,0.5,0.75,
1.0,1.25,1.5,1.75,
2.0,2.25,2.5,2.75,3.0}; */
- Float_t binLimitsEta[] = {-2.8,-2.4,-2,-1.6,-1.2,-0.8,-0.4,0,0.4,0.8,1.2,1.6,2.0,2.4,2.8};
+ //Float_t binLimitsEta[] = {-2.8,-2.4,-2,-1.6,-1.2,-0.8,-0.4,0,0.4,0.8,1.2,1.6,2.0,2.4,2.8};
//Float_t binLimitsEta[] = {-2.8,-2.4,-2,-1.6,-1.2,-0.8,-0.5,0,0.5,0.8,1.2,1.6,2.0,2.4,2.8};
//Float_t binLimitsEta[] = {-3,-2.6,-2.2,-1.8,-1.4,-1,-0.6,-0.2,0.2,0.6,1,1.4,1.8,2.2,2.6,3.0};
-// Float_t binLimitsEta[] = {-2,-1.9,-1.8,-1.7,-1.6,-1.5,-1.4,-1.3,-1.2,-1.1,-1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0};
+ Float_t binLimitsEta[] = {-2,-1.9,-1.8,-1.7,-1.6,-1.5,-1.4,-1.3,-1.2,-1.1,-1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0};
// Float_t binLimitsEta[] = { -7.0, -6.9, -6.8, -6.7, -6.6, -6.5, -6.4, -6.3, -6.2, -6.1, -6.0, -5.9, -5.8, -5.7, -5.6, -5.5, -5.4, -5.3, -5.2, -5.1, -5.0, -4.9, -4.8, -4.7, -4.6, -4.5, -4.4, -4.3, -4.2, -4.1, -4.0, -3.9, -3.8, -3.7, -3.6, -3.5, -3.4, -3.3, -3.2, -3.1, -3.0, -2.9, -2.8, -2.7, -2.6, -2.5, -2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, -0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5.0, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6.0, 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7.0 };
- fEventCorr = new AliCorrectionMatrix2D("EventCorrection", Form("%s EventCorrection", fTitle.Data()), 18, binLimitsVtx, nBinsN, binLimitsN);
+ fEventCorr = new AliCorrectionMatrix2D("EventCorrection", Form("%s EventCorrection", fTitle.Data()), 22, binLimitsVtx, nBinsN, binLimitsN);
if (analysisMode & AliPWG0Helper::kSPD)
{
- TH3F* dummyBinning = new TH3F("dummyBinning","dummyBinning",18, binLimitsVtx, 14, binLimitsEta, nBinsN2, binLimitsN2);
+ TH3F* dummyBinning = new TH3F("dummyBinning","dummyBinning",22, binLimitsVtx, 40, binLimitsEta, nBinsN2, binLimitsN2);
fTrackCorr = new AliCorrectionMatrix3D("TrackCorrection", Form("%s TrackCorrection", fTitle.Data()), dummyBinning);
fTrackCorr->SetAxisTitles("vtx-z (cm)", "#eta", title3);
delete dummyBinning;
}
else
{
- TH3F* dummyBinning = new TH3F("dummyBinning","dummyBinning",18, binLimitsVtx, 14, binLimitsEta , nBinsPt, binLimitsPt);
+ TH3F* dummyBinning = new TH3F("dummyBinning","dummyBinning",22, binLimitsVtx, 40, binLimitsEta , nBinsPt, binLimitsPt);
fTrackCorr = new AliCorrectionMatrix3D("TrackCorrection", Form("%s TrackCorrection", fTitle.Data()), dummyBinning);
fTrackCorr->SetAxisTitles("vtx-z (cm)", "#eta", "p_{T} (GeV/c)");
delete dummyBinning;
Float_t eventsM = measuredEvents->Integral(measuredEvents->GetXaxis()->FindBin(-zRange), measuredEvents->GetXaxis()->FindBin(zRange), 1, measuredEvents->GetNbinsY());
Float_t eventsG = generatedEvents->Integral(generatedEvents->GetXaxis()->FindBin(-zRange), generatedEvents->GetXaxis()->FindBin(zRange), 1, generatedEvents->GetNbinsY());
- Printf("tracks measured: %.1f tracks generated: %.1f, events measured: %.1f, events generated %.1f", tracksM, tracksG, eventsM, eventsG);
+ Float_t eventsM1 = measuredEvents->Integral(measuredEvents->GetXaxis()->FindBin(-zRange), measuredEvents->GetXaxis()->FindBin(zRange), 2, measuredEvents->GetNbinsY());
+ Float_t eventsG1 = generatedEvents->Integral(generatedEvents->GetXaxis()->FindBin(-zRange), generatedEvents->GetXaxis()->FindBin(zRange), 2, generatedEvents->GetNbinsY());
+
+ Printf("tracks measured: %.1f tracks generated: %.1f, events measured: %.1f, events generated %.1f, events (N>0) measured: %.1f, events (N>0) generated %.1f", tracksM, tracksG, eventsM, eventsG, eventsM1, eventsG1);
if (tracksM > 0)
Printf("Effective average correction factor for TRACKS: %.3f", tracksG / tracksM);
#include <AliESD.h>
#include <AliESDEvent.h>
#include <AliESDVertex.h>
+#include <AliESDRun.h>
#include <AliVertexerTracks.h>
#include <AliMultiplicity.h>
}
//____________________________________________________________________
-const AliESDVertex* AliPWG0Helper::GetVertex(AliESDEvent* aEsd, AnalysisMode analysisMode, Bool_t debug)
+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;
- if (analysisMode & kSPD || analysisMode & kTPCITS)
+ if (analysisMode & kSPD)
{
vertex = aEsd->GetPrimaryVertexSPD();
if (debug)
Printf("AliPWG0Helper::GetVertex: Returning SPD vertex");
}
+ else if (analysisMode & kTPCITS)
+ {
+ 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 tracks");
+ Printf("AliPWG0Helper::GetVertex: Returning vertex from TPC-only tracks");
}
else
Printf("AliPWG0Helper::GetVertex: ERROR: Invalid second argument %d", analysisMode);
}
- 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, AliTriggerAnalysis::Trigger trigger)
+void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, AliTriggerAnalysis::Trigger trigger, AliPWG0Helper::DiffTreatment diffTreatment)
{
//
// Prints the given configuration
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;
+}
class AliStack;
class TTree;
class AliOfflineTrigger;
+class AliMultiplicity;
class AliPWG0Helper : public TObject
{
public:
enum AnalysisMode { kInvalid = -1, kSPD = 0x1, kTPC = 0x2, kTPCITS = 0x4, kFieldOn = 0x8, kSPDOnlyL0 = 0x10 };
// in case we want to use bitmaps...
- enum MCProcessType { kInvalidProcess = -1, kND = 0x1, kDD = 0x2, kSD = 0x4 };
+ enum MCProcessType { kInvalidProcess = -1, kND = 0x1, kDD = 0x2, kSD = 0x4, kOnePart = 0x8 };
+ enum DiffTreatment { kMCFlags = 0, kUA5Cuts = 1, kE710Cuts, kALICEHadronLevel };
- static const AliESDVertex* GetVertex(AliESDEvent* aEsd, AnalysisMode analysisMethod, Bool_t debug = kFALSE);
+ static const AliESDVertex* GetVertex(const AliESDEvent* aEsd, AnalysisMode analysisMethod, Bool_t debug = kFALSE);
static Bool_t TestVertex(const AliESDVertex* vertex, AnalysisMode analysisMode, Bool_t debug = kFALSE);
static Bool_t IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimaries, Bool_t adebug = kFALSE);
- static AliPWG0Helper::MCProcessType GetEventProcessType(AliHeader* aHeader, Bool_t adebug = kFALSE);
- static AliPWG0Helper::MCProcessType GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug = kFALSE);
- static AliPWG0Helper::MCProcessType GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug = kFALSE);
+ static MCProcessType GetEventProcessType(AliESDEvent* esd, AliHeader* header, AliStack* stack, DiffTreatment diffTreatment);
+ static MCProcessType GetEventProcessType(AliHeader* aHeader, Bool_t adebug = kFALSE);
+ static MCProcessType GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug = kFALSE);
+ static MCProcessType GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug = kFALSE);
static Int_t GetLastProcessType() { return fgLastProcessType; }
static TParticle* FindPrimaryMother(AliStack* stack, Int_t label);
static void NormalizeToBinWidth(TH1* hist);
static void NormalizeToBinWidth(TH2* hist);
- static void PrintConf(AnalysisMode analysisMode, AliTriggerAnalysis::Trigger trigger);
+ static void PrintConf(AnalysisMode analysisMode, AliTriggerAnalysis::Trigger trigger, DiffTreatment diffTreatment);
+
+ static Double_t Rapidity(Double_t pt, Double_t pz, Double_t m);
+ static Bool_t IsHadronLevelSingleDiffractive(AliStack* stack, Float_t cms, Float_t xiMin, Float_t xiMax);
protected:
static Int_t fgLastProcessType; // stores the raw value of the last process type extracted
fVertexRecoCorrection(0),
fTriggerBiasCorrectionMBToINEL(0),
fTriggerBiasCorrectionMBToNSD(0),
- fTriggerBiasCorrectionMBToND(0)
+ fTriggerBiasCorrectionMBToND(0),
+ fTriggerBiasCorrectionMBToOnePart(0)
{
// default constructor
}
fVertexRecoCorrection(0),
fTriggerBiasCorrectionMBToINEL(0),
fTriggerBiasCorrectionMBToNSD(0),
- fTriggerBiasCorrectionMBToND(0)
+ fTriggerBiasCorrectionMBToND(0),
+ fTriggerBiasCorrectionMBToOnePart(0)
{
//
// constructor
fTriggerBiasCorrectionMBToINEL = new AliCorrection("TriggerBias_MBToINEL", "TriggerBias_MBToINEL", analysis);
fTriggerBiasCorrectionMBToNSD = new AliCorrection("TriggerBias_MBToNSD", "TriggerBias_MBToNSD", analysis);
fTriggerBiasCorrectionMBToND = new AliCorrection("TriggerBias_MBToND", "TriggerBias_MBToND", analysis);
+ fTriggerBiasCorrectionMBToOnePart = new AliCorrection("TriggerBias_MBToOnePart", "TriggerBias_MBToOnePart", analysis);
}
//____________________________________________________________________
delete fTriggerBiasCorrectionMBToND;
fTriggerBiasCorrectionMBToND = 0;
}
+
+ if (fTriggerBiasCorrectionMBToOnePart) {
+ delete fTriggerBiasCorrectionMBToOnePart;
+ fTriggerBiasCorrectionMBToOnePart = 0;
+ }
}
//____________________________________________________________________
fTriggerBiasCorrectionMBToINEL->Divide();
fTriggerBiasCorrectionMBToNSD->Divide();
fTriggerBiasCorrectionMBToND->Divide();
+ fTriggerBiasCorrectionMBToOnePart->Divide();
}
//____________________________________________________________________
TList* collectionTriggerBiasMBToINEL = new TList;
TList* collectionTriggerBiasMBToNSD = new TList;
TList* collectionTriggerBiasMBToND = new TList;
+ TList* collectionTriggerBiasMBToOnePart = new TList;
Int_t count = 0;
while ((obj = iter->Next())) {
collectionTriggerBiasMBToINEL->Add(entry->fTriggerBiasCorrectionMBToINEL);
collectionTriggerBiasMBToNSD ->Add(entry->fTriggerBiasCorrectionMBToNSD);
collectionTriggerBiasMBToND ->Add(entry->fTriggerBiasCorrectionMBToND);
+ collectionTriggerBiasMBToOnePart ->Add(entry->fTriggerBiasCorrectionMBToOnePart);
count++;
}
fTriggerBiasCorrectionMBToINEL ->Merge(collectionTriggerBiasMBToINEL);
fTriggerBiasCorrectionMBToNSD ->Merge(collectionTriggerBiasMBToNSD);
fTriggerBiasCorrectionMBToND ->Merge(collectionTriggerBiasMBToND);
+ fTriggerBiasCorrectionMBToOnePart->Merge(collectionTriggerBiasMBToOnePart);
delete collectionNtrackToNparticle;
delete collectionVertexReco;
delete collectionTriggerBiasMBToINEL;
delete collectionTriggerBiasMBToNSD;
delete collectionTriggerBiasMBToND;
+ delete collectionTriggerBiasMBToOnePart;
return count+1;
}
fTriggerBiasCorrectionMBToINEL ->Add(aCorrectionsToAdd->GetTriggerBiasCorrectionINEL(),c);
fTriggerBiasCorrectionMBToNSD ->Add(aCorrectionsToAdd->GetTriggerBiasCorrectionNSD() ,c);
fTriggerBiasCorrectionMBToND ->Add(aCorrectionsToAdd->GetTriggerBiasCorrectionND() ,c);
-
+ fTriggerBiasCorrectionMBToOnePart ->Add(aCorrectionsToAdd->GetTriggerBiasCorrectionOnePart() ,c);
+}
+
+//____________________________________________________________________
+void AlidNdEtaCorrection::Scale(Float_t c)
+{
+ //
+ // scales all contained corrections
+ //
+
+ fTrack2ParticleCorrection ->Scale(c);
+ fVertexRecoCorrection ->Scale(c);
+ fTriggerBiasCorrectionMBToINEL ->Scale(c);
+ fTriggerBiasCorrectionMBToNSD ->Scale(c);
+ fTriggerBiasCorrectionMBToND ->Scale(c);
+ fTriggerBiasCorrectionMBToOnePart ->Scale(c);
}
//____________________________________________________________________
fTriggerBiasCorrectionMBToINEL ->Reset();
fTriggerBiasCorrectionMBToNSD ->Reset();
fTriggerBiasCorrectionMBToND ->Reset();
-
+ fTriggerBiasCorrectionMBToOnePart ->Reset();
}
fTriggerBiasCorrectionMBToINEL ->LoadHistograms();
fTriggerBiasCorrectionMBToNSD ->LoadHistograms();
fTriggerBiasCorrectionMBToND ->LoadHistograms();
+ fTriggerBiasCorrectionMBToOnePart ->LoadHistograms();
gDirectory->cd("..");
fTriggerBiasCorrectionMBToINEL->SaveHistograms();
fTriggerBiasCorrectionMBToNSD ->SaveHistograms();
fTriggerBiasCorrectionMBToND ->SaveHistograms();
+ fTriggerBiasCorrectionMBToOnePart->SaveHistograms();
gDirectory->cd("..");
}
fTriggerBiasCorrectionMBToINEL->DrawHistograms();
fTriggerBiasCorrectionMBToNSD ->DrawHistograms();
fTriggerBiasCorrectionMBToND ->DrawHistograms();
+ fTriggerBiasCorrectionMBToOnePart ->DrawHistograms();
}
//____________________________________________________________________
fTriggerBiasCorrectionMBToINEL->DrawOverview(canvasName);
fTriggerBiasCorrectionMBToNSD ->DrawOverview(canvasName);
fTriggerBiasCorrectionMBToND ->DrawOverview(canvasName);
+ fTriggerBiasCorrectionMBToOnePart ->DrawOverview(canvasName);
}
//____________________________________________________________________
fTriggerBiasCorrectionMBToINEL->GetTrackCorrection()->FillGene(vtx, eta, pt);
- if (processType != AliPWG0Helper::kSD )
+ if ((processType & AliPWG0Helper::kSD) == 0)
fTriggerBiasCorrectionMBToNSD->GetTrackCorrection()->FillGene(vtx, eta, pt);
- if (processType == AliPWG0Helper::kND )
+ if (processType & AliPWG0Helper::kND )
fTriggerBiasCorrectionMBToND->GetTrackCorrection()->FillGene(vtx, eta, pt);
+ if (processType & AliPWG0Helper::kOnePart)
+ fTriggerBiasCorrectionMBToOnePart->GetTrackCorrection()->FillGene(vtx, eta, pt);
+
if (!trigger)
return;
fTriggerBiasCorrectionMBToINEL->GetTrackCorrection()->FillMeas(vtx, eta, pt);
fTriggerBiasCorrectionMBToNSD->GetTrackCorrection()->FillMeas(vtx, eta, pt);
fTriggerBiasCorrectionMBToND->GetTrackCorrection()->FillMeas(vtx, eta, pt);
+ fTriggerBiasCorrectionMBToOnePart->GetTrackCorrection()->FillMeas(vtx, eta, pt);
fVertexRecoCorrection->GetTrackCorrection()->FillGene(vtx, eta, pt);
if (!vertex)
fTriggerBiasCorrectionMBToINEL->GetEventCorrection()->FillGene(vtx, n);
- if ( processType != AliPWG0Helper::kSD )
+ if ((processType & AliPWG0Helper::kSD) == 0)
fTriggerBiasCorrectionMBToNSD->GetEventCorrection()->FillGene(vtx, n);
- if (processType == AliPWG0Helper::kND )
+ if (processType & AliPWG0Helper::kND )
fTriggerBiasCorrectionMBToND->GetEventCorrection()->FillGene(vtx, n);
+ if (processType & AliPWG0Helper::kOnePart)
+ fTriggerBiasCorrectionMBToOnePart->GetEventCorrection()->FillGene(vtx, n);
+
if (!trigger)
return;
fTriggerBiasCorrectionMBToINEL->GetEventCorrection()->FillMeas(vtx, n);
fTriggerBiasCorrectionMBToNSD->GetEventCorrection()->FillMeas(vtx, n);
fTriggerBiasCorrectionMBToND->GetEventCorrection()->FillMeas(vtx, n);
+ fTriggerBiasCorrectionMBToOnePart->GetEventCorrection()->FillMeas(vtx, n);
fVertexRecoCorrection->GetEventCorrection()->FillGene(vtx, n);
if (!vertex)
fTriggerBiasCorrectionMBToINEL ->ReduceInformation();
fTriggerBiasCorrectionMBToNSD ->ReduceInformation();
fTriggerBiasCorrectionMBToND ->ReduceInformation();
+ fTriggerBiasCorrectionMBToOnePart ->ReduceInformation();
}
//____________________________________________________________________
case kINEL : return fTriggerBiasCorrectionMBToINEL;
case kNSD : return fTriggerBiasCorrectionMBToNSD;
case kND : return fTriggerBiasCorrectionMBToND;
+ case kOnePart : return fTriggerBiasCorrectionMBToOnePart;
}
return 0;
kVertexReco, // MB sample
kINEL,
kNSD,
- kND
+ kND,
+ kOnePart
};
AlidNdEtaCorrection();
AliCorrection* GetTriggerBiasCorrectionINEL() {return fTriggerBiasCorrectionMBToINEL;}
AliCorrection* GetTriggerBiasCorrectionNSD() {return fTriggerBiasCorrectionMBToNSD;}
AliCorrection* GetTriggerBiasCorrectionND() {return fTriggerBiasCorrectionMBToND;}
+ AliCorrection* GetTriggerBiasCorrectionOnePart() {return fTriggerBiasCorrectionMBToOnePart;}
AliCorrection* GetCorrection(CorrectionType correctionType);
void Reset(void);
void Add(AlidNdEtaCorrection* aCorrectionsToAdd, Float_t c=1);
+ void Scale(Float_t c);
void SaveHistograms();
Bool_t LoadHistograms(const Char_t* dir = 0);
AliCorrection* fTriggerBiasCorrectionMBToINEL; //-> handles the trigger bias MB->INEL, function of n and vtx_z
AliCorrection* fTriggerBiasCorrectionMBToNSD; //-> handles the trigger bias MB->NSD, function of n and vtx_z
AliCorrection* fTriggerBiasCorrectionMBToND; //-> handles the trigger bias MB->ND, function of n and vtx_z
+ AliCorrection* fTriggerBiasCorrectionMBToOnePart; //-> handles the trigger bias MB->OnePart, function of n and vtx_z
private:
AlidNdEtaCorrection(const AlidNdEtaCorrection&);
AlidNdEtaCorrection& operator=(const AlidNdEtaCorrection&);
- ClassDef(AlidNdEtaCorrection, 1)
+ ClassDef(AlidNdEtaCorrection, 2)
};
#endif
#include <TParticle.h>
#include <TParticlePDG.h>
#include <TDatabasePDG.h>
+#include <TRandom.h>
#include <AliLog.h>
#include <AliESDVertex.h>
#include <AliESDEvent.h>
+#include <AliESDRun.h>
#include <AliStack.h>
#include <AliHeader.h>
#include <AliGenEventHeader.h>
#include "dNdEta/dNdEtaAnalysis.h"
#include "dNdEta/AlidNdEtaCorrection.h"
#include "AliTriggerAnalysis.h"
+#include "AliPhysicsSelection.h"
ClassImp(AlidNdEtaCorrectionTask)
fFillPhi(kFALSE),
fDeltaPhiCut(-1),
fSymmetrize(kFALSE),
+ fMultAxisEta1(kFALSE),
+ fDiffTreatment(AliPWG0Helper::kMCFlags),
fSignMode(0),
fOnlyPrimaries(kFALSE),
fStatError(0),
+ fSystSkipParticles(kFALSE),
fEsdTrackCuts(0),
fdNdEtaCorrection(0),
fdNdEtaAnalysisMC(0),
for (Int_t i=0; i<8; i++)
fDeltaPhi[i] = 0;
+
+ AliLog::SetClassDebugLevel("AlidNdEtaCorrectionTask", AliLog::kWarning);
}
AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
fFillPhi(kFALSE),
fDeltaPhiCut(0),
fSymmetrize(kFALSE),
+ fMultAxisEta1(kFALSE),
+ fDiffTreatment(AliPWG0Helper::kMCFlags),
fSignMode(0),
fOnlyPrimaries(kFALSE),
fStatError(0),
+ fSystSkipParticles(kFALSE),
fEsdTrackCuts(0),
fdNdEtaCorrection(0),
fdNdEtaAnalysisMC(0),
for (Int_t i=0; i<8; i++)
fDeltaPhi[i] = 0;
+
+ AliLog::SetClassDebugLevel("AlidNdEtaCorrectionTask", AliLog::kWarning);
}
AlidNdEtaCorrectionTask::~AlidNdEtaCorrectionTask()
}
- fTemp1 = new TH2F("fTemp1", "fTemp1", 200, -20, 20, 200, -0.5, 199.5);
+ //fTemp1 = new TH2F("fTemp1", "fTemp1", 4, 0.5, 4.5, 101, -1.5, 99.5); // nsd study
+ fTemp1 = new TH2F("fTemp1", "fTemp1", 300, -15, 15, 80, -2.0, 2.0);
fOutput->Add(fTemp1);
- /*
- fTemp2 = new TH1F("fTemp2", "fTemp2", 2000, -5, 5);
+
+ fTemp2 = new TH2F("fTemp2", "fTemp2", 300, -15, 15, 80, -2.0, 2.0);
fOutput->Add(fTemp2);
- */
fVertexCorrelation = new TH2F("fVertexCorrelation", "fVertexCorrelation;MC z-vtx;ESD z-vtx", 120, -30, 30, 120, -30, 30);
fOutput->Add(fVertexCorrelation);
fOutput->Add(fVertexProfile);
fVertexShift = new TH1F("fVertexShift", "fVertexShift;(MC z-vtx - ESD z-vtx);Entries", 201, -2, 2);
fOutput->Add(fVertexShift);
- fVertexShiftNorm = new TH2F("fVertexShiftNorm", "fVertexShiftNorm;(MC z-vtx - ESD z-vtx) / #sigma_{ESD z-vtx};rec. tracks;Entries", 200, -100, 100, 100, -0.5, 99.5);
+ fVertexShiftNorm = new TH2F("fVertexShiftNorm", "fVertexShiftNorm;(MC z-vtx - ESD z-vtx);rec. tracks;Entries", 200, -100, 100, 100, -0.5, 99.5);
fOutput->Add(fVertexShiftNorm);
fEtaCorrelation = new TH2F("fEtaCorrelation", "fEtaCorrelation;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);
fOutput->Add(fDeltaPhi[i]);
}
- fEventStats = new TH2F("fEventStats", "fEventStats;event type;status;count", 106, -5.5, 100.5, 4, -0.5, 3.5);
+ fEventStats = new TH2F("fEventStats", "fEventStats;event type;status;count", 109, -6.5, 102.5, 4, -0.5, 3.5);
fOutput->Add(fEventStats);
fEventStats->GetXaxis()->SetBinLabel(1, "INEL"); // x = -5
fEventStats->GetXaxis()->SetBinLabel(2, "NSD"); // x = -4
fEventStats->GetXaxis()->SetBinLabel(3, "ND"); // x = -3
fEventStats->GetXaxis()->SetBinLabel(4, "SD"); // x = -2
fEventStats->GetXaxis()->SetBinLabel(5, "DD"); // x = -1
+
+ fEventStats->GetXaxis()->SetBinLabel(108, "INEL=0"); // x = -101
+ fEventStats->GetXaxis()->SetBinLabel(109, "INEL>0"); // x = -102
- for (Int_t i=0; i<100; i++)
+ for (Int_t i=-1; i<100; i++)
fEventStats->GetXaxis()->SetBinLabel(7+i, Form("%d", i));
fEventStats->GetYaxis()->SetBinLabel(1, "nothing");
if (fEsdTrackCuts)
{
fEsdTrackCuts->SetName("fEsdTrackCuts");
- fOutput->Add(fEsdTrackCuts);
+ // TODO like this we send an empty object back...
+ fOutput->Add(fEsdTrackCuts->Clone());
}
}
if (fStatError > 0)
Printf("WARNING: Statistical error evaluation active!");
- // trigger definition
- static AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis;
- Bool_t eventTriggered = triggerAnalysis->IsTriggerFired(fESD, fTrigger);
- //Printf("Trigger mask: %lld", fESD->GetTriggerMask());
+ AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+ if (!inputHandler)
+ {
+ Printf("ERROR: Could not receive input handler");
+ return;
+ }
+
+ Bool_t eventTriggered = inputHandler->IsEventSelected();
+ static AliTriggerAnalysis* triggerAnalysis = 0;
+ if (!triggerAnalysis)
+ {
+ AliPhysicsSelection* physicsSelection = dynamic_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
+ if (physicsSelection)
+ triggerAnalysis = physicsSelection->GetTriggerAnalysis();
+ }
+
+ if (eventTriggered)
+ eventTriggered = triggerAnalysis->IsTriggerFired(fESD, fTrigger);
+
if (!eventTriggered)
Printf("No trigger");
return;
}
- // get process type;
- Int_t processType = AliPWG0Helper::GetEventProcessType(header);
+ // get process type
+ Int_t processType = AliPWG0Helper::GetEventProcessType(fESD, header, stack, fDiffTreatment);
+
AliDebug(AliLog::kDebug+1, Form("Found process type %d", processType));
if (processType == AliPWG0Helper::kInvalidProcess)
- AliDebug(AliLog::kError, "Unknown process type.");
+ {
+ AliDebug(AliLog::kWarning, "Unknown process type. Setting to ND");
+ processType = AliPWG0Helper::kND;
+ }
// get the MC vertex
AliGenEventHeader* genHeader = header->GenEventHeader();
const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
if (vtxESD && AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
{
- eventVertex = kTRUE;
vtxESD->GetXYZ(vtx);
+ eventVertex = kTRUE;
+
+ // remove vertices outside +- 15 cm
+ if (TMath::Abs(vtx[2]) > 15)
+ {
+ eventVertex = kFALSE;
+ vtxESD = 0;
+ }
}
else
vtxESD = 0;
- // fill process type
- Int_t biny = (Int_t) eventTriggered + 2 * (Int_t) eventVertex;
- // INEL
- fEventStats->Fill(-5, biny);
- // NSD
- if (processType != AliPWG0Helper::kSD)
- fEventStats->Fill(-4, biny);
- // SD, ND, DD
- if (processType == AliPWG0Helper::kND)
- fEventStats->Fill(-3, biny);
- if (processType == AliPWG0Helper::kSD)
- fEventStats->Fill(-2, biny);
- if (processType == AliPWG0Helper::kDD)
- fEventStats->Fill(-1, biny);
-
// create list of (label, eta, pt) tuples
Int_t inputCount = 0;
Int_t* labelArr = 0;
Float_t* deltaPhiArr = 0;
if (fAnalysisMode & AliPWG0Helper::kSPD)
{
- // get tracklets
- const AliMultiplicity* mult = fESD->GetMultiplicity();
- if (!mult)
- {
- AliDebug(AliLog::kError, "AliMultiplicity not available");
- return;
- }
-
- labelArr = new Int_t[mult->GetNumberOfTracklets()];
- labelArr2 = new Int_t[mult->GetNumberOfTracklets()];
- etaArr = new Float_t[mult->GetNumberOfTracklets()];
- thirdDimArr = new Float_t[mult->GetNumberOfTracklets()];
- deltaPhiArr = new Float_t[mult->GetNumberOfTracklets()];
-
- // get multiplicity from SPD tracklets
- for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
+ if (vtxESD)
{
- //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
-
- Float_t phi = mult->GetPhi(i);
- if (phi < 0)
- phi += TMath::Pi() * 2;
- Float_t deltaPhi = mult->GetDeltaPhi(i);
-
- if (TMath::Abs(deltaPhi) > 1)
- printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
-
- if (fOnlyPrimaries)
- if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
+ // get tracklets
+ const AliMultiplicity* mult = fESD->GetMultiplicity();
+ if (!mult)
+ {
+ AliDebug(AliLog::kError, "AliMultiplicity not available");
+ return;
+ }
+
+ labelArr = new Int_t[mult->GetNumberOfTracklets()];
+ labelArr2 = new Int_t[mult->GetNumberOfTracklets()];
+ etaArr = new Float_t[mult->GetNumberOfTracklets()];
+ thirdDimArr = new Float_t[mult->GetNumberOfTracklets()];
+ deltaPhiArr = new Float_t[mult->GetNumberOfTracklets()];
+
+ Bool_t foundInEta10 = kFALSE;
+
+ // get multiplicity from SPD tracklets
+ for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
+ {
+ //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
+
+ Float_t phi = mult->GetPhi(i);
+ if (phi < 0)
+ phi += TMath::Pi() * 2;
+ Float_t deltaPhi = mult->GetDeltaPhi(i);
+
+ if (TMath::Abs(deltaPhi) > 1)
+ printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
+
+ if (fOnlyPrimaries)
+ if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
+ continue;
+
+ if (fDeltaPhiCut > 0 && (TMath::Abs(deltaPhi) > fDeltaPhiCut || TMath::Abs(mult->GetDeltaTheta(i)) > fDeltaPhiCut / 0.08 * 0.025))
continue;
-
- if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut)
- continue;
+
+ if (fSystSkipParticles && gRandom->Uniform() < 0.0153)
+ {
+ Printf("Skipped tracklet!");
+ continue;
+ }
+
+ // TEST exclude potentially inefficient phi region
+ //if (phi > 5.70 || phi < 0.06)
+ // continue;
+
+ // we have to repeat the trigger here, because the tracklet might have been kicked out fSystSkipParticles
+ if (TMath::Abs(mult->GetEta(i)) < 1)
+ foundInEta10 = kTRUE;
+
+ etaArr[inputCount] = mult->GetEta(i);
+ if (fSymmetrize)
+ etaArr[inputCount] = TMath::Abs(etaArr[inputCount]);
+ labelArr[inputCount] = mult->GetLabel(i, 0);
+ labelArr2[inputCount] = mult->GetLabel(i, 1);
+ thirdDimArr[inputCount] = phi;
+ deltaPhiArr[inputCount] = deltaPhi;
+ ++inputCount;
+ }
- etaArr[inputCount] = mult->GetEta(i);
- if (fSymmetrize)
- etaArr[inputCount] = TMath::Abs(etaArr[inputCount]);
- labelArr[inputCount] = mult->GetLabel(i, 0);
- labelArr2[inputCount] = mult->GetLabel(i, 1);
- thirdDimArr[inputCount] = phi;
- deltaPhiArr[inputCount] = deltaPhi;
- ++inputCount;
+ /*
+ for (Int_t i=0; i<mult->GetNumberOfSingleClusters(); ++i)
+ {
+ if (TMath::Abs(TMath::Log(TMath::Tan(mult->GetThetaSingle(i)/2.))) < 1);
+ {
+ foundInEta10 = kTRUE;
+ break;
+ }
+ }
+ */
+
+ if (fSystSkipParticles && (fTrigger & AliTriggerAnalysis::kOneParticle) && !foundInEta10)
+ eventTriggered = kFALSE;
}
}
else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
return;
}
+ Bool_t foundInEta10 = kFALSE;
+
if (vtxESD)
{
// get multiplicity from ESD tracks
continue;
}
- //Printf("status is: %u", esdTrack->GetStatus());
+ if (esdTrack->Pt() < 0.15)
+ continue;
if (fOnlyPrimaries)
{
if (stack->IsPhysicalPrimary(label) == kFALSE)
continue;
- }
+ }
+
+ // INEL>0 trigger
+ if (TMath::Abs(esdTrack->Eta()) < 1 && esdTrack->Pt() > 0.15)
+ foundInEta10 = kTRUE;
etaArr[inputCount] = esdTrack->Eta();
if (fSymmetrize)
delete list;
- if (eventTriggered)
+ // TODO this code crashes for TPCITS because particles are requested from the stack for some labels that are out of bound
+ if (0 && eventTriggered)
{
// collect values for primaries and secondaries
for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++)
}
}
}
+
+ if (!foundInEta10)
+ eventTriggered = kFALSE;
}
else
return;
+ // fill process type
+ Int_t biny = (Int_t) eventTriggered + 2 * (Int_t) eventVertex;
+ // INEL
+ fEventStats->Fill(-6, biny);
+ // NSD
+ if (processType != AliPWG0Helper::kSD)
+ fEventStats->Fill(-5, biny);
+ // SD, ND, DD
+ if (processType == AliPWG0Helper::kND)
+ fEventStats->Fill(-4, biny);
+ if (processType == AliPWG0Helper::kSD)
+ fEventStats->Fill(-3, biny);
+ if (processType == AliPWG0Helper::kDD)
+ fEventStats->Fill(-2, biny);
+
// loop over mc particles
Int_t nPrim = stack->GetNprimary();
Int_t nAccepted = 0;
+ Bool_t oneParticleEvent = kFALSE;
+ for (Int_t iMc = 0; iMc < nPrim; ++iMc)
+ {
+ //Printf("Getting particle %d", iMc);
+ TParticle* particle = stack->Particle(iMc);
+
+ if (!particle)
+ continue;
+
+ if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
+ continue;
+
+ if (TMath::Abs(particle->Eta()) < 1.0)
+ {
+ oneParticleEvent = kTRUE;
+ break;
+ }
+ }
+
+ if (TMath::Abs(vtxMC[2]) < 5.5)
+ {
+ if (oneParticleEvent)
+ fEventStats->Fill(102, biny);
+ else
+ fEventStats->Fill(101, biny);
+ }
+
for (Int_t iMc = 0; iMc < nPrim; ++iMc)
{
//Printf("Getting particle %d", iMc);
if (SignOK(particle->GetPDG()) == kFALSE)
continue;
-
+
if (fPIDParticles && TMath::Abs(particle->Eta()) < 1.0)
fPIDParticles->Fill(particle->GetPdgCode());
//fTemp1->Fill(eta);
//fTemp2->Fill(y);
- fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
+ Int_t processType2 = processType;
+ if (oneParticleEvent)
+ processType2 |= AliPWG0Helper::kOnePart;
+ fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
if (fOption.Contains("process-types"))
{
// non diffractive
if (processType==AliPWG0Helper::kND)
- fdNdEtaCorrectionSpecial[0]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
+ fdNdEtaCorrectionSpecial[0]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
// single diffractive
if (processType==AliPWG0Helper::kSD)
- fdNdEtaCorrectionSpecial[1]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
+ fdNdEtaCorrectionSpecial[1]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
// double diffractive
if (processType==AliPWG0Helper::kDD)
- fdNdEtaCorrectionSpecial[2]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
+ fdNdEtaCorrectionSpecial[2]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
}
if (fOption.Contains("particle-species"))
case 2212: id = 2; break;
default: id = 3; break;
}
- fdNdEtaCorrectionSpecial[id]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
+ fdNdEtaCorrectionSpecial[id]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
}
if (eventTriggered)
nAccepted++;
}
- fEventStats->Fill(AliPWG0Helper::GetLastProcessType(), biny);
+ if (AliPWG0Helper::GetLastProcessType() >= -1)
+ fEventStats->Fill(AliPWG0Helper::GetLastProcessType(), biny);
fMultAll->Fill(nAccepted);
if (eventTriggered) {
}
Int_t nEta05 = 0;
+ Int_t nEta10 = 0;
for (Int_t i=0; i<inputCount; ++i)
{
if (TMath::Abs(etaArr[i]) < 0.5)
nEta05++;
+ if (TMath::Abs(etaArr[i]) < 1.0)
+ nEta10++;
+ }
+ for (Int_t i=0; i<inputCount; ++i)
+ {
Int_t label = labelArr[i];
Int_t label2 = labelArr2[i];
continue;
}
- fPIDTracks->Fill(particle->GetPdgCode());
+ if (TMath::Abs(particle->Eta()) < 1.0)
+ fPIDTracks->Fill(particle->GetPdgCode());
// find particle that is filled in the correction map
// this should be the particle that has been reconstructed
{
if (fAnalysisMode & AliPWG0Helper::kSPD && !fFillPhi)
{
- thirdDim = inputCount;
+ thirdDim = (fMultAxisEta1) ? nEta10 : inputCount;
}
else
thirdDim = thirdDimArr[i];
}
if (fillTrack)
+ {
fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], eta, thirdDim);
-
+ fTemp2->Fill(vtxMC[2], eta);
+ }
+
// eta comparison for tracklets with the same label (others are background)
if (label == label2)
{
{
Double_t diff = vtxMC[2] - vtx[2];
fVertexShift->Fill(diff);
- if (vtxESD->GetZRes() > 0)
- fVertexShiftNorm->Fill(diff / vtxESD->GetZRes(), inputCount);
-
+
fVertexCorrelation->Fill(vtxMC[2], vtx[2]);
fVertexCorrelationShift->Fill(vtxMC[2], vtxMC[2] - vtx[2], inputCount);
fVertexProfile->Fill(vtxMC[2], vtxMC[2] - vtx[2]);
-
- fTemp1->Fill(vtxMC[2], nEta05);
+
+ if (vtxESD->IsFromVertexerZ() && inputCount > 0)
+ fVertexShiftNorm->Fill(diff, vtxESD->GetNContributors());
}
+ Int_t multAxis = inputCount;
+ if (fMultAxisEta1)
+ multAxis = nEta10;
+
if (eventTriggered && eventVertex)
{
- fdNdEtaAnalysisMC->FillEvent(vtxMC[2], inputCount);
- fdNdEtaAnalysisESD->FillEvent(vtxMC[2], inputCount);
+ fdNdEtaAnalysisMC->FillEvent(vtxMC[2], multAxis);
+ fdNdEtaAnalysisESD->FillEvent(vtxMC[2], multAxis);
}
- // stuff regarding the vertex reco correction and trigger bias correction
- fdNdEtaCorrection->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
+ Int_t processType2 = processType;
+ if (oneParticleEvent)
+ processType2 |= AliPWG0Helper::kOnePart;
+
+ // stuff regarding the vertex reco correction and trigger bias correction
+ fdNdEtaCorrection->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2);
if (fOption.Contains("process-types"))
{
// non diffractive
if (processType == AliPWG0Helper::kND)
- fdNdEtaCorrectionSpecial[0]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
+ fdNdEtaCorrectionSpecial[0]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2);
// single diffractive
if (processType == AliPWG0Helper::kSD)
- fdNdEtaCorrectionSpecial[1]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
+ fdNdEtaCorrectionSpecial[1]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2);
// double diffractive
if (processType == AliPWG0Helper::kDD)
- fdNdEtaCorrectionSpecial[2]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
+ fdNdEtaCorrectionSpecial[2]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2);
}
if (fOption.Contains("particle-species"))
for (Int_t id=0; id<4; id++)
- fdNdEtaCorrectionSpecial[id]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
+ fdNdEtaCorrectionSpecial[id]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2);
if (etaArr)
delete[] etaArr;
void SetFillPhi(Bool_t flag = kTRUE) { fFillPhi = flag; }
void SetDeltaPhiCut(Float_t cut) { fDeltaPhiCut = cut; }
void SetSymmetrize(Bool_t flag = kTRUE) { fSymmetrize = flag; }
+ void SetMultAxisEta1(Bool_t flag = kTRUE) { fMultAxisEta1 = flag; }
+ void SetDiffTreatment(AliPWG0Helper::DiffTreatment diffTreatment) { fDiffTreatment = diffTreatment; }
+ void SetSkipParticles(Bool_t flag = kTRUE) { fSystSkipParticles = flag; }
void SetOption(const char* opt) { fOption = opt; }
Bool_t fFillPhi; // if true phi is filled as 3rd coordinate in all maps
Float_t fDeltaPhiCut; // cut in delta phi (only SPD)
Bool_t fSymmetrize; // move all negative to positive eta
+ Bool_t fMultAxisEta1; // restrict multiplicity count to |eta| < 1
+ AliPWG0Helper::DiffTreatment fDiffTreatment; // how to identify SD events (see AliPWG0Helper::GetEventProcessType)
Int_t fSignMode; // if 0 process all particles, if +-1 process only particles with that sign
Bool_t fOnlyPrimaries; // only process primaries (syst. studies)
Int_t fStatError; // statistical error evaluation: if set to 1 we only count unique primaries (binomial errors are valid), for 2 all the rest
+ Bool_t fSystSkipParticles; // if true skips particles (systematic study)
AliESDtrackCuts* fEsdTrackCuts; // Object containing the parameters of the esd track cuts
fUseMCKine(kFALSE),
fCheckEventType(kFALSE),
fSymmetrize(kFALSE),
+ fMultAxisEta1(kFALSE),
+ fDiffTreatment(AliPWG0Helper::kMCFlags),
fEsdTrackCuts(0),
- fPhysicsSelection(0),
fTriggerAnalysis(0),
fdNdEtaAnalysisESD(0),
fMult(0),
fdNdEtaAnalysis(0),
fdNdEtaAnalysisND(0),
fdNdEtaAnalysisNSD(0),
+ fdNdEtaAnalysisOnePart(0),
fdNdEtaAnalysisTr(0),
fdNdEtaAnalysisTrVtx(0),
fdNdEtaAnalysisTracks(0),
fFiredChips(0),
fTrackletsVsClusters(0),
fTrackletsVsUnassigned(0),
- fTriggerVsTime(0),
fStats(0),
fStats2(0)
{
for (Int_t i=0; i<3; ++i)
{
- fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -6, 6);
+ fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -3, 3);
fPartEta[i]->Sumw2();
fOutput->Add(fPartEta[i]);
}
fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 800, -40, 40);
fOutput->Add(fEvents);
- Float_t resMax = (fAnalysisMode & AliPWG0Helper::kSPD) ? 0.2 : 2;
- fVertexResolution = new TH1F("dndeta_vertex_resolution_z", "dndeta_vertex_resolution_z", 1000, 0, resMax);
+ fVertexResolution = new TH2F("dndeta_vertex_resolution_z", ";z resolution;#Delta #phi", 1000, 0, 2, 100, 0, 0.2);
fOutput->Add(fVertexResolution);
fPhi = new TH1F("fPhi", "fPhi;#phi in rad.;count", 720, 0, 2 * TMath::Pi());
fOutput->Add(fPhi);
- fEtaPhi = new TH2F("fEtaPhi", "fEtaPhi;#eta;#phi in rad.;count", 80, -4, 4, 18*5, 0, 2 * TMath::Pi());
+ fEtaPhi = new TH2F("fEtaPhi", "fEtaPhi;#eta;#phi in rad.;count", 80, -4, 4, 18*20, 0, 2 * TMath::Pi());
fOutput->Add(fEtaPhi);
- fTriggerVsTime = new TGraph; //TH1F("fTriggerVsTime", "fTriggerVsTime;trigger time;count", 100, 0, 100);
- fTriggerVsTime->SetName("fTriggerVsTime");
- fTriggerVsTime->GetXaxis()->SetTitle("trigger time");
- fTriggerVsTime->GetYaxis()->SetTitle("count");
- fOutput->Add(fTriggerVsTime);
-
fStats = new TH1F("fStats", "fStats", 5, 0.5, 5.5);
fStats->GetXaxis()->SetBinLabel(1, "vertexer 3d");
fStats->GetXaxis()->SetBinLabel(2, "vertexer z");
fStats2 = new TH2F("fStats2", "fStats2", 7, -0.5, 6.5, 10, -0.5, 9.5);
fStats2->GetXaxis()->SetBinLabel(1, "No trigger");
- fStats2->GetXaxis()->SetBinLabel(2, "Splash identification");
- fStats2->GetXaxis()->SetBinLabel(3, "No Vertex");
- fStats2->GetXaxis()->SetBinLabel(4, "|z-vtx| > 10");
- fStats2->GetXaxis()->SetBinLabel(5, "0 tracklets");
- fStats2->GetXaxis()->SetBinLabel(6, "V0 Veto");
+ fStats2->GetXaxis()->SetBinLabel(2, "No Vertex");
+ fStats2->GetXaxis()->SetBinLabel(3, "Vtx res. too bad");
+ fStats2->GetXaxis()->SetBinLabel(4, "|z-vtx| > 15");
+ fStats2->GetXaxis()->SetBinLabel(5, "|z-vtx| > 10");
+ fStats2->GetXaxis()->SetBinLabel(6, "0 tracklets");
fStats2->GetXaxis()->SetBinLabel(7, "Selected");
fStats2->GetYaxis()->SetBinLabel(1, "n/a");
if (fAnalysisMode & AliPWG0Helper::kSPD)
{
- fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 500, -0.2, 0.2);
+ fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 500, -0.16, 0.16);
fOutput->Add(fDeltaPhi);
- fDeltaTheta = new TH1F("fDeltaTheta", "fDeltaTheta;#Delta #theta;Entries", 500, -0.2, 0.2);
+ fDeltaTheta = new TH1F("fDeltaTheta", "fDeltaTheta;#Delta #theta;Entries", 500, -0.05, 0.05);
fOutput->Add(fDeltaTheta);
fFiredChips = new TH2F("fFiredChips", "fFiredChips;Chips L1 + L2;tracklets", 1201, -0.5, 1201, 50, -0.5, 49.5);
fOutput->Add(fFiredChips);
fdNdEtaAnalysisNSD = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD", fAnalysisMode);
fOutput->Add(fdNdEtaAnalysisNSD);
+ fdNdEtaAnalysisOnePart = new dNdEtaAnalysis("dndetaOnePart", "dndetaOnePart", fAnalysisMode);
+ fOutput->Add(fdNdEtaAnalysisOnePart);
+
fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr", fAnalysisMode);
fOutput->Add(fdNdEtaAnalysisTr);
fOutput->Add(fEsdTrackCuts);
}
- if (!fPhysicsSelection)
- fPhysicsSelection = new AliPhysicsSelection;
-
- fPhysicsSelection->SetName("AliPhysicsSelection_outputlist"); // to prevent conflict with object that is automatically streamed back
//AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kDebug);
- //fOutput->Add(fPhysicsSelection);
-
- fTriggerAnalysis = new AliTriggerAnalysis;
- fTriggerAnalysis->EnableHistograms();
- fTriggerAnalysis->SetSPDGFOThreshhold(2);
- fOutput->Add(fTriggerAnalysis);
}
void AlidNdEtaTask::Exec(Option_t*)
return;
}
-// if (fCheckEventType)
-// eventTriggered = fPhysicsSelection->IsCollisionCandidate(fESD);
-
- // check event type (should be PHYSICS = 7)
- if (fCheckEventType)
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
+ if (!inputHandler)
{
- UInt_t eventType = esdHeader->GetEventType();
- if (eventType != 7)
- {
- Printf("Skipping event because it is of type %d", eventType);
- return;
- }
-
- //Printf("Trigger classes: %s:", fESD->GetFiredTriggerClasses().Data());
-
- Bool_t accept = kTRUE;
- if (fRequireTriggerClass.Length() > 0 && !fESD->IsTriggerClassFired(fRequireTriggerClass))
- accept = kFALSE;
- if (fRejectTriggerClass.Length() > 0 && fESD->IsTriggerClassFired(fRejectTriggerClass))
- accept = kFALSE;
-
- if (!accept)
- {
- Printf("Skipping event because it does not have the correct trigger class(es): %s", fESD->GetFiredTriggerClasses().Data());
- return;
- }
-
- fStats->Fill(4);
+ Printf("ERROR: Could not receive input handler");
+ return;
}
- fTriggerAnalysis->FillHistograms(fESD);
+ eventTriggered = inputHandler->IsEventSelected();
+
+ if (!fTriggerAnalysis)
+ {
+ AliPhysicsSelection* physicsSelection = dynamic_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
+ if (physicsSelection)
+ fTriggerAnalysis = physicsSelection->GetTriggerAnalysis();
+ }
+
+ if (eventTriggered)
+ eventTriggered = fTriggerAnalysis->IsTriggerFired(fESD, fTrigger);
AliTriggerAnalysis::V0Decision v0A = fTriggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kASide, kFALSE);
AliTriggerAnalysis::V0Decision v0C = fTriggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kCSide, kFALSE);
if (v0A == AliTriggerAnalysis::kV0BB && v0C == AliTriggerAnalysis::kV0BG) vZero = 8;
if (v0A == AliTriggerAnalysis::kV0BG && v0C == AliTriggerAnalysis::kV0BB) vZero = 9;
}
+
+ if (vZero == 0)
+ Printf("VZERO: %d %d %d %d", vZero, eventTriggered, v0A, v0C);
Bool_t filled = kFALSE;
- // trigger
- eventTriggered = fTriggerAnalysis->IsTriggerFired(fESD, fTrigger);
-
if (!eventTriggered)
{
fStats2->Fill(0.0, vZero);
filled = kTRUE;
}
- if (v0A == AliTriggerAnalysis::kV0BG || v0C == AliTriggerAnalysis::kV0BG)
- eventTriggered = kFALSE;
-
if (eventTriggered)
- {
fStats->Fill(3);
- }
- if (fCheckEventType)
+ /*
+ // ITS cluster tree
+ AliESDInputHandlerRP* handlerRP = dynamic_cast<AliESDInputHandlerRP*> (AliAnalysisManager::GetAnalys
+isManager()->GetInputEventHandler());
+ if (handlerRP)
{
- /*const Int_t kMaxEvents = 1;
- UInt_t maskedEvents[kMaxEvents][2] = {
- {-1, -1}
- };
-
- for (Int_t i=0; i<kMaxEvents; i++)
- {
- if (fESD->GetPeriodNumber() == maskedEvents[i][0] && fESD->GetOrbitNumber() == maskedEvents[i][1])
- {
- Printf("Skipping event because it is masked: period: %d orbit: %x", fESD->GetPeriodNumber(), fESD->GetOrbitNumber());
- if (!filled)
+ TTree* itsClusterTree = handlerRP->GetTreeR("ITS");
+ if (!itsClusterTree)
+ return;
+
+ TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
+ TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
+
+ itsClusterBranch->SetAddress(&itsClusters);
+
+ Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
+
+ Int_t totalClusters = 0;
+
+ // loop over the its subdetectors
+ for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
+
+ if (!itsClusterTree->GetEvent(iIts))
+ continue;
+
+ Int_t nClusters = itsClusters->GetEntriesFast();
+ totalClusters += nClusters;
+
+ #ifdef FULLALIROOT
+ if (fAnalysisMode & AliPWG0Helper::kSPD)
{
- fStats2->Fill(1, vZero);
- filled = kTRUE;
- }
- return;
- }
- }*/
-
- // ITS cluster tree
- AliESDInputHandlerRP* handlerRP = dynamic_cast<AliESDInputHandlerRP*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
- if (handlerRP)
- {
- TTree* itsClusterTree = handlerRP->GetTreeR("ITS");
- if (!itsClusterTree)
- return;
-
- TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
- TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
-
- itsClusterBranch->SetAddress(&itsClusters);
-
- Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
-
- Int_t totalClusters = 0;
-
- // loop over the its subdetectors
- for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
-
- if (!itsClusterTree->GetEvent(iIts))
- continue;
-
- Int_t nClusters = itsClusters->GetEntriesFast();
- totalClusters += nClusters;
-
- #ifdef FULLALIROOT
- if (fAnalysisMode & AliPWG0Helper::kSPD)
- {
- // loop over clusters
- while (nClusters--) {
- AliITSRecPoint* cluster = (AliITSRecPoint*) itsClusters->UncheckedAt(nClusters);
-
- Int_t layer = cluster->GetLayer();
-
- if (layer > 1)
- continue;
-
- Float_t xyz[3] = {0., 0., 0.};
- cluster->GetGlobalXYZ(xyz);
-
- Float_t phi = TMath::Pi() + TMath::ATan2(-xyz[1], -xyz[0]);
- Float_t z = xyz[2];
-
- fZPhi[layer]->Fill(z, phi);
- fModuleMap->Fill(layer * 80 + cluster->GetDetectorIndex());
- }
+ // loop over clusters
+ while (nClusters--) {
+ AliITSRecPoint* cluster = (AliITSRecPoint*) itsClusters->UncheckedAt(nClusters);
+
+ Int_t layer = cluster->GetLayer();
+
+ if (layer > 1)
+ continue;
+
+ Float_t xyz[3] = {0., 0., 0.};
+ cluster->GetGlobalXYZ(xyz);
+
+ Float_t phi = TMath::Pi() + TMath::ATan2(-xyz[1], -xyz[0]);
+ Float_t z = xyz[2];
+
+ fZPhi[layer]->Fill(z, phi);
+ fModuleMap->Fill(layer * 80 + cluster->GetDetectorIndex());
}
- #endif
- }
-
- const AliMultiplicity* mult = fESD->GetMultiplicity();
- if (!mult)
- return;
-
- fTrackletsVsClusters->Fill(mult->GetNumberOfTracklets(), totalClusters);
-
- Int_t limit = 80 + mult->GetNumberOfTracklets() * 220/20;
-
- if (totalClusters > limit)
- {
- Printf("Skipping event because %d clusters is above limit of %d from %d tracklets: Period number: %d Orbit number: %x", totalClusters, limit, mult->GetNumberOfTracklets(), fESD->GetPeriodNumber(), fESD->GetOrbitNumber());
- if (!filled)
- {
- fStats2->Fill(1, vZero);
- filled = kTRUE;
}
- return; // TODO we skip this also for the MC. not good...
- }
+ #endif
}
}
-
+ */
+
vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
if (!vtxESD)
{
if (!filled)
{
- fStats2->Fill(2, vZero);
+ fStats2->Fill(1, vZero);
filled = kTRUE;
}
}
else
{
+ if (!AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
+ {
+ if (!filled)
+ {
+ fStats2->Fill(2, vZero);
+ filled = kTRUE;
+ }
+ }
+
Double_t vtx[3];
vtxESD->GetXYZ(vtx);
- if (TMath::Abs(vtx[2]) > 10)
+ // try to compare spd vertex and vertexer from tracks
+ // remove vertices outside +- 15 cm
+ if (TMath::Abs(vtx[2]) > 15)
{
if (!filled)
{
filled = kTRUE;
}
}
+
+ if (TMath::Abs(vtx[2]) > 10)
+ {
+ if (!filled)
+ {
+ fStats2->Fill(4, vZero);
+ filled = kTRUE;
+ }
+ }
const AliMultiplicity* mult = fESD->GetMultiplicity();
if (!mult)
{
if (!filled)
{
- fStats2->Fill(4, vZero);
+ fStats2->Fill(5, vZero);
filled = kTRUE;
}
}
}
- if (fCheckEventType)
- {
- if (vZero >= 5)
- {
- if (!filled)
- fStats2->Fill(5, vZero);
- return;
- }
- }
-
if (!filled)
+ {
fStats2->Fill(6, vZero);
-
- //Printf("Skipping event %d: Period number: %d Orbit number: %x", decision, fESD->GetPeriodNumber(), fESD->GetOrbitNumber());
+ //Printf("File: %s, IEV: %d, TRG: ---, Orbit: 0x%x, Period: %d, BC: %d", ((TTree*) GetInputData(0))->GetCurrentFile()->GetName(), fESD->GetEventNumberInFile(), fESD->GetOrbitNumber(),fESD->GetPeriodNumber(),fESD->GetBunchCrossNumber());
+ }
if (eventTriggered)
fStats->Fill(3);
// get the ESD vertex
vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
- //Printf("vertex is: %s", fESD->GetPrimaryVertex()->GetTitle());
-
Double_t vtx[3];
// fill z vertex resolution
if (vtxESD)
{
- fVertexResolution->Fill(vtxESD->GetZRes());
+ if (strcmp(vtxESD->GetTitle(), "vertexer: Z") == 0)
+ fVertexResolution->Fill(vtxESD->GetZRes(), vtxESD->GetDispersion());
+ else
+ fVertexResolution->Fill(vtxESD->GetZRes(), 0);
+
//if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0)
{
fVertex->Fill(vtxESD->GetXv(), vtxESD->GetYv(), vtxESD->GetZv());
if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0)
{
fStats->Fill(1);
- if (fCheckEventType && TMath::Abs(vtx[0] > 0.3))
+ if (fCheckEventType && TMath::Abs(vtx[0]) > 0.3)
{
Printf("Suspicious x-vertex x=%f y=%f z=%f (period: %d orbit %x)", vtx[0], vtx[1], vtx[2], fESD->GetPeriodNumber(), fESD->GetOrbitNumber());
}
Printf("Suspicious y-vertex x=%f y=%f z=%f (period: %d orbit %x)", vtx[0], vtx[1], vtx[2], fESD->GetPeriodNumber(), fESD->GetOrbitNumber());
}
}
- else if (strcmp(vtxESD->GetTitle(), "vertexer: Z") == 0)
- {
+
+ if (strcmp(vtxESD->GetTitle(), "vertexer: Z") == 0)
fStats->Fill(2);
- }
+
+ // remove vertices outside +- 15 cm
+ if (TMath::Abs(vtx[2]) > 15)
+ vtxESD = 0;
}
else
vtxESD = 0;
return;
}
}
-
+
// create list of (label, eta, pt) tuples
Int_t inputCount = 0;
Int_t* labelArr = 0;
Float_t* thirdDimArr = 0;
if (fAnalysisMode & AliPWG0Helper::kSPD)
{
- // get tracklets
- const AliMultiplicity* mult = fESD->GetMultiplicity();
- if (!mult)
- {
- AliDebug(AliLog::kError, "AliMultiplicity not available");
- return;
- }
-
- labelArr = new Int_t[mult->GetNumberOfTracklets()];
- etaArr = new Float_t[mult->GetNumberOfTracklets()];
- thirdDimArr = new Float_t[mult->GetNumberOfTracklets()];
-
- // get multiplicity from SPD tracklets
- for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
+ if (vtxESD)
{
- //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
-
- if (fOnlyPrimaries)
- if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
- continue;
-
- Float_t deltaPhi = mult->GetDeltaPhi(i);
- // prevent values to be shifted by 2 Pi()
- if (deltaPhi < -TMath::Pi())
- deltaPhi += TMath::Pi() * 2;
- if (deltaPhi > TMath::Pi())
- deltaPhi -= TMath::Pi() * 2;
-
- if (TMath::Abs(deltaPhi) > 1)
- printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
-
- Int_t label = mult->GetLabel(i, 0);
- Float_t eta = mult->GetEta(i);
-
- // control histograms
- Float_t phi = mult->GetPhi(i);
- if (phi < 0)
- phi += TMath::Pi() * 2;
- fPhi->Fill(phi);
- fEtaPhi->Fill(eta, phi);
-
-// if (deltaPhi < 0.01)
-// {
-// // layer 0
-// Float_t z = vtx[2] + 3.9 / TMath::Tan(2 * TMath::ATan(TMath::Exp(- eta)));
-// fZPhi[0]->Fill(z, phi);
-// // layer 1
-// z = vtx[2] + 7.6 / TMath::Tan(2 * TMath::ATan(TMath::Exp(- eta)));
-// fZPhi[1]->Fill(z, phi);
-// }
-
- if (vtxESD && TMath::Abs(vtx[2]) < 10)
+ // get tracklets
+ const AliMultiplicity* mult = fESD->GetMultiplicity();
+ if (!mult)
{
- fDeltaPhi->Fill(deltaPhi);
- fDeltaTheta->Fill(mult->GetDeltaTheta(i));
+ AliDebug(AliLog::kError, "AliMultiplicity not available");
+ return;
+ }
+
+ Int_t arrayLength = mult->GetNumberOfTracklets();
+ if (fAnalysisMode & AliPWG0Helper::kSPDOnlyL0)
+ arrayLength += mult->GetNumberOfSingleClusters();
+
+ labelArr = new Int_t[arrayLength];
+ etaArr = new Float_t[arrayLength];
+ thirdDimArr = new Float_t[arrayLength];
+
+ // get multiplicity from SPD tracklets
+ for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
+ {
+ //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
+
+ if (fOnlyPrimaries)
+ if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
+ continue;
+
+ Float_t deltaPhi = mult->GetDeltaPhi(i);
+ // prevent values to be shifted by 2 Pi()
+ if (deltaPhi < -TMath::Pi())
+ deltaPhi += TMath::Pi() * 2;
+ if (deltaPhi > TMath::Pi())
+ deltaPhi -= TMath::Pi() * 2;
+
+ if (TMath::Abs(deltaPhi) > 1)
+ printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
+
+ Int_t label = mult->GetLabel(i, 0);
+ Float_t eta = mult->GetEta(i);
+
+ // control histograms
+ Float_t phi = mult->GetPhi(i);
+ if (phi < 0)
+ phi += TMath::Pi() * 2;
+
+ // TEST exclude probably inefficient phi region
+ //if (phi > 5.70 || phi < 0.06)
+ // continue;
+
+ fPhi->Fill(phi);
+
+ if (vtxESD && TMath::Abs(vtx[2]) < 10)
+ {
+ fEtaPhi->Fill(eta, phi);
+ fDeltaPhi->Fill(deltaPhi);
+ fDeltaTheta->Fill(mult->GetDeltaTheta(i));
+ }
+
+ if (fDeltaPhiCut > 0 && (TMath::Abs(deltaPhi) > fDeltaPhiCut || TMath::Abs(mult->GetDeltaTheta(i)) > fDeltaPhiCut / 0.08 * 0.025))
+ continue;
+
+ if (fUseMCKine)
+ {
+ if (label > 0)
+ {
+ TParticle* particle = stack->Particle(label);
+ eta = particle->Eta();
+ phi = particle->Phi();
+ }
+ else
+ Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found");
+ }
+
+ if (fSymmetrize)
+ eta = TMath::Abs(eta);
+
+ etaArr[inputCount] = eta;
+ labelArr[inputCount] = label;
+ thirdDimArr[inputCount] = phi;
+ ++inputCount;
}
- if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut)
- continue;
-
- if (fUseMCKine)
+ if (fAnalysisMode & AliPWG0Helper::kSPDOnlyL0)
{
- if (label > 0)
+ // get additional clusters from L0
+ for (Int_t i=0; i<mult->GetNumberOfSingleClusters(); ++i)
{
- TParticle* particle = stack->Particle(label);
- eta = particle->Eta();
- phi = particle->Phi();
+ etaArr[inputCount] = -TMath::Log(TMath::Tan(mult->GetThetaSingle(i)/2.));
+ labelArr[inputCount] = -1;
+ thirdDimArr[inputCount] = mult->GetPhiSingle(i);
+
+ ++inputCount;
}
- else
- Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found");
}
- if (fSymmetrize)
- eta = TMath::Abs(eta);
-
- etaArr[inputCount] = eta;
- labelArr[inputCount] = label;
- thirdDimArr[inputCount] = phi;
- ++inputCount;
- }
-
- if (!fFillPhi)
- {
- // fill multiplicity in third axis
- for (Int_t i=0; i<inputCount; ++i)
- thirdDimArr[i] = inputCount;
+ if (!fFillPhi)
+ {
+ // fill multiplicity in third axis
+ for (Int_t i=0; i<inputCount; ++i)
+ thirdDimArr[i] = inputCount;
+ }
+
+ Int_t firedChips = mult->GetNumberOfFiredChips(0) + mult->GetNumberOfFiredChips(1);
+ fFiredChips->Fill(firedChips, inputCount);
+ Printf("Accepted %d tracklets (%d fired chips)", inputCount, firedChips);
+
+ fTrackletsVsUnassigned->Fill(inputCount, mult->GetNumberOfSingleClusters());
}
-
- Int_t firedChips = mult->GetNumberOfFiredChips(0) + mult->GetNumberOfFiredChips(1);
- fFiredChips->Fill(firedChips, inputCount);
- Printf("Accepted %d tracklets (%d fired chips)", inputCount, firedChips);
-
- fTrackletsVsUnassigned->Fill(inputCount, mult->GetNumberOfSingleClusters());
}
else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
{
return;
}
+ Bool_t foundInEta10 = kFALSE;
+
if (vtxESD)
{
// get multiplicity from ESD tracks
if (eventTriggered && vtxESD)
fRawPt->Fill(pT);
+ if (esdTrack->Pt() < 0.15)
+ continue;
+
+ // INEL>0 trigger
+ if (TMath::Abs(esdTrack->Eta()) < 1 && esdTrack->Pt() > 0.15)
+ foundInEta10 = kTRUE;
+
if (fOnlyPrimaries)
{
if (label == 0)
++inputCount;
}
- if (inputCount > 30)
- Printf("Event with %d accepted TPC tracks. Period number: %d Orbit number: %x Bunch crossing number: %d", inputCount, fESD->GetPeriodNumber(), fESD->GetOrbitNumber(), fESD->GetBunchCrossNumber());
+ //if (inputCount > 30)
+ // Printf("Event with %d accepted TPC tracks. Period number: %d Orbit number: %x Bunch crossing number: %d", inputCount, fESD->GetPeriodNumber(), fESD->GetOrbitNumber(), fESD->GetBunchCrossNumber());
// TODO restrict inputCount used as measure for the multiplicity to |eta| < 1
delete list;
}
+
+ if (!foundInEta10)
+ eventTriggered = kFALSE;
}
else
return;
fMult->Fill(inputCount);
fdNdEtaAnalysisESD->FillTriggeredEvent(inputCount);
- fTriggerVsTime->SetPoint(fTriggerVsTime->GetN(), fESD->GetTimeStamp(), 1);
-
if (vtxESD)
{
// control hist
if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0)
fVertexVsMult->Fill(vtxESD->GetXv(), vtxESD->GetYv(), inputCount);
- fMultVtx->Fill(inputCount);
-
+ Int_t part05 = 0;
+ Int_t part10 = 0;
for (Int_t i=0; i<inputCount; ++i)
{
Float_t eta = etaArr[i];
else
fPartEta[2]->Fill(eta);
}
+
+ if (TMath::Abs(eta) < 0.5)
+ part05++;
+ if (TMath::Abs(eta) < 1.0)
+ part10++;
}
+
+ Int_t multAxis = inputCount;
+ if (fMultAxisEta1)
+ multAxis = part10;
+
+ fMultVtx->Fill(multAxis);
+ //if (TMath::Abs(vtx[2]) < 10)
+ // fMultVtx->Fill(part05);
// for event count per vertex
- fdNdEtaAnalysisESD->FillEvent(vtx[2], inputCount);
+ fdNdEtaAnalysisESD->FillEvent(vtx[2], multAxis);
// control hist
- if (inputCount > 0)
+ if (multAxis > 0)
fEvents->Fill(vtx[2]);
if (fReadMC)
thirdDim = particle->Phi();
}
else
- thirdDim = inputCount;
+ thirdDim = multAxis;
}
else
thirdDim = particle->Pt();
} // end of track loop
// for event count per vertex
- fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], inputCount);
+ fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], multAxis);
}
}
}
-
if (etaArr)
delete[] etaArr;
if (labelArr)
TArrayF vtxMC(3);
genHeader->PrimaryVertex(vtxMC);
-
+
// get process type
- Int_t processType = AliPWG0Helper::GetEventProcessType(header);
+ Int_t processType = AliPWG0Helper::GetEventProcessType(fESD, header, stack, fDiffTreatment);
AliDebug(AliLog::kDebug+1, Form("Found process type %d", processType));
- if (processType==AliPWG0Helper::kInvalidProcess)
+ if (processType == AliPWG0Helper::kInvalidProcess)
AliDebug(AliLog::kError, Form("Unknown process type %d.", processType));
// loop over mc particles
Int_t nPrim = stack->GetNprimary();
Int_t nAcceptedParticles = 0;
+ Bool_t oneParticleEvent = kFALSE;
// count particles first, then fill
for (Int_t iMc = 0; iMc < nPrim; ++iMc)
//AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
Float_t eta = particle->Eta();
- // make a rough eta cut (so that nAcceptedParticles is not too far off the true measured value (NB: this histograms are only gathered for comparison))
+ if (TMath::Abs(eta) < 1.0)
+ oneParticleEvent = kTRUE;
+
+ // make a rough eta cut (so that nAcceptedParticles is not too far off the true measured value (NB: these histograms are only gathered for comparison))
if (TMath::Abs(eta) < 1.5) // && pt > 0.3)
nAcceptedParticles++;
}
if (processType == AliPWG0Helper::kND)
fdNdEtaAnalysisND->FillTrack(vtxMC[2], eta, thirdDim);
+
+ if (oneParticleEvent)
+ fdNdEtaAnalysisOnePart->FillTrack(vtxMC[2], eta, thirdDim);
if (eventTriggered)
{
fdNdEtaAnalysisNSD->FillEvent(vtxMC[2], nAcceptedParticles);
if (processType == AliPWG0Helper::kND)
fdNdEtaAnalysisND->FillEvent(vtxMC[2], nAcceptedParticles);
+ if (oneParticleEvent)
+ fdNdEtaAnalysisOnePart->FillEvent(vtxMC[2], nAcceptedParticles);
if (eventTriggered)
{
for (Int_t i=0; i<3; ++i)
fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
- fVertexResolution = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
+ fVertexResolution = dynamic_cast<TH2F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("vertex_check"));
fVertexVsMult = dynamic_cast<TH3F*> (fOutput->FindObject("fVertexVsMult"));
fFiredChips = dynamic_cast<TH2F*> (fOutput->FindObject("fFiredChips"));
fTrackletsVsClusters = dynamic_cast<TH2F*> (fOutput->FindObject("fTrackletsVsClusters"));
fTrackletsVsUnassigned = dynamic_cast<TH2F*> (fOutput->FindObject("fTrackletsVsUnassigned"));
- fTriggerVsTime = dynamic_cast<TGraph*> (fOutput->FindObject("fTriggerVsTime"));
fStats = dynamic_cast<TH1F*> (fOutput->FindObject("fStats"));
fStats2 = dynamic_cast<TH2F*> (fOutput->FindObject("fStats2"));
fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCuts"));
- fPhysicsSelection = dynamic_cast<AliPhysicsSelection*> (fOutput->FindObject("AliPhysicsSelection_outputlist"));
- fTriggerAnalysis = dynamic_cast<AliTriggerAnalysis*> (fOutput->FindObject("AliTriggerAnalysis"));
}
if (!fdNdEtaAnalysisESD)
if (fPartEta[0])
{
- Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-19.9), fEvents->GetXaxis()->FindBin(-0.001));
- Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(19.9));
+ Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-9.9), fEvents->GetXaxis()->FindBin(-0.001));
+ Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(9.9));
Printf("%d events with vertex used", events1 + events2);
+ Printf("%d total events with vertex, %d total triggered events, %d triggered events with 0 multiplicity", (Int_t) fMultVtx->Integral(), (Int_t) fMult->Integral(), (Int_t) fMult->GetBinContent(1));
if (events1 > 0 && events2 > 0)
{
if (fEsdTrackCuts)
fEsdTrackCuts->SaveHistograms("esd_track_cuts");
- if (fPhysicsSelection)
- {
- fPhysicsSelection->SaveHistograms("physics_selection");
- fPhysicsSelection->Print();
- }
-
- if (fTriggerAnalysis)
- fTriggerAnalysis->SaveHistograms();
-
if (fMult)
fMult->Write();
fEvents->Write();
if (fVertexResolution)
+ {
fVertexResolution->Write();
+ fVertexResolution->ProjectionX();
+ fVertexResolution->ProjectionY();
+ }
if (fDeltaPhi)
fDeltaPhi->Write();
if (fTrackletsVsUnassigned)
fTrackletsVsUnassigned->Write();
- if (fTriggerVsTime)
- fTriggerVsTime->Write();
-
if (fStats)
fStats->Write();
fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
fdNdEtaAnalysisND = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaND"));
fdNdEtaAnalysisNSD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaNSD"));
+ fdNdEtaAnalysisOnePart = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaOnePart"));
fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
fdNdEtaAnalysisND->Finish(0, -1, AlidNdEtaCorrection::kNone);
fdNdEtaAnalysisNSD->Finish(0, -1, AlidNdEtaCorrection::kNone);
+ fdNdEtaAnalysisOnePart->Finish(0, -1, AlidNdEtaCorrection::kNone);
fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);
fdNdEtaAnalysis->SaveHistograms();
fdNdEtaAnalysisND->SaveHistograms();
fdNdEtaAnalysisNSD->SaveHistograms();
+ fdNdEtaAnalysisOnePart->SaveHistograms();
fdNdEtaAnalysisTr->SaveHistograms();
fdNdEtaAnalysisTrVtx->SaveHistograms();
fdNdEtaAnalysisTracks->SaveHistograms();
class TH2F;
class TH3F;
class AliESDEvent;
-class TGraph;
-class AliPhysicsSelection;
class AliTriggerAnalysis;
class AlidNdEtaTask : public AliAnalysisTask {
void SetDeltaPhiCut(Float_t cut) { fDeltaPhiCut = cut; }
void SetCheckEventType(Bool_t flag = kTRUE) { fCheckEventType = flag; }
void SetSymmetrize(Bool_t flag = kTRUE) { fSymmetrize = flag; }
+ void SetMultAxisEta1(Bool_t flag = kTRUE) { fMultAxisEta1 = flag; }
+ void SetDiffTreatment(AliPWG0Helper::DiffTreatment diffTreatment) { fDiffTreatment = diffTreatment; }
void SetOption(const char* opt) { fOption = opt; }
Bool_t fUseMCKine; // use the MC values for each found track/tracklet (for syst. check)
Bool_t fCheckEventType; // check if event type is physics (for real data)
Bool_t fSymmetrize; // move all negative to positive eta
+ Bool_t fMultAxisEta1; // restrict multiplicity count to |eta| < 1
+ AliPWG0Helper::DiffTreatment fDiffTreatment; // how to identify SD events (see AliPWG0Helper::GetEventProcessType)
AliESDtrackCuts* fEsdTrackCuts; // Object containing the parameters of the esd track cuts
- AliPhysicsSelection* fPhysicsSelection; // Event Selection object
AliTriggerAnalysis* fTriggerAnalysis;
// Gathered from ESD
dNdEtaAnalysis* fdNdEtaAnalysisESD; //! contains the dndeta from the ESD
// control hists
- TH1F* fMult; //! raw multiplicity histogram (control histogram)
+ TH1F* fMult; //! raw multiplicity histogram
TH1F* fMultVtx; //! raw multiplicity histogram of evts with vtx (control histogram)
TH1F* fPartEta[3]; //! counted particles as function of eta (full vertex range, below 0 range, above 0 range)
TH1F* fEvents; //! events counted as function of vtx
- TH1F* fVertexResolution; //! z resolution of the vertex
+ TH2F* fVertexResolution; //! z resolution of the vertex
// Gathered from MC (when fReadMC is set)
dNdEtaAnalysis* fdNdEtaAnalysis; //! contains the dndeta from the full sample
dNdEtaAnalysis* fdNdEtaAnalysisND; //! contains the dndeta for the ND sample
dNdEtaAnalysis* fdNdEtaAnalysisNSD; //! contains the dndeta for the NSD sample
+ dNdEtaAnalysis* fdNdEtaAnalysisOnePart; //! contains the dndeta for the one particle sample
dNdEtaAnalysis* fdNdEtaAnalysisTr; //! contains the dndeta from the triggered events
dNdEtaAnalysis* fdNdEtaAnalysisTrVtx; //! contains the dndeta from the triggered events with vertex
dNdEtaAnalysis* fdNdEtaAnalysisTracks; //! contains the dndeta from the triggered events with vertex counted from the mc particles associated to the tracks (comparing this to the raw values from the esd shows the effect of the detector resolution)
TH2F* fFiredChips; //! fired chips l1+l2 vs. number of tracklets (only for SPD analysis)
TH2F* fTrackletsVsClusters; //! number of tracklets vs. clusters in all ITS detectors (only for SPD analysis)
TH2F* fTrackletsVsUnassigned; //! number of tracklets vs. number of unassigned clusters in L1 (only for SPD analysis)
- TGraph* fTriggerVsTime; //! trigger as function of event time
TH1F* fStats; //! further statistics : bin 1 = vertexer 3d, bin 2 = vertexer z, etc (see CreateOutputObjects)
TH2F* fStats2; //! V0 vs SPD statistics
{
loadlibs();
- AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
- if (!TFile::Open(correctionMapFile))
- return;
- dNdEtaCorrection->LoadHistograms();
-
+ AlidNdEtaCorrection* dNdEtaCorrection = 0;
+
+ if (correctionMapFile)
+ {
+ dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
+ if (!TFile::Open(correctionMapFile))
+ return;
+ dNdEtaCorrection->LoadHistograms();
+ }
+
TFile* file = TFile::Open(dataInput);
if (!file)
cout << "Error. File not found" << endl;
return;
}
-
- // Note: the last parameter does not define which analysis is going to happen, the histograms will be overwritten when loading from the f
+
+ Int_t backgroundEvents = 0;
+
+ //backgroundEvents = 1162+434; // Michele for MB1, run 104892, 15.02.10
+ //backgroundEvents = 842; // JF estimate for V0 systematic case 1
+ backgroundEvents = 6; // Michele for V0AND, run 104892, 15.02.10
+
+ //backgroundEvents = 1758+434; // MB1, run 104867-92
+
+ //backgroundEvents = 4398+961; // Michele for MB1, run 104824-52, 16.02.10
+ //backgroundEvents = 19; // Michele for V0AND, run 104824-52, 16.02.10
+
+ //backgroundEvents = -1; // use 0 bin from MC! for 2.36 TeV
+ //backgroundEvents = 900; // my estimate for 2.36 TeV
+
+ //backgroundEvents = 918; // V0OR for run 114786 w/o bunch intensities w/o proper 0 checking!
+ //backgroundEvents = 723; // V0OR for run 114798 w/o bunch intensities w/o proper 0 checking!
+
+ Printf("Subtracting %d background events!!!", backgroundEvents);
+ gSystem->Sleep(1000);
+
+ TH1* combinatoricsCorrection = 0;
+ if (1)
+ {
+ TFile::Open("corrComb.root");
+ combinatoricsCorrection = (TH1*) gFile->Get("ratioofratios");
+ file->cd();
+ }
+
+ // Note: the last parameter does not define which analysis is going to happen, the histograms will be overwritten when loading from the file
dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD");
fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
- fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.21, AlidNdEtaCorrection::kNSD, "ESD -> NSD");
+ fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.151, AlidNdEtaCorrection::kNSD, "ESD -> NSD", backgroundEvents, combinatoricsCorrection);
//fdNdEtaAnalysis->DrawHistograms(kTRUE);
TFile* file2 = TFile::Open(dataOutput, "RECREATE");
fdNdEtaAnalysis->SaveHistograms();
-
+
file->cd();
dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
- fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.21, AlidNdEtaCorrection::kINEL, "ESD -> full inelastic");
+ fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.151, AlidNdEtaCorrection::kINEL, "ESD -> full inelastic", backgroundEvents, combinatoricsCorrection);
//fdNdEtaAnalysis->DrawHistograms(kTRUE);
file2->cd();
fdNdEtaAnalysis->SaveHistograms();
file->cd();
fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
- fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.21, AlidNdEtaCorrection::kVertexReco, "ESD -> minimum bias");
+ fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.151, AlidNdEtaCorrection::kVertexReco, "ESD -> minimum bias", backgroundEvents, combinatoricsCorrection);
//fdNdEtaAnalysis->DrawHistograms(kTRUE);
file2->cd();
fdNdEtaAnalysis->SaveHistograms();
file->cd();
fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
- fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.21, AlidNdEtaCorrection::kTrack2Particle, "ESD -> MB with vertex");
+ fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.151, AlidNdEtaCorrection::kTrack2Particle, "ESD -> MB with vertex", backgroundEvents, combinatoricsCorrection);
//fdNdEtaAnalysis->DrawHistograms(kTRUE);
file2->cd();
fdNdEtaAnalysis->SaveHistograms();
file->cd();
fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
- fdNdEtaAnalysis->Finish(0, 0.21, AlidNdEtaCorrection::kNone, "ESD raw with pt cut");
+ fdNdEtaAnalysis->Finish(0, 0.151, AlidNdEtaCorrection::kNone, "ESD raw with pt cut", backgroundEvents, combinatoricsCorrection);
//fdNdEtaAnalysis->DrawHistograms(kTRUE);
file2->cd();
fdNdEtaAnalysis->SaveHistograms();
file->cd();
fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracksAll", "dndetaTracksAll");
fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
- fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone, "ESD raw");
+ fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone, "ESD raw", backgroundEvents, combinatoricsCorrection);
+ //fdNdEtaAnalysis->DrawHistograms(kTRUE);
+ file2->cd();
+ fdNdEtaAnalysis->SaveHistograms();
+
+ file->cd();
+ fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaOnePart", "dndetaOnePart");
+ fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
+ fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.151, AlidNdEtaCorrection::kOnePart, "ESD -> OnePart", backgroundEvents, combinatoricsCorrection);
//fdNdEtaAnalysis->DrawHistograms(kTRUE);
file2->cd();
fdNdEtaAnalysis->SaveHistograms();
TFile::Open(correctionMapFile);
dNdEtaCorrection->LoadHistograms();
- fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.21, AlidNdEtaCorrection::kINEL);
+ fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.151, AlidNdEtaCorrection::kINEL);
//fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0, AlidNdEtaCorrection::kINEL);
//fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0, AlidNdEtaCorrection::kTrack2Particle);
}
else
- fdNdEtaAnalysis->Finish(0, 0.21, AlidNdEtaCorrection::kNone);
+ fdNdEtaAnalysis->Finish(0, 0.151, AlidNdEtaCorrection::kNone);
+ return fdNdEtaAnalysis;
+
fdNdEtaAnalysis->DrawHistograms(simple);
TH1* hist = fdNdEtaAnalysis->GetdNdEtaHistogram(1);
return fdNdEtaAnalysis;
}
-void correct(Bool_t onlyESD = kFALSE)
+void FinishMC()
{
- FinishAnalysisAll();
+ loadlibs();
+
+ result1 = (dNdEtaAnalysis*) FinishAnalysis("analysis_mc.root", "dndeta", 0, 0, 1);
+ result2 = (dNdEtaAnalysis*) FinishAnalysis("analysis_mc.root", "dndetaNSD", 0, 0, 1);
+
+ file = TFile::Open("out.root", "RECREATE");
+ result1->SaveHistograms();
+ result2->SaveHistograms();
+ file->Close();
+}
+
+void correct(Bool_t onlyESD = kFALSE, Bool_t mergedXSections = kTRUE)
+{
+ gSystem->Unlink("analysis_esd.root");
+ if (mergedXSections)
+ FinishAnalysisAll("analysis_esd_raw.root", "analysis_esd.root", "correction_map2.root", "dndeta_correction_ua5");
+ else
+ FinishAnalysisAll("analysis_esd_raw.root", "analysis_esd.root", "correction_map.root", "dndeta_correction");
+
gROOT->ProcessLine(".L $ALICE_ROOT/PWG0/dNdEta/drawPlots.C");
dNdEta(onlyESD);
}
}
//____________________________________________________________________
-void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType, const char* tag)
+void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType, const char* tag, Int_t backgroundSubtraction, TH1* combinatoricsCorrection)
{
//
// correct with the given correction values and calculate dNdEta and pT distribution
fTag.Form("Correcting dN/deta spectrum (data: %d) >>> %s <<<. Correction type: %d, pt cut: %.2f.", (Int_t) fAnalysisMode, tag, (Int_t) correctionType, ptCut);
Printf("\n\n%s", fTag.Data());
+
+ if (combinatoricsCorrection)
+ Printf("Combinatorics subtraction active!");
// set corrections to 1
fData->SetCorrectionToUnity();
eventCorr->Multiply(correction->GetTriggerBiasCorrectionND()->GetEventCorrection()->GetCorrectionHistogram());
break;
}
+ case AlidNdEtaCorrection::kOnePart :
+ {
+ trackCorr->Multiply(correction->GetTriggerBiasCorrectionOnePart()->GetTrackCorrection()->GetCorrectionHistogram());
+ eventCorr->Multiply(correction->GetTriggerBiasCorrectionOnePart()->GetEventCorrection()->GetCorrectionHistogram());
+ break;
+ }
default : break;
}
}
fData->ResetErrorsOnCorrections();
fData->Multiply();
-
+
if (correctionType >= AlidNdEtaCorrection::kVertexReco)
{
// There are no events with vertex that have 0 multiplicity, therefore
//TH2* measuredEvents = fData->GetEventCorrection()->GetMeasuredHistogram();
TH2* correctedEvents = fData->GetEventCorrection()->GetGeneratedHistogram();
+ TH2* eAll = correction->GetCorrection(correctionType)->GetEventCorrection()->GetGeneratedHistogram();
TH2* eTrig = correction->GetVertexRecoCorrection()->GetEventCorrection()->GetGeneratedHistogram();
TH2* eTrigVtx = correction->GetVertexRecoCorrection()->GetEventCorrection()->GetMeasuredHistogram();
- TH1* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 2, rawMeasured->GetNbinsY()+1);
+ TH1* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 2, eTrigVtx->GetNbinsY()+1);
//new TCanvas; correctedEvents->DrawCopy("TEXT");
Int_t triggeredEventsWith0Mult = (Int_t) fMult->GetBinContent(1);
Printf("%d triggered events with 0 mult. -- %d triggered events with vertex", triggeredEventsWith0Mult, allEventsWithVertex);
+
+ if (backgroundSubtraction > 0)
+ {
+ triggeredEventsWith0Mult -= backgroundSubtraction;
+ Printf("Subtracted %d background events from 0 mult. bin", backgroundSubtraction);
+ }
TH1* kineBias = (TH1*) vertexDist->Clone("kineBias");
kineBias->Reset();
-
+
// loop over vertex bins
for (Int_t i = 1; i <= rawMeasured->GetNbinsX(); i++)
{
if (eTrigVtx_projx->GetBinContent(i) == 0)
continue;
- Double_t fZ = eTrigVtx_projx->Integral(0, eTrigVtx_projx->GetNbinsX()+1) / eTrigVtx_projx->GetBinContent(i) *
- eTrig->GetBinContent(i, 1) / eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1);
- kineBias->SetBinContent(i, fZ);
-
- events *= fZ;
-
- // multiply with trigger correction if set above
- events *= fData->GetEventCorrection()->GetCorrectionHistogram()->GetBinContent(i, 1);
-
// calculate how many events we would have got with a pure MC-based correction
// in the given bin: measured events with vertex (mult > 0) * triggered events with mult 0 (mc) / events with vertex and mult > 0 (mc) * trigger correction for bin 0
- Double_t mcEvents = vertexDist->GetBinContent(i) * eTrig->GetBinContent(i, 1) / eTrigVtx_projx->GetBinContent(i) * fData->GetEventCorrection()->GetCorrectionHistogram()->GetBinContent(i, 1);
+ Printf("+++ 0-Bin Correction for bin %d +++", i);
+ Printf(" Events: %f", vertexDist->GetBinContent(i));
+ Printf(" Ratio triggered N==0 / triggered vertex N>0: %f", eTrig->GetBinContent(i, 1) / eTrigVtx_projx->GetBinContent(i));
+ Printf(" Ratio all N==0 / triggered vertex N>0: %f", eAll->GetBinContent(i, 1) / eTrigVtx_projx->GetBinContent(i));
+ Printf(" Correction factor: %f", fData->GetEventCorrection()->GetCorrectionHistogram()->GetBinContent(i, 1));
+
+ //Double_t mcEvents = vertexDist->GetBinContent(i) * eTrig->GetBinContent(i, 1) / eTrigVtx_projx->GetBinContent(i) * fData->GetEventCorrection()->GetCorrectionHistogram()->GetBinContent(i, 1);
+ Double_t mcEvents = vertexDist->GetBinContent(i) * eAll->GetBinContent(i, 1) / eTrigVtx_projx->GetBinContent(i);
+ if (backgroundSubtraction == -1)
+ {
+ Printf("Using MC value for 0-bin correction!");
+ events = mcEvents;
+ }
+ else
+ {
+ Double_t fZ = eTrigVtx_projx->Integral(0, eTrigVtx_projx->GetNbinsX()+1) / eTrigVtx_projx->GetBinContent(i) *
+ eTrig->GetBinContent(i, 1) / eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1);
+
+ kineBias->SetBinContent(i, fZ);
- Printf("Bin %d, alpha is %.2f, fZ is %.3f, number of events with 0 mult.: %.2f (MC comparison: %.2f)", i, alpha * 100., fZ, events, mcEvents);
- Printf("Using MC value for 0-bin correction!");
- events = mcEvents;
+ events *= fZ;
+
+ // multiply with trigger correction if set above
+ events *= fData->GetEventCorrection()->GetCorrectionHistogram()->GetBinContent(i, 1);
+
+ Printf(" Bin %d, alpha is %.2f%%, fZ is %.3f, number of events with 0 mult.: %.2f (MC comparison: %.2f)", i, alpha * 100., fZ, events, mcEvents);
+ }
correctedEvents->SetBinContent(i, 1, events);
}
}
//new TCanvas(tag, tag, 800, 600); vtxVsEta->DrawCopy("COLZ");
-
+
// clear result histograms
for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
{
Int_t vertexBinEnd = vertexHist->GetXaxis()->FindBin(vertexRangeEnd[vertexPos]);
const Int_t *binBegin = 0;
- const Int_t maxBins = 14;
+ const Int_t maxBins = 40;
if (dataHist->GetNbinsY() != maxBins)
AliFatal(Form("Binning of acceptance is different from data histogram: data=%d, acceptance=%d", dataHist->GetNbinsY(), maxBins));
// produce with drawPlots.C: DetermineAcceptance(...)
if (fAnalysisMode & AliPWG0Helper::kSPD)
{
- const Int_t binBeginSPD[maxBins] = {-1, 16, 13, 9, 7, 5, 4, 4, 3, 3, 2, 2, 2, -1};
-
+ //const Int_t binBeginSPD[maxBins] = {15, 14, 13, 12, 11, 10, 9, 9, 8, 7, 7, 6, 6, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4};
+ const Int_t binBeginSPD[maxBins] = {19, 18, 17, 15, 14, 12, 10, 9, 8, 7, 6, 6, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4};
+
binBegin = binBeginSPD;
}
else if (fAnalysisMode & AliPWG0Helper::kTPC)
{
- //const Int_t binBeginTPC[30] = {-1, -1, -1, -1, -1, -1, -1, -1, -1, 4, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1}; // limit 5, pt cut off 0.2 mev/c
- //const Int_t binBeginTPC[maxBins] = {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 9, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}; // limit 5
+ const Int_t binBeginTPC[maxBins] = {-1, -1, -1, -1, -1, -1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1};
- //binBegin = binBeginTPC;
+ binBegin = binBeginTPC;
}
else if (fAnalysisMode & AliPWG0Helper::kTPCITS)
{
- // TODO create map
+ const Int_t binBeginTPCITS[maxBins] = {-1, -1, -1, -1, -1, -1, -1, 14, 10, 8, 7, 6, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1, -1, -1};
+
+ binBegin = binBeginTPCITS;
}
Int_t vtxBegin = 1;
Float_t dndeta = sum / totalEvents;
Float_t error = TMath::Sqrt(sumError2) / totalEvents;
+ // correct for additional combinatorics
+ Float_t combCorr = 0;
+ if (combinatoricsCorrection)
+ {
+ combCorr = combinatoricsCorrection->GetBinContent(combinatoricsCorrection->GetXaxis()->FindBin(vtxVsEta->GetYaxis()->GetBinCenter(iEta)));
+ dndeta *= combCorr;
+ error *= combCorr;
+ }
+
dndeta = dndeta / fdNdEta[vertexPos]->GetBinWidth(bin);
error = error / fdNdEta[vertexPos]->GetBinWidth(bin);
fdNdEtaPtCutOffCorrected[vertexPos]->SetBinContent(bin, dndeta);
fdNdEtaPtCutOffCorrected[vertexPos]->SetBinError(bin, error);
- //Printf("Bin %d has dN/deta = %f +- %f; %.2f tracks %.2f events (outside acceptance: %.2f tracks, %.2f events)", bin, dndeta, error, sum, totalEvents, unusedTracks, unusedEvents);
+ //Printf("Bin %d has dN/deta = %f +- %f; %.2f tracks %.2f events (outside acceptance: %.2f tracks, %.2f events) (comb. corr: %f)", bin, dndeta, error, sum, totalEvents, unusedTracks, unusedEvents, combCorr);
}
}
}
void FillEvent(Float_t vtx, Float_t n);
void FillTriggeredEvent(Float_t n);
- void Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType, const char* tag = "");
+ void Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType, const char* tag = "", Int_t backgroundSubtraction = 0, TH1* combinatoricsCorrection = 0);
void DrawHistograms(Bool_t simple = kFALSE);
void LoadHistograms(const Char_t* dir = 0);
fdNdEtaAnalysis->LoadHistograms("dndeta");
TH1* histESD = (TH1*) file->Get("dndeta/dNdEta_corrected");
+ TH1* histESD1 = (TH1*) file->Get("dndeta/dNdEta_corrected_1");
+ TH1* histESD2 = (TH1*) file->Get("dndeta/dNdEta_corrected_2");
TH1* histESDnsd = (TH1*) file->Get("dndetaNSD/dNdEta_corrected");
TH1* histESDnsdNoPt = (TH1*) file->Get("dndetaNSD/dNdEta");
+ TH1* histESDonePart = (TH1*) file->Get("dndetaOnePart/dNdEta_corrected");
TH1* histESDNoPt = (TH1*) file->Get("dndeta/dNdEta");
TH1* histESDMB = (TH1*) file->Get("dndetaTr/dNdEta_corrected");
TH1* histESDMBNoPt = (TH1*) file->Get("dndetaTr/dNdEta");
TH1* histESDMBTracksNoPt = (TH1*) file->Get("dndetaTracks/dNdEta");
Prepare1DPlot(histESD);
+ Prepare1DPlot(histESD1);
+ Prepare1DPlot(histESD2);
Prepare1DPlot(histESDnsd);
+ Prepare1DPlot(histESDonePart);
Prepare1DPlot(histESDMB);
Prepare1DPlot(histESDMBVtx);
histESD->SetLineWidth(0);
histESDnsd->SetLineWidth(0);
+ histESDonePart->SetLineWidth(0);
histESDMB->SetLineWidth(0);
histESDMBVtx->SetLineWidth(0);
histESD->SetMarkerColor(1);
histESDnsd->SetMarkerColor(6);
+ histESDonePart->SetMarkerColor(3);
histESDMB->SetMarkerColor(2);
histESDMBVtx->SetMarkerColor(4);
histESD->SetLineColor(1);
histESDnsd->SetLineColor(6);
+ histESDonePart->SetLineColor(3);
histESDMB->SetLineColor(2);
histESDMBVtx->SetLineColor(4);
histESD->SetMarkerStyle(20);
histESDnsd->SetMarkerStyle(29);
+ histESDonePart->SetMarkerStyle(24);
histESDMB->SetMarkerStyle(21);
histESDMBVtx->SetMarkerStyle(22);
histESDMBVtxNoPt->SetMarkerStyle(22);
histESDMBTracksNoPt->SetMarkerStyle(23);
- Float_t etaLimit = (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC) ? 0.89 : 1.79;
+ Float_t etaLimit = (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC) ? 0.89 : 1.99;
Float_t etaPlotLimit = (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC) ? 1.2 : 2.3;
//Float_t etaLimit = (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC) ? 0.89 : 1.39;
//Float_t etaPlotLimit = (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC) ? 1.2 : 1.9;
histESDMB->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
histESD->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
histESDnsd->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
+ histESDonePart->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
histESDNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
histESDMBNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
legend->AddEntry(histESDMB, "Triggered");
legend->AddEntry(histESD, "All INEL events");
legend->AddEntry(histESDnsd, "All NSD events");
+ legend->AddEntry(histESDonePart, "One Particle");
TH2F* dummy = new TH2F("dummy", "", 100, -etaPlotLimit, etaPlotLimit, 1000, 0, max * 1.1);
dummy->GetYaxis()->SetRangeUser(2.1, max * 1.1);
histESDMB->Draw("SAME");
histESD->Draw("SAME");
histESDnsd->Draw("SAME");
+ histESDonePart->Draw("SAME");
legend->Draw();
-
+
if (save)
{
- canvas->SaveAs("dNdEta1.gif");
- canvas->SaveAs("dNdEta1.eps");
+ canvas->SaveAs("dNdEta1.png");
+ //canvas->SaveAs("dNdEta1.eps");
}
- histESD->Fit("pol0", "0", "", -0.45, 0.45);
- histESDnsd->Fit("pol0", "0", "", -0.45, 0.45);
-
+ histESD->Fit("pol0", "0", "", -0.49, 0.49);
+ histESDnsd->Fit("pol0", "0", "", -0.49, 0.49);
+ histESDonePart->Fit("pol0", "0", "", -0.49, 0.49);
+ histESDonePart->Fit("pol0", "0", "", -0.99, 0.99);
+
+ canvas = new TCanvas("dNdEta1_mirrored", "dNdEta1_mirrored", 500, 500);
+ canvas->SetGridx();
+ canvas->SetGridy();
+
+ dummy->DrawCopy()->GetXaxis()->SetRangeUser(0, 100);
+ histESD->DrawCopy("SAME")->SetMarkerStyle(24);
+ histESDnsd->DrawCopy("SAME")->SetMarkerStyle(24);
+
+ graph = new TGraphErrors(histESD);
+ for (Int_t i=0; i<graph->GetN(); i++)
+ graph->GetX()[i] *= -1;
+ graph->SetMarkerStyle(5);
+ graph->Draw("P SAME");
+
+ graph = new TGraphErrors(histESDnsd);
+ for (Int_t i=0; i<graph->GetN(); i++)
+ graph->GetX()[i] *= -1;
+ graph->SetMarkerStyle(5);
+ graph->SetMarkerColor(histESDnsd->GetMarkerColor());
+ graph->Draw("P SAME");
+
+ canvas = new TCanvas("dNdEta1_ratio", "dNdEta1_ratio", 500, 500);
+ canvas->SetGridx();
+ canvas->SetGridy();
+
+ dummy_clone = dummy->DrawCopy();
+ dummy_clone->GetXaxis()->SetRangeUser(0, 100);
+ dummy_clone->GetYaxis()->SetRangeUser(0.5, 1.5);
+
+ graph = new TGraphErrors(histESD);
+ for (Int_t i=0; i<graph->GetN(); i++)
+ {
+ Int_t bin = histESD->GetXaxis()->FindBin(-graph->GetX()[i]);
+ if (histESD->GetBinContent(bin) > 0 && graph->GetY()[i] > 0)
+ {
+ graph->GetEY()[i] = sqrt(graph->GetEY()[i] * graph->GetEY()[i] / graph->GetY()[i] / graph->GetY()[i]
+ + histESD->GetBinError(bin) * histESD->GetBinError(bin) / histESD->GetBinContent(bin) / histESD->GetBinContent(bin));
+ graph->GetY()[i] /= histESD->GetBinContent(bin);
+ graph->GetEY()[i] *= graph->GetY()[i];
+ }
+ else
+ graph->GetY()[i] = 0;
+ }
+ graph->SetMarkerStyle(5);
+ graph->Draw("P SAME");
+
+ graph = new TGraphErrors(histESDnsd);
+ for (Int_t i=0; i<graph->GetN(); i++)
+ {
+ Int_t bin = histESDnsd->GetXaxis()->FindBin(-graph->GetX()[i]);
+ if (histESDnsd->GetBinContent(bin) > 0 && graph->GetY()[i] > 0)
+ {
+ graph->GetEY()[i] = sqrt(graph->GetEY()[i] * graph->GetEY()[i] / graph->GetY()[i] / graph->GetY()[i]
+ + histESDnsd->GetBinError(bin) * histESDnsd->GetBinError(bin) / histESDnsd->GetBinContent(bin) / histESDnsd->GetBinContent(bin));
+ graph->GetY()[i] /= histESDnsd->GetBinContent(bin);
+ graph->GetEY()[i] *= graph->GetY()[i];
+ graph->GetY()[i] += 0.2;
+ }
+ }
+ graph->SetMarkerStyle(5);
+ graph->SetMarkerColor(histESDnsd->GetMarkerColor());
+ graph->Draw("P SAME");
+
+ canvas = new TCanvas("dNdEta1_vertex", "dNdEta1_vertex", 500, 500);
+ dummy->DrawCopy();
+ histESD->DrawCopy("SAME");
+ histESD1->SetLineColor(2);
+ histESD1->DrawCopy("SAME");
+ histESD2->SetLineColor(4);
+ histESD2->DrawCopy("SAME");
+
if (onlyESD)
return;
TFile* file2 = TFile::Open("analysis_mc.root");
- TH1* histMCTrVtx = (TH1*) GetMCHist("dndetaTrVtx", -1, "MC: MB with trigger")->Clone("histMCTrVtx");
+ TH1* histMCTrVtx = (TH1*) GetMCHist("dndetaTrVtx", -1, "MC: MB with vertex")->Clone("histMCTrVtx");
TH1* ratioTrVtx = (TH1*) DrawdNdEtaRatio(histESDMBVtx, histMCTrVtx, "triggered_vertex", etaPlotLimit)->Clone();
TH1* histMC = (TH1*) GetMCHist("dndeta", -1, "MC: full inelastic")->Clone("histMC");
TH1* histMCTr = (TH1*) GetMCHist("dndetaTr", -1, "MC: minimum bias")->Clone("histMCTr");
TH1* histMCnsd = (TH1*) GetMCHist("dndetaNSD", -1, "MC: NSD")->Clone("histMCnsd");
+ TH1* histMConePart = (TH1*) GetMCHist("dndetaOnePart", -1, "MC: OnePart")->Clone("histMConePart");
- TH1* histMCPtCut = (TH1*) GetMCHist("dndeta", 0.21, "MC: full inelastic, pt cut")->Clone("histMCPtCut");
- TH1* histMCTrPtCut = (TH1*) GetMCHist("dndetaTr", 0.21, "MC: minimum bias, pt cut")->Clone("histMCTrPtCut");
- TH1* histMCTrVtxPtCut = (TH1*) GetMCHist("dndetaTrVtx", 0.21, "MC: MB with trigger, pt cut")->Clone("histMCTrVtxPtCut");
- TH1* histMCnsdNoPt = (TH1*) GetMCHist("dndetaNSD", 0.21, "MC: NSD, put cut")->Clone("histMCnsdNoPt");
- TH1* histMCTracksPtCut = (TH1*) GetMCHist("dndetaTracks", 0.21, "MC: Tracks w/o resolution effect, pt cut")->Clone("histMCTracksPtCut");
+ TH1* histMCPtCut = (TH1*) GetMCHist("dndeta", 0.151, "MC: full inelastic, pt cut")->Clone("histMCPtCut");
+ TH1* histMCTrPtCut = (TH1*) GetMCHist("dndetaTr", 0.151, "MC: minimum bias, pt cut")->Clone("histMCTrPtCut");
+ TH1* histMCTrVtxPtCut = (TH1*) GetMCHist("dndetaTrVtx", 0.151, "MC: MB with vertex, pt cut")->Clone("histMCTrVtxPtCut");
+ TH1* histMCnsdNoPt = (TH1*) GetMCHist("dndetaNSD", 0.151, "MC: NSD, put cut")->Clone("histMCnsdNoPt");
+ TH1* histMCTracksPtCut = (TH1*) GetMCHist("dndetaTracks", 0.151, "MC: Tracks w/o resolution effect, pt cut")->Clone("histMCTracksPtCut");
Prepare1DPlot(histMC);
Prepare1DPlot(histMCnsd);
TH1* ratioTrVtx = (TH1*) DrawdNdEtaRatio(histESDMBVtx, histMCTrVtx, "triggered_vertex", etaPlotLimit)->Clone();
TH1* ratioTrVtxNoPt = (TH1*) DrawdNdEtaRatio(histESDMBVtxNoPt, histMCTrVtxPtCut, "triggered_vertex_nopt", etaPlotLimit)->Clone();
TH1* ratioNSD = (TH1*) DrawdNdEtaRatio(histESDnsd, histMCnsd, "NSD", etaPlotLimit)->Clone();
+ TH1* ratioOnePart = (TH1*) DrawdNdEtaRatio(histESDonePart, histMConePart, "OnePart", etaPlotLimit)->Clone();
// draw ratios of single steps
c7 = new TCanvas("all_ratios", "all_ratios", 600, 600);
ratioNSD->Draw("HIST EP SAME");
legend7->Draw();
- c7->SaveAs("ratios.eps");
+ //c7->SaveAs("ratios.eps");
new TCanvas;
dummy2->DrawCopy();
PrintIntegratedDeviation(histMC, histESD, "all events");
PrintIntegratedDeviation(histMCnsd, histESDnsd, "all events (NSD)");
+ PrintIntegratedDeviation(histMConePart, histESDonePart, "all events (INEL>0)");
PrintIntegratedDeviation(histMCTr, histESDMB, "triggered");
PrintIntegratedDeviation(histMCTrVtx, histESDMBVtx, "trigger, vertex");
PrintIntegratedDeviation(histMCPtCut, histESDNoPt, "all events (no pt corr)");
legend->Draw();
}
+void CompareTwodNdEta(const char* fileName1, const char* fileName2, Bool_t errorsCorrelated = kFALSE)
+{
+ c = new TCanvas;
+
+ c->SetGridx();
+ c->SetGridy();
+
+ hist = new TH2F("dummy", ";#eta;dN_{ch}/d#eta", 100, -2.5, 2.5, 100, 0, 8);
+ hist->SetStats(0);
+ hist->DrawCopy();//->GetYaxis()->SetRangeUser(2, 4.5);
+
+ l = new TLegend(0.2, 0.13, 0.8, 0.35);
+ l->SetNColumns(2);
+ l->SetFillColor(0);
+
+ TH1* histESD[2];
+ TH1* histESDnsd[2];
+
+ for (Int_t i=0; i<2; i++)
+ {
+ if (i == 0)
+ file = TFile::Open(fileName1);
+ if (i == 1)
+ {
+ if (fileName2 == 0)
+ break;
+ file = TFile::Open(fileName2);
+ }
+
+ histESD[i] = (TH1*) file->Get("dndeta/dNdEta_corrected");
+ histESDnsd[i] = (TH1*) file->Get("dndetaNSD/dNdEta_corrected");
+
+ histESD[i]->SetMarkerStyle(20 + i*4);
+ histESDnsd[i]->SetMarkerStyle(21 + i*4);
+
+ histESD[i]->SetMarkerColor(i+1);
+ histESD[i]->SetLineColor(i+1);
+ histESDnsd[i]->SetMarkerColor(i+1);
+ histESDnsd[i]->SetLineColor(i+1);
+
+ histESD[i]->DrawCopy("SAME");
+ histESDnsd[i]->DrawCopy("SAME");
+
+ l->AddEntry(histESD[i], Form("Data %d INEL", i), "P");
+ l->AddEntry(histESDnsd[i], Form("Data %d NSD", i), "P");
+ }
+
+ if (0)
+ {
+ TGraphErrors *gre = new TGraphErrors(16);
+ gre->SetFillColor(4);
+ gre->SetMarkerColor(4);
+ gre->SetMarkerStyle(26);
+ gre->SetPoint(0,0.125,3.14);
+ gre->SetPointError(0,0,0.07);
+ gre->SetPoint(1,0.375,3.04);
+ gre->SetPointError(1,0,0.07);
+ gre->SetPoint(2,0.625,3.17);
+ gre->SetPointError(2,0,0.07);
+ gre->SetPoint(3,0.875,3.33);
+ gre->SetPointError(3,0,0.07);
+ gre->SetPoint(4,1.125,3.33);
+ gre->SetPointError(4,0,0.07);
+ gre->SetPoint(5,1.375,3.53);
+ gre->SetPointError(5,0,0.07);
+ gre->SetPoint(6,1.625,3.46);
+ gre->SetPointError(6,0,0.07);
+ gre->SetPoint(7,1.875,3.41);
+ gre->SetPointError(7,0,0.07);
+ gre->SetPoint(8,-0.125,3.14);
+ gre->SetPointError(8,0,0.07);
+ gre->SetPoint(9,-0.375,3.04);
+ gre->SetPointError(9,0,0.07);
+ gre->SetPoint(10,-0.625,3.17);
+ gre->SetPointError(10,0,0.07);
+ gre->SetPoint(11,-0.875,3.33);
+ gre->SetPointError(11,0,0.07);
+ gre->SetPoint(12,-1.125,3.33);
+ gre->SetPointError(12,0,0.07);
+ gre->SetPoint(13,-1.375,3.53);
+ gre->SetPointError(13,0,0.07);
+ gre->SetPoint(14,-1.625,3.46);
+ gre->SetPointError(14,0,0.07);
+ gre->SetPoint(15,-1.875,3.41);
+ gre->SetPointError(15,0,0.07);
+ gre->Draw("p");
+
+ l->AddEntry(gre, "UA5 INEL", "P");
+
+ gre = new TGraphErrors(16);
+ gre->SetMarkerColor(4);
+ gre->SetFillColor(4);
+ gre->SetMarkerStyle(22);
+ gre->SetPoint(0,0.125,3.48);
+ gre->SetPointError(0,0,0.07);
+ gre->SetPoint(1,0.375,3.38);
+ gre->SetPointError(1,0,0.07);
+ gre->SetPoint(2,0.625,3.52);
+ gre->SetPointError(2,0,0.07);
+ gre->SetPoint(3,0.875,3.68);
+ gre->SetPointError(3,0,0.07);
+ gre->SetPoint(4,1.125,3.71);
+ gre->SetPointError(4,0,0.07);
+ gre->SetPoint(5,1.375,3.86);
+ gre->SetPointError(5,0,0.07);
+ gre->SetPoint(6,1.625,3.76);
+ gre->SetPointError(6,0,0.07);
+ gre->SetPoint(7,1.875,3.66);
+ gre->SetPointError(7,0,0.07);
+ gre->SetPoint(8,-0.125,3.48);
+ gre->SetPointError(8,0,0.07);
+ gre->SetPoint(9,-0.375,3.38);
+ gre->SetPointError(9,0,0.07);
+ gre->SetPoint(10,-0.625,3.52);
+ gre->SetPointError(10,0,0.07);
+ gre->SetPoint(11,-0.875,3.68);
+ gre->SetPointError(11,0,0.07);
+ gre->SetPoint(12,-1.125,3.71);
+ gre->SetPointError(12,0,0.07);
+ gre->SetPoint(13,-1.375,3.86);
+ gre->SetPointError(13,0,0.07);
+ gre->SetPoint(14,-1.625,3.76);
+ gre->SetPointError(14,0,0.07);
+ gre->SetPoint(15,-1.875,3.66);
+ gre->SetPointError(15,0,0.07);
+ gre->Draw("p");
+
+ l->AddEntry(gre, "UA5 NSD", "P");
+ }
+
+ l->Draw();
+
+ if (fileName2 == 0)
+ return;
+
+ new TCanvas;
+ gPad->SetGridx();
+ gPad->SetGridy();
+
+ if (errorsCorrelated)
+ {
+ for (Int_t i=1; i<=histESD[1]->GetNbinsX(); i++)
+ {
+ histESD[1]->SetBinError(i, 0);
+ histESDnsd[1]->SetBinError(i, 0);
+ }
+ }
+
+ histESD[0]->Divide(histESD[0], histESD[1]);
+ histESDnsd[0]->Divide(histESDnsd[0], histESDnsd[1]);
+
+ for (Int_t i=1; i<=histESD[1]->GetNbinsX(); i++)
+ histESDnsd[0]->SetBinContent(i, histESDnsd[0]->GetBinContent(i) + 0.2);
+
+ hist->DrawCopy()->GetYaxis()->SetRangeUser(0.8, 1.4);
+ histESD[0]->Draw("SAME");
+ histESDnsd[0]->Draw("SAME");
+}
+
TH1* DrawdNdEtaRatio(TH1* corr, TH1* mc, const char* name, Float_t etaPlotLimit)
{
TCanvas* canvas3 = new TCanvas(name, name, 600, 600);
AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
dNdEtaCorrection->LoadHistograms();
- TH1* hist = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("x");
- TH1* hist2 = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
+ TH1* hist = dNdEtaCorrection->GetTriggerBiasCorrectionOnePart()->GetEventCorrection()->Get1DCorrection("x");
+ TH1* hist2 = dNdEtaCorrection->GetTriggerBiasCorrectionOnePart()->GetEventCorrection()->Get1DCorrection("y", -5, 5);
TCanvas* canvas = new TCanvas("TriggerBias1D", "TriggerBias1D", 1000, 500);
canvas->Divide(2, 1);
TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
pave->SetFillColor(0);
- pave->AddText("|z| < 10 cm");
+ pave->AddText("|z| < 5 cm");
pave->Draw();
+
+ Float_t triggerEff = 100.0 / hist2->GetBinContent(1);
+ Printf("trigger eff in 0 bin is: %.2f %%", triggerEff);
+
+ return;
- canvas->SaveAs("TriggerBias1D.eps");
+ TH1* hist2 = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
+ //new TCanvas;
+ //hist2->Draw();
+
+ Printf("vertex reco eff in 0 bin is: %.2f %%", 100.0 / hist2->GetBinContent(1));
+
+ Printf("combined efficiency is %.2f %%", triggerEff / hist2->GetBinContent(1));
}
void VtxRecon()
c->SaveAs(Form("%s.eps", canvas->GetName()));
}
-void Track2Particle2DCreatePlots(const char* fileName = "correction_map.root")
+void Track2Particle2DCreatePlots(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
{
TFile::Open(fileName);
- AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
+ AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folder, folder);
dNdEtaCorrection->LoadHistograms();
TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
{
gSystem->Load("libPWG0base");
- Track2Particle2DCreatePlots(fileName);
+ Track2Particle2DCreatePlots(fileName, folder);
TH2* corrYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx"));
TH2* corrZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx"));
TH3* hist1 = (TH3*) dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("mc");
hist1->SetTitle("mc");
Printf("mc contains %f entries", hist1->Integral());
- Printf("mc contains %f entries in |vtx-z| < 10, pt > 0.3", hist1->Integral(hist1->GetXaxis()->FindBin(-9.9), hist1->GetXaxis()->FindBin(9.9), hist1->GetYaxis()->FindBin(-0.99), hist1->GetYaxis()->FindBin(0.99), hist1->GetZaxis()->FindBin(ptMin), hist1->GetNbinsZ()));
+ Printf("mc contains %f entries in |vtx-z| < 10, |eta| < 1, pt > 0.3", hist1->Integral(hist1->GetXaxis()->FindBin(-9.9), hist1->GetXaxis()->FindBin(9.9), hist1->GetYaxis()->FindBin(-0.99), hist1->GetYaxis()->FindBin(0.99), hist1->GetZaxis()->FindBin(ptMin), hist1->GetNbinsZ()));
TH3* hist2 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd");
hist2->SetTitle("esd");
Printf("esd contains %f entries", hist2->Integral());
- Printf("esd contains %f entries in |vtx-z| < 10, pt > 0.3", hist2->Integral(hist2->GetXaxis()->FindBin(-9.9), hist2->GetXaxis()->FindBin(9.9), hist2->GetYaxis()->FindBin(-0.99), hist2->GetYaxis()->FindBin(0.99), hist2->GetZaxis()->FindBin(ptMin), hist2->GetNbinsZ()));
+ Printf("esd contains %f entries in |vtx-z| < 10, |eta| < 1, pt > 0.3", hist2->Integral(hist2->GetXaxis()->FindBin(-9.9), hist2->GetXaxis()->FindBin(9.9), hist2->GetYaxis()->FindBin(-0.99), hist2->GetYaxis()->FindBin(0.99), hist2->GetZaxis()->FindBin(ptMin), hist2->GetNbinsZ()));
AliPWG0Helper::CreateDividedProjections(hist1, hist2);
AliPWG0Helper::CreateDividedProjections(hist1, hist2, "x");
PrintHist2D(stats);
PrintHist(proj);
+ Float_t ua5_SD = 0.153;
+ Float_t ua5_DD = 0.080;
+ Float_t ua5_ND = 0.767;
+
+ Printf("+++ FRACTIONS +++");
+
+ Printf("ND: %f", proj->GetBinContent(3) / proj->GetBinContent(1));
+ Printf("SD: %f", proj->GetBinContent(4) / proj->GetBinContent(1));
+ Printf("DD: %f", proj->GetBinContent(5) / proj->GetBinContent(1));
+
Printf("+++ TRIGGER EFFICIENCIES +++");
Printf("INEL = %.1f", 100. * (proj->GetBinContent(1) - stats->GetBinContent(1, 1) - stats->GetBinContent(1, 3)) / proj->GetBinContent(1));
Printf("NSD = %.1f", 100. * (proj->GetBinContent(2) - stats->GetBinContent(2, 1) - stats->GetBinContent(2, 3)) / proj->GetBinContent(2));
- Printf("ND = %.1f", 100. * (proj->GetBinContent(3) - stats->GetBinContent(3, 1) - stats->GetBinContent(3, 3)) / proj->GetBinContent(3));
- Printf("SD = %.1f", 100. * (proj->GetBinContent(4) - stats->GetBinContent(4, 1) - stats->GetBinContent(4, 3)) / proj->GetBinContent(4));
- Printf("DD = %.1f", 100. * (proj->GetBinContent(5) - stats->GetBinContent(5, 1) - stats->GetBinContent(5, 3)) / proj->GetBinContent(5));
+
+
+ Float_t trigND = 100. * (proj->GetBinContent(3) - stats->GetBinContent(3, 1) - stats->GetBinContent(3, 3)) / proj->GetBinContent(3);
+ Float_t trigSD = 100. * (proj->GetBinContent(4) - stats->GetBinContent(4, 1) - stats->GetBinContent(4, 3)) / proj->GetBinContent(4);
+ Float_t trigDD = 100. * (proj->GetBinContent(5) - stats->GetBinContent(5, 1) - stats->GetBinContent(5, 3)) / proj->GetBinContent(5);
+
+ Printf("ND = %.1f", trigND);
+ Printf("SD = %.1f", trigSD);
+ Printf("DD = %.1f", trigDD);
+
+ Float_t trigINELUA5 = ua5_SD * trigSD + ua5_DD * trigDD + ua5_ND * trigND;
+ Float_t trigNSDUA5 = (ua5_DD * trigDD + ua5_ND * trigND) / (ua5_DD + ua5_ND);
+ Printf("INEL (UA5) = %.1f", trigINELUA5);
+ Printf("NSD (UA5) = %.1f", trigNSDUA5);
Printf("+++ VERTEX EFFICIENCIES +++");
Printf("SD = %.1f", vtxSD);
Printf("DD = %.1f", vtxDD);
- Float_t ua5_SD = 0.153;
- Float_t ua5_DD = 0.080;
- Float_t ua5_ND = 0.767;
-
Float_t vtxINELUA5 = ua5_SD * vtxSD + ua5_DD * vtxDD + ua5_ND * vtxND;
Float_t vtxNSDUA5 = (ua5_DD * vtxDD + ua5_ND * vtxND) / (ua5_DD + ua5_ND);
Printf("INEL (UA5) = %.1f", vtxINELUA5);
}
}
+void CompareDiamond(const char* mc, const char* data)
+{
+ TFile::Open(mc);
+
+ hist = (TH3*) gFile->Get("vertex_check");
+
+ gStyle->SetOptFit(1);
+
+ TH1* proj[3];
+ proj[0] = hist->ProjectionX("vertex_check_px");
+ proj[1] = hist->ProjectionY("vertex_check_py");
+ proj[2] = hist->ProjectionZ("vertex_check_pz");
+
+ TFile::Open(data);
+
+ hist = (TH3*) gFile->Get("vertex_check");
+
+ TH1* proj2[3];
+ proj2[0] = hist->ProjectionX("vertex_check_px2");
+ proj2[1] = hist->ProjectionY("vertex_check_py2");
+ proj2[2] = hist->ProjectionZ("vertex_check_pz2");
+
+ for (Int_t i=0; i<3; i++)
+ {
+ CompareQualityHists(proj[i], proj2[i], 1, 1);
+ }
+}
+
void FitDiamondVsMult()
{
TFile::Open("analysis_esd_raw.root");
file2 = TFile::Open(fileName2);
hist2 = (TH1*) file2->Get(plotName);
+ hist1->SetStats(0);
+
+ Printf("Entries in histograms: %f %f", hist1->Integral(), hist2->Integral());
+
if (exec)
{
hist1 = (TH1*) gROOT->ProcessLine(Form(exec, hist1, "hist1a"));
hist2 = (TH1*) gROOT->ProcessLine(Form(exec, hist2, "hist2a"));
+ hist1->Sumw2();
+ hist2->Sumw2();
+ Printf("Entries in histograms: %f %f", hist1->Integral(), hist2->Integral());
}
CompareQualityHists(hist1, hist2, rebin1, rebin2);
hist2->Rebin(TMath::Abs(rebin2));
}
+ //hist2 = hist2->Rebin(hist1->GetNbinsX(), Form("%s_rebinned", hist2->GetName()), hist1->GetXaxis()->GetXbins()->GetArray());
+
+ //hist1->Scale(1.0 / 0.83);
+
+//hist1->GetXaxis()->SetRangeUser(0, 50);
+/* hist1->GetYaxis()->SetRangeUser(0.9, 1.2);
+ hist1->Scale(1.0 / 0.808751);*/
+
+ //hist1->Scale(1.0 / 1.24632);
+ //hist1->Scale(1.0 / 1.23821);
+ //hist1->Scale(1.0 / 1.26213);
+
if (rebin1 > 0 && rebin2 > 0)
{
- hist1->Scale(hist2->Integral() / hist1->Integral() / rebin1 * rebin2);
-
+ hist1->Scale(hist2->Integral() / hist1->Integral() * hist2->GetXaxis()->GetBinWidth(1) / hist1->GetXaxis()->GetBinWidth(1) / rebin1 * rebin2);
+
//hist1->Scale(0.5);
//hist2->Scale(0.5);
}
c = new TCanvas;
- hist1->GetYaxis()->SetRangeUser(0, hist1->GetMaximum() * 1.3);
- hist1->DrawCopy();
- hist2->DrawCopy("SAME");
- c->SaveAs(Form("%s_1.png", hist1->GetName()));
+ if (strcmp(hist1->GetName(), "fDeltaTheta") == 0 || strcmp(hist1->GetName(), "fDeltaPhi") == 0 || strcmp(hist1->GetName(), "fMultVtx") == 0 || TString(hist1->GetName()).BeginsWith("vertex_check"))
+ c->SetLogy();
+
+ if (TString(hist1->GetName()).BeginsWith("fMultiplicityESD"))
+ {
+ c->SetLogy();
+ loadlibs();
+ AliPWG0Helper::NormalizeToBinWidth(hist1);
+ AliPWG0Helper::NormalizeToBinWidth(hist2);
+ }
+
+ Printf("Means: %f %f %e", hist1->GetMean(), hist2->GetMean(), 1.0 - hist2->GetMean() / hist1->GetMean());
- for (Int_t i=1; i<=hist1->GetNbinsX(); i++)
- if (hist1->GetBinContent(i) == 0 && hist2->GetBinContent(i) > 0 || hist1->GetBinContent(i) > 0 && hist2->GetBinContent(i) == 0)
- Printf("Inconsistent bin %d: %f %f", i, hist1->GetBinContent(i), hist2->GetBinContent(i));
+ //hist1->GetYaxis()->SetRangeUser(0.01, hist1->GetMaximum() * 1.3);
+ hist1->DrawCopy("HISTE");
+ hist2->DrawCopy("HISTE SAME");
+ gPad->SetGridx();
+ gPad->SetGridy();
+ //gPad->SetLogy();
+ c->SaveAs(Form("%s_1.png", hist1->GetName()));
if (rebin1 == rebin2)
{
+ for (Int_t i=1; i<=hist1->GetNbinsX(); i++)
+ if (hist1->GetBinContent(i) == 0 && hist2->GetBinContent(i) > 0 || hist1->GetBinContent(i) > 0 && hist2->GetBinContent(i) == 0)
+ Printf("Inconsistent bin %d: %f %f", i, hist1->GetBinContent(i), hist2->GetBinContent(i));
+
c2 = new TCanvas;
+ hist1->GetYaxis()->SetRangeUser(0.5, 1.5);
hist1->Divide(hist2);
- hist1->DrawCopy();
+ hist1->DrawCopy("HIST");
+ gPad->SetGridx();
+ gPad->SetGridy();
c2->SaveAs(Form("%s_2.png", hist1->GetName()));
+ /*
for (Int_t i=1; i<=hist1->GetNbinsX(); i++)
if (hist1->GetBinContent(i) > 0.9 && hist1->GetBinContent(i) < 1.1)
hist1->SetBinContent(i, 0);
new TCanvas;
hist1->SetMarkerStyle(20);
hist1->DrawCopy("P");
+ */
}
}
c->SaveAs("clusters_vs_tracklets.eps");
}
-
+
void VertexPlotBackgroundNote()
{
TFile::Open("all.root");
proj->SetLineColor(2);
proj->Draw("SAME");
-
-
}
void BackgroundAnalysis(const char* signal, const char* background)
c->SaveAs(Form("%s/stats.png", list[i]));
}
}
+
+void CompareMCDataTrigger(const char* mcFile, const char* dataFile)
+{
+ TH1* stat[2];
+
+ TFile::Open(mcFile);
+ mc = (TH1*) gFile->Get("trigger_histograms_/fHistFiredBitsSPD");
+ stat[0] = (TH1*) gFile->Get("fHistStatistics");
+
+ TFile::Open(dataFile);
+ data = (TH1*) gFile->Get("trigger_histograms_+CINT1B-ABCE-NOPF-ALL/fHistFiredBitsSPD");
+ if (!data)
+ data = (TH1*) gFile->Get("trigger_histograms_+CSMBB-ABCE-NOPF-ALL/fHistFiredBitsSPD");
+
+ stat[1] = (TH1*) gFile->Get("fHistStatistics");
+
+ CompareQualityHists(mc, data);
+
+ for (Int_t i=0; i<2; i++)
+ {
+ Float_t total = stat[i]->GetBinContent(stat[i]->GetXaxis()->FindBin("Trigger class"), 1);
+ Float_t spd = stat[i]->GetBinContent(stat[i]->GetXaxis()->FindBin("FO >= 2"), 1);
+ Float_t v0A = stat[i]->GetBinContent(stat[i]->GetXaxis()->FindBin("V0A"), 1);
+ Float_t v0C = stat[i]->GetBinContent(stat[i]->GetXaxis()->FindBin("V0C"), 1);
+
+ Printf("%s:\nSPD / V0A: %.3f\nSPD / V0C: %.3f\nV0A / V0C: %.3f", (i == 0) ? "MC " : "Data", spd / v0A, spd / v0C, v0A / v0C);
+ Printf("SPD / Total: %.3f\nV0A / Total: %.3f\nV0C / Total: %.3f\n", spd / total, v0A / total, v0C / total);
+ }
+}
+
+void CompareMCDatadNdEta(const char* mcFile, const char* dataFile)
+{
+ //CompareQualityHists(mcFile, dataFile, "fEtaPhi", 4, 4, "((TH2*)%p)->ProjectionY(\"%s\", 1, 40)");
+ //CompareQualityHists(mcFile, dataFile, "fEtaPhi", 4, 4, "((TH2*)%p)->ProjectionY(\"%s\", 41, 80)");
+
+ CompareQualityHists(mcFile, dataFile, "fEtaPhi", 1, 1, "((TH2*)%p)->ProjectionX(\"%s\", 271, 360)");
+}
+
+void TrigVsTrigVtx(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
+{
+ loadlibs();
+ if (!TFile::Open(fileName))
+ return;
+
+ AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
+ if (!dNdEtaCorrection->LoadHistograms())
+ return;
+
+ TH2* eTrig = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetGeneratedHistogram();
+ TH2* eTrigVtx = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetMeasuredHistogram();
+
+ eTrig_proj = eTrig->ProjectionY("y1", eTrig->GetYaxis()->FindBin(-9.9), eTrig->GetYaxis()->FindBin(9.9));
+ eTrigVtx_proj = eTrigVtx->ProjectionY("y2", eTrig->GetYaxis()->FindBin(-9.9), eTrig->GetYaxis()->FindBin(9.9));
+
+ new TCanvas;
+ eTrig_proj->Draw();
+ eTrig_proj->GetXaxis()->SetRangeUser(0, 20);
+ eTrigVtx_proj->SetLineColor(2);
+ eTrigVtx_proj->Draw("SAME");
+
+ gPad->SetLogy();
+}
+
+void PrintAverageNSDCorrectionFactors()
+{
+ // factors estimated from MC, can be slighly different with data b/c correction is applies as function of measured multiplicity
+
+ loadlibs();
+
+ if (!TFile::Open("correction_mapprocess-types.root"))
+ return;
+
+ AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction_ND", "dndeta_correction_ND");
+ if (!dNdEtaCorrection->LoadHistograms())
+ return;
+
+ AlidNdEtaCorrection* dNdEtaCorrection2 = new AlidNdEtaCorrection("dndeta_correction_DD", "dndeta_correction_DD");
+ if (!dNdEtaCorrection2->LoadHistograms())
+ return;
+
+ // for scaling factors see drawSystematics.C; GetRelativeFractions()
+ // 900 GeV
+ //dNdEtaCorrection->Scale(1.06);
+ //dNdEtaCorrection->Add(dNdEtaCorrection2, 9.5 / 12.3);
+ // 2.36 TeV
+ dNdEtaCorrection->Scale(1.036);
+ dNdEtaCorrection->Add(dNdEtaCorrection2, 0.075 * 1.43 / 0.127);
+
+ Printf("event adding: %f", dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral() / dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral());
+
+ Printf("track adding: %f", dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetGeneratedHistogram()->Integral() / dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetMeasuredHistogram()->Integral());
+
+ AlidNdEtaCorrection* dNdEtaCorrection3 = new AlidNdEtaCorrection("dndeta_correction_SD", "dndeta_correction_SD");
+ if (!dNdEtaCorrection3->LoadHistograms())
+ return;
+
+ // 900 GeV
+ //dNdEtaCorrection3->Scale(0.153 / 0.189);
+ // 2.36 TeV
+ dNdEtaCorrection3->Scale(0.159 / 0.166);
+ dNdEtaCorrection->Add(dNdEtaCorrection3);
+
+ Printf("event subtraction: %f", dNdEtaCorrection3->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() / dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral());
+
+ Printf("track subtraction: %f", dNdEtaCorrection3->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetMeasuredHistogram()->Integral() / dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetMeasuredHistogram()->Integral());
+
+ dNdEtaCorrection->GetTriggerBiasCorrectionNSD()->PrintInfo(0.0);
+}
+
\ No newline at end of file
break;
// one species enhanced / reduced
- case 2: // + 50% kaons
- case 3: // - 50% kaons
- case 4: // + 50% protons
- case 5: // - 50% protons
- case 6: // + 50% kaons + 50% protons
- case 7: // - 50% kaons - 50% protons
- case 8: // + 50% kaons - 50% protons
- case 9: // - 50% kaons + 50% protons
+ case 2: // + 30% kaons
+ case 3: // - 30% kaons
+ case 4: // + 30% protons
+ case 5: // - 30% protons
+ case 6: // + 30% kaons + 30% protons
+ case 7: // - 30% kaons - 30% protons
+ case 8: // + 30% kaons - 30% protons
+ case 9: // - 30% kaons + 30% protons
+ case 10: // + 30% others
+ case 11: // - 30% others
TString* str = new TString;
if (index < 6)
{
Int_t correctionIndex = index / 2;
- Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5;
+ Double_t scaleFactor = (index % 2 == 0) ? 1.3 : 0.7;
fdNdEtaCorrection[correctionIndex]->GetTrack2ParticleCorrection()->GetTrackCorrection()->Scale(scaleFactor);
str->Form("%s%s", (correctionIndex == 0) ? "Pi" : ((correctionIndex == 1) ? "K" : (correctionIndex == 2) ? "p" : "others"), (index % 2 == 0) ? "Boosted" : "Reduced");
}
- else
+ else if (index < 10)
{
- Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5;
+ Double_t scaleFactor = (index % 2 == 0) ? 1.3 : 0.7;
fdNdEtaCorrection[1]->GetTrack2ParticleCorrection()->GetTrackCorrection()->Scale(scaleFactor);
str->Form("%s%s", "K", (scaleFactor > 1) ? "Boosted" : "Reduced");
if (index >= 8)
- scaleFactor = (index % 2 == 0) ? 0.5 : 1.5;
+ scaleFactor = (index % 2 == 0) ? 0.3 : 1.7;
fdNdEtaCorrection[2]->GetTrack2ParticleCorrection()->GetTrackCorrection()->Scale(scaleFactor);
*str += Form("%s%s", "p", (scaleFactor > 1) ? "Boosted" : "Reduced");
}
+ else
+ {
+ Double_t scaleFactor = (index % 2 == 0) ? 1.3 : 0.7;
+ fdNdEtaCorrection[3]->GetTrack2ParticleCorrection()->GetTrackCorrection()->Scale(scaleFactor);
+ str->Form("%s%s", "others", (scaleFactor > 1) ? "Boosted" : "Reduced");
+ }
return str->Data();
break;
const char* names[] = { "pi", "K", "p", "other" };
TH1* hRatios[20];
- Int_t nCompositions = 10;
+ //backgroundEvents = 1162+434; // Michele for MB1, run 104892, 15.02.10
+ backgroundEvents = -1; // use 0 bin from MC! for 2.36 TeV
+
+ Printf("Subtracting %d background events!!!", backgroundEvents);
+ gSystem->Sleep(1000);
+
+ Int_t nCompositions = 12;
Int_t counter = 0;
for (Int_t comp = 1; comp < nCompositions; ++comp)
{
// skip "other" particle correction here
// with them has to be dealt differently, maybe just increasing the neutral particles...
- for (Int_t i=1; i<3; ++i)
+ for (Int_t i=1; i<4; ++i)
collection->Add(fdNdEtaCorrection[i]);
fdNdEtaCorrection[0]->Merge(collection);
canvas2->SaveAs("particlecomposition_result.eps");
}
-void mergeCorrectionsWithDifferentCrosssections(Int_t correctionTarget = 3, Char_t* correctionFileName="correction_mapprocess-types.root", const char* analysisFileName = "analysis_esd_raw.root", const Char_t* outputFileName="systematics_vtxtrigger_compositions.root") {
+void mergeCorrectionsWithDifferentCrosssections(Int_t origin, Int_t correctionTarget = 3, Char_t* correctionFileName="correction_mapprocess-types.root", const char* analysisFileName = "analysis_esd_raw.root", const Char_t* outputFileName=0) {
//
// Function used to merge standard corrections with vertex
// reconstruction corrections obtained by a certain mix of ND, DD
// correctionTarget is of type AlidNdEtaCorrection::CorrectionType
// kINEL = 3
// kNSD = 4
+ // kOnePart = 6
+
+ if (outputFileName == 0)
+ {
+ if (correctionTarget == 3)
+ outputFileName = "systematics_vtxtrigger_compositions_inel.root";
+ if (correctionTarget == 4)
+ outputFileName = "systematics_vtxtrigger_compositions_nsd.root";
+ if (correctionTarget == 6)
+ outputFileName = "systematics_vtxtrigger_compositions_onepart.root";
+ }
loadlibs();
//Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.5, 0.5, 1.5, 0.5, 0.5, 1.5, 1.0, 1.0, 1.25, 0.75, 1.25, 0.75, 0.75, 1.25};
//Float_t scalesDD[] = {1.0, 1.4, 0.6, 1.0, 1.0, 1.4, 0.6, 1.4, 0.6, 1.25, 0.75, 1.0, 1.0, 1.25, 0.75, 1.25, 0.75};
//Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.4, 0.6, 1.4, 0.6, 0.4, 1.6, 1.0, 1.0, 1.25, 0.75, 1.25, 0.75, 0.75, 1.25};
- Int_t nChanges = 9;
+/* Int_t nChanges = 9;
const Char_t* changes[] = { "ua5","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdlessddmore", "sdmoreddless" };
Float_t scalesND[] = {0.767, 0.767, 0.767, 0.767, 0.767, 0.767, 0.767, 0.767, 0.767};
Float_t scalesDD[] = {0.080, 0.130, 0.030, 0.080, 0.080, 0.130, 0.030, 0.130, 0.030};
- Float_t scalesSD[] = {0.153, 0.153, 0.153, 0.203, 0.103, 0.203, 0.103, 0.103, 0.203};
+ Float_t scalesSD[] = {0.153, 0.153, 0.153, 0.203, 0.103, 0.203, 0.103, 0.103, 0.203};*/
+
+ Float_t ref_SD = -1;
+ Float_t ref_DD = -1;
+ Float_t ref_ND = -1;
+
+ Float_t error_SD = -1;
+ Float_t error_DD = -1;
+ Float_t error_ND = -1;
+
+ GetRelativeFractions(origin, ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND);
+
+ Printf("Using x-sections:\n SD: %f +- %f\n DD: %f +- %f\n ND: %f +- %f", ref_SD, error_SD, ref_DD, error_DD, ref_ND, error_ND);
+
+ const Char_t* changes[] = { "default","sdless","sdmore","ddless","ddmore", "dless", "dmore", "sdlessddmore", "sdmoreddless" };
+ Int_t nChanges = 9;
+ Float_t scalesSD[9];
+ Float_t scalesDD[9];
+ Float_t scalesND[9];
- for (Int_t i=0; i<9; i++)
- scalesND[i] = 1.0 - scalesDD[i] - scalesSD[i];
+ if (1)
+ {
+ // sample 8 points on the error ellipse
+ for (Int_t i=0; i<9; i++)
+ {
+ Float_t factorSD = 0;
+ Float_t factorDD = 0;
+
+ if (i > 0 && i < 3)
+ factorSD = (i % 2 == 0) ? 1 : -1;
+ else if (i >= 3 && i < 5)
+ factorDD = (i % 2 == 0) ? 1 : -1;
+ else if (i >= 5 && i < 9)
+ {
+ factorSD = ((i % 2 == 0) ? 1.0 : -1.0) / TMath::Sqrt(2);
+ if (i == 5 || i == 6)
+ factorDD = factorSD;
+ else
+ factorDD = -factorSD;
+ }
+
+ scalesSD[i] = ref_SD + factorSD * error_SD;
+ scalesDD[i] = ref_DD + factorDD * error_DD;
+ scalesND[i] = 1.0 - scalesDD[i] - scalesSD[i];
+
+ Printf("Case %d: SD: %f DD: %f ND: %f", i, scalesSD[i], scalesDD[i], scalesND[i]);
+ }
+ }
+ else
+ {
+ Printf("WARNING: Special treatment for ratios active");
+ gSystem->Sleep(1000);
+
+ // constrained values by allowed changing of cross-sections
+ Float_t pythiaScaling = 0.224 / 0.189;
+
+ if (origin == 10)
+ {
+ // 900 GeV
+ for (Int_t i=0; i<9; i++)
+ {
+ scalesSD[i] = 15.3;
+ scalesDD[i] = 9.5;
+ }
+
+ scalesSD[1] = 15.7;
+ scalesSD[2] = 17.6;
+ scalesSD[3] = 13.5;
+ scalesSD[4] = 17.6;
+
+ scalesDD[5] = 15.5;
+ scalesDD[6] = 8.8;
+ scalesDD[7] = 13.8;
+ scalesDD[8] = 7.6;
+ }
+ else if (origin == 20)
+ {
+ // 2.36 TeV
+ pythiaScaling = 0.217 / 0.167;
+
+ for (Int_t i=0; i<9; i++)
+ {
+ scalesSD[i] = 15.9;
+ scalesDD[i] = 10.7;
+ }
+
+ scalesSD[1] = 13.5;
+ scalesSD[2] = 15.2;
+ scalesSD[3] = 13.5;
+ scalesSD[4] = 17.6;
+
+ scalesDD[5] = 13.8;
+ scalesDD[6] = 7.6;
+ scalesDD[7] = 13.8;
+ scalesDD[8] = 7.6;
+ }
+ else
+ AliFatal("Not supported");
+ for (Int_t i=0; i<9; i++)
+ {
+ scalesSD[i] /= 100;
+ scalesSD[i] *= pythiaScaling;
+ scalesDD[i] /= 100;
+ scalesND[i] = 1.0 - scalesDD[i] - scalesSD[i];
+ Printf("Case %d: SD: %f DD: %f ND: %f", i, scalesSD[i], scalesDD[i], scalesND[i]);
+ }
+ }
+
+ Int_t backgroundEvents = 0;
+
+ //backgroundEvents = 1162+434; // Michele for MB1, run 104892, 15.02.10
+ //backgroundEvents = 6; // Michele for V0AND, run 104892, 15.02.10
+
+ //backgroundEvents = 4398+961; // Michele for MB1, run 104824-52, 16.02.10
+ //backgroundEvents = 19; // Michele for V0AND, run 104824-52, 16.02.10
+
+ backgroundEvents = -1; // use 0 bin from MC! for 2.36 TeV
+
+ Printf("Subtracting %d background events!!!", backgroundEvents);
+ gSystem->Sleep(1000);
+
/*
const Char_t* changes[] = { "pythia", "qgsm", "phojet"};
Float_t scalesND[] = {1.0, 1.10, 1.11};
// cross section from Pythia
// 14 TeV!
- Float_t sigmaND = 55.2;
- Float_t sigmaDD = 9.78;
- Float_t sigmaSD = 14.30;
+// Float_t sigmaND = 55.2;
+// Float_t sigmaDD = 9.78;
+// Float_t sigmaSD = 14.30;
// standard correction
TFile::Open(correctionFileName);
correctionStandard->GetTriggerBiasCorrectionINEL()->Reset();
correctionStandard->GetTriggerBiasCorrectionNSD()->Reset();
correctionStandard->GetTriggerBiasCorrectionND()->Reset();
+ correctionStandard->GetTriggerBiasCorrectionOnePart()->Reset();
AlidNdEtaCorrection* corrections[100];
TH1F* hRatios[100];
if (j == 1 || j == 2)
{
dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->Scale(scaleND);
- dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->Scale(scalesDD[i]);
- dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->Scale(scalesSD[i]);
+ dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->Scale(scaleDD);
+ dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->Scale(scaleSD);
dNdEtaCorrectionND->GetTriggerBiasCorrectionNSD()->Scale(scaleND);
dNdEtaCorrectionDD->GetTriggerBiasCorrectionNSD()->Scale(scaleDD);
dNdEtaCorrectionND->GetTriggerBiasCorrectionND()->Scale(scaleND);
dNdEtaCorrectionDD->GetTriggerBiasCorrectionND()->Scale(scaleDD);
dNdEtaCorrectionSD->GetTriggerBiasCorrectionND()->Scale(scaleSD);
+
+ dNdEtaCorrectionND->GetTriggerBiasCorrectionOnePart()->Scale(scaleND);
+ dNdEtaCorrectionDD->GetTriggerBiasCorrectionOnePart()->Scale(scaleDD);
+ dNdEtaCorrectionSD->GetTriggerBiasCorrectionOnePart()->Scale(scaleSD);
}
//clear track in correction
current->Merge(&collection);
current->Finish();
+
+ // print 0 bin efficiency
+ TH1* hist2 = current->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
+ if (hist2->GetBinContent(1))
+ {
+ Float_t triggerEff = 100.0 / hist2->GetBinContent(1);
+ Printf("trigger eff in 0 bin is: %.2f %%", triggerEff);
+ }
corrections[counter] = current;
dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD");
fdNdEtaAnalysis->LoadHistograms();
- fdNdEtaAnalysis->Finish(current, 0.3, correctionTarget, Form("%d %d", j, i));
+ fdNdEtaAnalysis->Finish(current, 0.151, correctionTarget, Form("%d %d", j, i), backgroundEvents);
name = "ratio";
if (j==0) name.Append("_vetexReco_");
name.Form("ND #times %0.2f DD #times %0.2f, SD #times %0.2f",scaleND,scaleDD,scaleSD);
hRatios[counter]->SetTitle(name.Data());
hRatios[counter]->SetYTitle("dN_{ch}/d#eta ratio #frac{default cross-section}{modified cross-sections}");
-
+
+ TF1* pol0 = new TF1("pol0", "[0]", -0.49, 0.49);
+ pol0->SetParameter(0, 0);
+ hRatios[counter]->Fit(pol0, "RN");
+ Printf("Case %d: %f", i, pol0->GetParameter(0));
+
if (counter > 0)
hRatios[counter]->Divide(hRatios[0],hRatios[counter],1,1);
fout->Close();
}
-void CreateCorrectionsWithUA5CrossSections(Int_t origin, const Char_t* correctionFileName="correction_mapprocess-types.root", const Char_t* outputFileName="correction_map2.root") {
- //
- // Function used to merge standard corrections with vertex
- // reconstruction corrections obtained by a certain mix of ND, DD
- // and SD events.
- //
+void GetRelativeFractions(Int_t origin, Float_t& ref_SD, Float_t& ref_DD, Float_t& ref_ND, Float_t& error_SD, Float_t& error_DD, Float_t& error_ND)
+{
// origin:
// -1 = Pythia (test)
// 0 = UA5
// 3 = Durham
//
- loadlibs();
-
- const Char_t* typeName[] = { "vertexreco", "trigger", "vtxtrigger" };
-
- Float_t ref_SD = -1;
- Float_t ref_DD = -1;
- Float_t ref_ND = -1;
-
switch (origin)
{
- case -1: // Pythia, as test
+ case -10: // Pythia default at 900 GeV, 50% error
+ Printf("PYTHIA x-sections");
+ ref_SD = 0.223788; error_SD = ref_SD * 0.5;
+ ref_DD = 0.123315; error_DD = ref_DD * 0.5;
+ ref_ND = 0.652897; error_ND = 0;
+ break;
+
+ case -1: // Pythia default at 900 GeV, as test
+ Printf("PYTHIA x-sections");
ref_SD = 0.223788;
ref_DD = 0.123315;
ref_ND = 0.652897;
break;
case 0: // UA5
- ref_SD = 0.153;
- ref_DD = 0.080;
- ref_ND = 0.767;
+ Printf("UA5 x-sections a la first paper");
+ ref_SD = 0.153; error_SD = 0.05;
+ ref_DD = 0.080; error_DD = 0.05;
+ ref_ND = 0.767; error_ND = 0;
+ break;
+
+ case 10: // UA5
+ Printf("UA5 x-sections hadron level definition for Pythia");
+ // Fractions in Pythia with UA5 cuts selection for SD
+ // ND: 0.688662
+ // SD: 0.188588 --> this should be 15.3
+ // DD: 0.122750
+ ref_SD = 0.224 * 0.153 / 0.189; error_SD = 0.023 * 0.224 / 0.189;
+ ref_DD = 0.095; error_DD = 0.06;
+ ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
+ break;
+
+ case 11: // UA5
+ Printf("UA5 x-sections hadron level definition for Phojet");
+ // Fractions in Phojet with UA5 cuts selection for SD
+ // ND: 0.783573
+ // SD: 0.151601 --> this should be 15.3
+ // DD: 0.064827
+ ref_SD = 0.191 * 0.153 / 0.152; error_SD = 0.023 * 0.191 / 0.152;
+ ref_DD = 0.095; error_DD = 0.06;
+ ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
break;
+ case 20: // E710, 1.8 TeV
+ Printf("E710 x-sections hadron level definition for Pythia");
+ // ND: 0.705709
+ // SD: 0.166590 --> this should be 15.9
+ // DD: 0.127701
+ ref_SD = 0.217 * 0.159 / 0.167; error_SD = 0.024 * 0.217 / 0.167;
+ ref_DD = 0.075 * 1.43; error_DD = 0.02 * 1.43;
+ ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
+ break;
+
+ case 21: // E710, 1.8 TeV
+ Printf("E710 x-sections hadron level definition for Phojet");
+ // ND: 0.817462
+ // SD: 0.125506 --> this should be 15.9
+ // DD: 0.057032
+ ref_SD = 0.161 * 0.159 / 0.126; error_SD = 0.024 * 0.161 / 0.126;
+ ref_DD = 0.075 * 1.43; error_DD = 0.02 * 1.43;
+ ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
+ break;
+
case 1: // data 1.8 TeV
+ Printf("??? x-sections");
ref_SD = 0.152;
ref_DD = 0.092;
ref_ND = 1 - ref_SD - ref_DD;
break;
case 2: // tel-aviv model
+ Printf("Tel-aviv model x-sections");
ref_SD = 0.171;
ref_DD = 0.094;
ref_ND = 1 - ref_SD - ref_DD;
break;
case 3: // durham model
+ Printf("Durham model x-sections");
ref_SD = 0.190;
ref_DD = 0.125;
ref_ND = 1 - ref_SD - ref_DD;
break;
default:
- return;
+ AliFatal(Form("Unknown origin %d", origin));
}
-
- //Karel (UA5):
+}
+
+void CreateCorrectionsWithUA5CrossSections(Int_t origin, const Char_t* correctionFileName="correction_mapprocess-types.root", const Char_t* outputFileName="correction_map2.root") {
+ //
+ // Function used to merge standard corrections with vertex
+ // reconstruction corrections obtained by a certain mix of ND, DD
+ // and SD events.
+ //
+ loadlibs();
+
+ const Char_t* typeName[] = { "vertexreco", "trigger", "vtxtrigger" };
+
+ Float_t ref_SD = -1;
+ Float_t ref_DD = -1;
+ Float_t ref_ND = -1;
+
+ Float_t error_SD = -1;
+ Float_t error_DD = -1;
+ Float_t error_ND = -1;
+
+ GetRelativeFractions(origin, ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND);
+
+ Printf("Using x-sections:\n SD: %f +- %f\n DD: %f +- %f\n ND: %f +- %f", ref_SD, error_SD, ref_DD, error_DD, ref_ND, error_ND);
+
+//Karel (UA5):
// fsd = 0.153 +- 0.031
// fdd = 0.080 +- 0.050
// fnd = 0.767 +- 0.059
// Durham model Sd/Inel = 0.190 Dd/Inel = 0.125
// Data Sd/Inel = 0.152 +- 0.030 Dd/Inel = 0.092 +- 0.45
-
-
// standard correction
TFile::Open(correctionFileName);
AlidNdEtaCorrection* correctionStandard = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
correctionStandard->GetTriggerBiasCorrectionINEL()->Reset();
correctionStandard->GetTriggerBiasCorrectionNSD()->Reset();
correctionStandard->GetTriggerBiasCorrectionND()->Reset();
+ correctionStandard->GetTriggerBiasCorrectionOnePart()->Reset();
AlidNdEtaCorrection* corrections[100];
TH1F* hRatios[100];
dNdEtaCorrectionDD->GetTriggerBiasCorrectionND()->Scale(scaleDD);
dNdEtaCorrectionSD->GetTriggerBiasCorrectionND()->Scale(scaleSD);
+ dNdEtaCorrectionND->GetTriggerBiasCorrectionOnePart()->Scale(scaleND);
+ dNdEtaCorrectionDD->GetTriggerBiasCorrectionOnePart()->Scale(scaleDD);
+ dNdEtaCorrectionSD->GetTriggerBiasCorrectionOnePart()->Scale(scaleSD);
+
//clear track in correction
dNdEtaCorrectionND->GetTrack2ParticleCorrection()->Reset();
dNdEtaCorrectionDD->GetTrack2ParticleCorrection()->Reset();
current->Merge(&collection);
current->Finish();
+ // print 0 bin efficiency
+ TH1* hist2 = current->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
+ if (hist2->GetBinContent(1) > 0)
+ {
+ Float_t triggerEff = 100.0 / hist2->GetBinContent(1);
+ Printf("trigger eff in 0 bin is: %.2f %%", triggerEff);
+ }
+
TFile* fout = new TFile(outputFileName,"RECREATE");
current->SaveHistograms();
Printf(">>>> After");
corr->PrintInfo(0);
}
+
+Float_t FitAverage(TH1* hist, Int_t n, Float_t* begin, Float_t *end, Int_t color, Int_t& totalBins)
+{
+ Float_t average = 0;
+ totalBins = 0;
+
+ for (Int_t i=0; i<n; i++)
+ {
+ func = new TF1("func", "[0]", hist->GetXaxis()->GetBinLowEdge(hist->GetXaxis()->FindBin(begin[i])), hist->GetXaxis()->GetBinUpEdge(hist->GetXaxis()->FindBin(end[i])));
+ Int_t bins = hist->GetXaxis()->FindBin(end[i]) - hist->GetXaxis()->FindBin(begin[i]) + 1;
+ func->SetParameter(0, 1);
+ func->SetLineColor(color);
+
+ hist->Fit(func, "RNQ");
+ func->Draw("SAME");
+
+ average += func->GetParameter(0) * bins;
+ totalBins += bins;
+ }
+
+ return average / totalBins;
+}
+
+void SPDIntegrateGaps(Bool_t all, const char* mcFile = "../../../LHC10b2/v0or/spd/analysis_esd_raw.root")
+{
+ Float_t eta = 1.29;
+ Int_t binBegin = ((TH2*) gFile->Get("fEtaPhi"))->GetXaxis()->FindBin(-eta);
+ Int_t binEnd = ((TH2*) gFile->Get("fEtaPhi"))->GetXaxis()->FindBin(eta);
+
+ Printf("eta range: %f bins: %d %d", eta, binBegin, binEnd);
+
+ if (!all)
+ Printf("Eta smaller than 0 side");
+
+ c = new TCanvas;
+ TFile::Open("analysis_esd_raw.root");
+ hist = ((TH2*) gFile->Get("fEtaPhi"))->ProjectionY("hist", binBegin, (all) ? binEnd : 40);
+ hist->Rebin(2);
+ hist->SetStats(0);
+ hist->Sumw2();
+ hist->Draw("HIST");
+ gPad->SetGridx();
+ gPad->SetGridy();
+
+ TFile::Open(mcFile);
+ mcHist = ((TH2*) gFile->Get("fEtaPhi"))->ProjectionY("mcHist", binBegin, (all) ? binEnd : 40);
+ mcHist->Rebin(2);
+ mcHist->SetLineColor(2);
+ mcHist->Scale(hist->Integral() / mcHist->Integral());
+ mcHist->Draw("SAME");
+
+ Float_t add = 0;
+ Int_t bins;
+
+ Float_t okRangeBegin[] = { 0.04, 0.67, 1.34 };
+ Float_t okRangeEnd[] = { 0.55, 1.24, 1.63 };
+ Float_t gapRangeBegin[] = { 0.6, 1.27 };
+ Float_t gapRangeEnd[] = { 0.65, 1.32 };
+ Float_t averageOK = FitAverage(hist, 3, okRangeBegin, okRangeEnd, 1, bins);
+ Float_t averageGap = FitAverage(hist, 2, gapRangeBegin, gapRangeEnd, 2, bins);
+ add += bins * (averageOK - averageGap);
+ Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
+
+ Float_t okRangeBegin2[] = { 2.4, 2.65 };
+ Float_t okRangeEnd2[] = { 2.55, 3.2 };
+ Float_t gapRangeBegin2[] = { 2.59, 3.3 };
+ Float_t gapRangeEnd2[] = { 2.61, 3.3 };
+ averageOK = FitAverage(hist, 2, okRangeBegin2, okRangeEnd2, 1, bins);
+ averageGap = FitAverage(hist, 2, gapRangeBegin2, gapRangeEnd2, 2, bins);
+ add += bins * (averageOK - averageGap);
+ Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
+
+ Float_t okRangeBegin3[] = { 3.55, 3.9 };
+ Float_t okRangeEnd3[] = { 3.8, 4.15 };
+ Float_t gapRangeBegin3[] = { 3.83 };
+ Float_t gapRangeEnd3[] = { 3.86 };
+ averageOK = FitAverage(hist, 2, okRangeBegin3, okRangeEnd3, 1, bins);
+ averageGap = FitAverage(hist, 1, gapRangeBegin3, gapRangeEnd3, 2, bins);
+ add += bins * (averageOK - averageGap);
+ Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
+
+ Float_t okRangeBegin4[] = { 4.2, 4.5 };
+ Float_t okRangeEnd4[] = { 4.4, 4.7 };
+ Float_t gapRangeBegin4[] = { 4.45 };
+ Float_t gapRangeEnd4[] = { 4.45 };
+ averageOK = FitAverage(hist, 2, okRangeBegin4, okRangeEnd4, 1, bins);
+ averageGap = FitAverage(hist, 1, gapRangeBegin4, gapRangeEnd4, 2, bins);
+ add += bins * (averageOK - averageGap);
+ Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
+
+ Float_t okRangeBegin5[] = { 5.4, 5.7 };
+ Float_t okRangeEnd5[] = { 5.6, 5.8 };
+ Float_t gapRangeBegin5[] = { 5.63 };
+ Float_t gapRangeEnd5[] = { 5.67 };
+ averageOK = FitAverage(hist, 2, okRangeBegin5, okRangeEnd5, 1, bins);
+ averageGap = FitAverage(hist, 1, gapRangeBegin5, gapRangeEnd5, 2, bins);
+ add += bins * (averageOK - averageGap);
+ Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
+
+ Printf("This adds %.2f %% to the total number of tracklets (%f)", 100.0 * add / hist->Integral(), hist->Integral());
+ c->SaveAs("gap1.png");
+
+ add1 = add;
+ total1 = hist->Integral();
+
+ if (all)
+ return;
+
+ Printf("\nEta larger than 0 side");
+
+ c = new TCanvas;
+ TFile::Open("analysis_esd_raw.root");
+ hist = ((TH2*) gFile->Get("fEtaPhi"))->ProjectionY("hist2", 41, binEnd);
+ hist->Rebin(2);
+ hist->SetStats(0);
+ hist->Sumw2();
+ hist->Draw("HIST");
+ gPad->SetGridx();
+ gPad->SetGridy();
+
+ TFile::Open(mcFile);
+ mcHist = ((TH2*) gFile->Get("fEtaPhi"))->ProjectionY("mcHist", 41, binEnd);
+ mcHist->Rebin(2);
+ mcHist->SetLineColor(2);
+ mcHist->Scale(hist->Integral() / mcHist->Integral());
+ mcHist->Draw("SAME");
+
+ add = 0;
+
+ Float_t okRangeBegin[] = { 0.04, 0.67, 1.34 };
+ Float_t okRangeEnd[] = { 0.55, 1.24, 1.63 };
+ Float_t gapRangeBegin[] = { 0.6, 1.27 };
+ Float_t gapRangeEnd[] = { 0.65, 1.32 };
+ Float_t averageOK = FitAverage(hist, 3, okRangeBegin, okRangeEnd, 1, bins);
+ Float_t averageGap = FitAverage(hist, 2, gapRangeBegin, gapRangeEnd, 2, bins);
+ add += bins * (averageOK - averageGap);
+ Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
+
+ Float_t okRangeBegin2[] = { 2.32, 2.65 };
+ Float_t okRangeEnd2[] = { 2.55, 3.2 };
+ Float_t gapRangeBegin2[] = { 2.59, 3.3 };
+ Float_t gapRangeEnd2[] = { 2.61, 3.3 };
+ averageOK = FitAverage(hist, 2, okRangeBegin2, okRangeEnd2, 1, bins);
+ averageGap = FitAverage(hist, 2, gapRangeBegin2, gapRangeEnd2, 2, bins);
+ add += bins * (averageOK - averageGap);
+ Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
+
+ Float_t okRangeBegin3[] = { 3.55, 3.9 };
+ Float_t okRangeEnd3[] = { 3.8, 4.15 };
+ Float_t gapRangeBegin3[] = { 3.83 };
+ Float_t gapRangeEnd3[] = { 3.86 };
+ averageOK = FitAverage(hist, 2, okRangeBegin3, okRangeEnd3, 1, bins);
+ averageGap = FitAverage(hist, 1, gapRangeBegin3, gapRangeEnd3, 2, bins);
+ add += bins * (averageOK - averageGap);
+ Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
+
+ Float_t okRangeBegin4[] = { 4.2, 4.5 };
+ Float_t okRangeEnd4[] = { 4.4, 4.7 };
+ Float_t gapRangeBegin4[] = { 4.45 };
+ Float_t gapRangeEnd4[] = { 4.45 };
+ averageOK = FitAverage(hist, 2, okRangeBegin4, okRangeEnd4, 1, bins);
+ averageGap = FitAverage(hist, 1, gapRangeBegin4, gapRangeEnd4, 2, bins);
+ add += bins * (averageOK - averageGap);
+ Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
+
+ Printf("This adds %.2f %% to the total number of tracklets (%f)", 100.0 * add / hist->Integral(), hist->Integral());
+ c->SaveAs("gap2.png");
+
+ Printf("In total we add %.2f %%.", 100.0 * (add1 + add) / (total1 + hist->Integral()));
+}
{
Float_t etaMax = 1.05;
if (spd)
- etaMax = 1.79;
+ etaMax = 1.49;
TH1F* hRatios[100];
for(Int_t i=0; i<nChanges; i++) {
hRatios[i]->SetMarkerSize(0.8);
Float_t average = hRatios[i]->Integral(hRatios[i]->FindBin(-1), hRatios[i]->FindBin(1)) / (hRatios[i]->FindBin(1) - hRatios[i]->FindBin(-1) + 1);
- Printf("%s: %.2f %%" , hRatios[i]->GetTitle(), (average - 1) * 100);
+ Printf("%s %s: %.2f %%" , changes[i], hRatios[i]->GetTitle(), (average - 1) * 100);
}
TPad* p = DrawCanvasAndPad("syst_changeInXsection",700,400);
TFile::Open(fileName);
// const Char_t* changes[] = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdmoreddless", "sdlessddmore", "ddmore25","ddless25","sdmore25","sdless25", "dmore25", "dless25", "sdmoreddless25", "sdlessddmore25" };
- const Char_t* changes[] = { "ua5","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdlessddmore", "sdmoreddless" };
+ const Char_t* changes[] = { "default","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdlessddmore", "sdmoreddless" };
//const Char_t* changes[] = { "pythia", "qgsm", "phojet" };
//const Int_t nChanges = 3;
Int_t colors[] = {1,1,4,1,2,2,4,2,1};
{
TFile::Open(fileName);
- const Char_t* changes[] = { "PythiaRatios", "KBoosted", "KReduced", "pBoosted", "pReduced", "KBoostedpBoosted", "KReducedpReduced", "KBoostedpReduced", "KReducedpBoosted"};
- const char* names[] = { "", "K #times 1.5", "K #times 0.5", "p #times 1.5", "p #times 0.5", "K #times 1.5, p #times 1.5", "K #times 0.5, p #times 0.5", "K #times 1.5, p #times 0.5", "K #times 0.5, p #times 1.5" };
+ const Char_t* changes[] = { "KBoosted", "KReduced", "pBoosted", "pReduced", "KBoostedpBoosted", "KReducedpReduced", "KBoostedpReduced", "KReducedpBoosted", "othersBoosted", "othersReduced" };
+ const char* names[] = { "K #times 1.3", "K #times 0.7", "p #times 1.3", "p #times 0.7", "K #times 1.3, p #times 1.3", "K #times 0.7, p #times 0.7", "K #times 1.3, p #times 0.7", "K #times 0.7, p #times 1.3", "O #times 1.3", "O #times 0.7" };
//const Char_t* changes[] = { "PythiaRatios", "PiBoosted", "PiReduced", "KBoosted", "KReduced", "pBoosted", "pReduced", "othersBoosted", "othersReduced" };
//const char* names[] = { "", "#pi #times 1.5", "#pi #times 0.5", "K #times 1.5", "K #times 0.5", "p #times 1.5", "p #times 0.5", "others #times 1.5", "others #times 0.5" };
- Int_t colors[] = {1,1,2,2,1,2,1,4,4};
+ Int_t colors[] = {1,1,2,2,1,2,1,4,4,1,1};
- c = DrawChange(spd, "", changes, 9, 9, colors, names, 0.03);
- c->SaveAs("compositions.eps");
+ c = DrawChange(spd, "", changes, 10, 10, colors, names, 0.03);
+ //c->SaveAs("compositions.eps");
}
TPad* DrawCanvasAndPad(const Char_t* name, Int_t sizeX=600, Int_t sizeY=500) {
void Load(const char* taskName, Bool_t debug)
{
TString compileTaskName;
- compileTaskName.Form("%s.cxx++", taskName);
+ compileTaskName.Form("%s.cxx+", taskName);
if (debug)
compileTaskName += "g";
if (aProof)
{
- TProof::Mgr("alicecaf")->SetROOTVersion("v5-24-00a");
TProof::Open("alicecaf");
- //gProof->SetParallel(2);
- //gProof->SetParameter("PROOF_Packetizer", "TPacketizer");
Bool_t fullAliroot = kFALSE;
// Enable the needed package
else
{
// needed if ITS recpoints are accessed, see AlidNdEtaTask, FULLALIROOT define statement
- gProof->UploadPackage("$ALICE_ROOT/v4-18-12-AN-all.par");
- gProof->EnablePackage("v4-18-12-AN-all");
+ gProof->UploadPackage("$ALICE_ROOT/v4-18-15-AN-all.par");
+ gProof->EnablePackage("v4-18-15-AN-all");
gProof->Exec("TGrid::Connect(\"alien://\")", kTRUE);
if (fullAliroot)
AliESDInputHandler* esdH = new AliESDInputHandlerRP; // for RecPoints
else
- AliESDInputHandler* esdH = new AliESDInputHandlerRP;
+ AliESDInputHandler* esdH = new AliESDInputHandler;
- esdH->SetInactiveBranches("AliESDACORDE FMD AliRawDataErrorLogs CaloClusters Cascades EMCALCells EMCALTrigger ESDfriend Kinks Kinks Cascades ALIESDACORDE MuonTracks TrdTracks CaloClusters");
+ esdH->SetInactiveBranches("FMD AliRawDataErrorLogs CaloClusters Cascades EMCALCells EMCALTrigger ESDfriend Kinks MuonTracks TrdTracks");
mgr->SetInputEventHandler(esdH);
AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kSPD | AliPWG0Helper::kFieldOn;
- AliTriggerAnalysis::Trigger trigger = AliTriggerAnalysis::kSPDGFOBits | AliTriggerAnalysis::kOfflineFlag; // AcceptAll;
+ //AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kSPD | AliPWG0Helper::kFieldOn | AliPWG0Helper::kSPDOnlyL0;
+ //AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kTPCITS | AliPWG0Helper::kFieldOn;
+
+ AliTriggerAnalysis::Trigger trigger = AliTriggerAnalysis::kAcceptAll | AliTriggerAnalysis::kOfflineFlag;
+ //AliTriggerAnalysis::Trigger trigger = AliTriggerAnalysis::kAcceptAll | AliTriggerAnalysis::kOfflineFlag | AliTriggerAnalysis::kOneParticle;
+
+ //AliTriggerAnalysis::Trigger trigger = AliTriggerAnalysis::kSPDGFOBits | AliTriggerAnalysis::kOfflineFlag;
+ //AliTriggerAnalysis::Trigger trigger = AliTriggerAnalysis::kSPDGFOBits | AliTriggerAnalysis::kOfflineFlag | AliTriggerAnalysis::kOneParticle;
+
+ //AliTriggerAnalysis::Trigger trigger = AliTriggerAnalysis::kV0AND | AliTriggerAnalysis::kOfflineFlag;
+
+ //AliTriggerAnalysis::Trigger trigger = AliTriggerAnalysis::kV0OR | AliTriggerAnalysis::kOfflineFlag;
+ //AliTriggerAnalysis::Trigger trigger = AliTriggerAnalysis::kV0OR | AliTriggerAnalysis::kOfflineFlag | AliTriggerAnalysis::kOneParticle;
- AliPWG0Helper::PrintConf(analysisMode, trigger);
+ AliPWG0Helper::DiffTreatment diffTreatment = AliPWG0Helper::kMCFlags;
+ //AliPWG0Helper::DiffTreatment diffTreatment = AliPWG0Helper::kE710Cuts;
+
+ AliPWG0Helper::PrintConf(analysisMode, trigger, diffTreatment);
AliESDtrackCuts* esdTrackCuts = 0;
if (!(analysisMode & AliPWG0Helper::kSPD))
save = kTRUE;
}
+ // physics selection
+ gROOT->ProcessLine(".L $ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+ physicsSelectionTask = AddTaskPhysicsSelection((requiredData == 2) ? kFALSE : kTRUE);
+
+ // 900 GeV
+ if (0 && requiredData == 2)
+ {
+ physicsSelectionTask->GetPhysicsSelection()->AddCollisionTriggerClass("+CINT1B-ABCE-NOPF-ALL #769 #3119");
+ physicsSelectionTask->GetPhysicsSelection()->AddBGTriggerClass("+CINT1A-ABCE-NOPF-ALL #446 #2554");
+ physicsSelectionTask->GetPhysicsSelection()->AddBGTriggerClass("+CINT1C-ABCE-NOPF-ALL #1334 #2228");
+ physicsSelectionTask->GetPhysicsSelection()->AddBGTriggerClass("+CINT1-E-NOPF-ALL #790");
+ }
+
+ // 7 TeV, run 114783
+ if (0 && requiredData == 2)
+ {
+ physicsSelectionTask->GetPhysicsSelection()->AddCollisionTriggerClass("+CINT1B-ABCE-NOPF-ALL #345");
+ physicsSelectionTask->GetPhysicsSelection()->AddBGTriggerClass("+CINT1A-ABCE-NOPF-ALL #2130");
+ physicsSelectionTask->GetPhysicsSelection()->AddBGTriggerClass("+CINT1C-ABCE-NOPF-ALL #3018");
+ physicsSelectionTask->GetPhysicsSelection()->AddBGTriggerClass("+CINT1-E-NOPF-ALL #1238");
+ }
+
+ // 7 TeV, run 114786,98
+ if (0 && requiredData == 2)
+ {
+ physicsSelectionTask->GetPhysicsSelection()->AddCollisionTriggerClass("+CINT1B-ABCE-NOPF-ALL #346");
+ physicsSelectionTask->GetPhysicsSelection()->AddBGTriggerClass("+CINT1A-ABCE-NOPF-ALL #2131");
+ physicsSelectionTask->GetPhysicsSelection()->AddBGTriggerClass("+CINT1C-ABCE-NOPF-ALL #3019");
+ physicsSelectionTask->GetPhysicsSelection()->AddBGTriggerClass("+CINT1-E-NOPF-ALL #1238");
+ //physicsSelectionTask->GetPhysicsSelection()->Initialize(114786);
+ }
+
+ // FO efficiency (for MC)
+ if (1 && requiredData != 2)
+ {
+ //const char* fastORFile = "spdFOEff_run104824_52.root";
+ const char* fastORFile = "spdFOEff_run104867_92.root";
+ //const char* fastORFile = "spdFOEff_run105054_7.root";
+ //const char* fastORFile = "spdFOEff_run114931.root";
+
+ Printf("NOTE: Simulating FAST-OR efficiency on the analysis level using file %s", fastORFile);
+ TFile::Open(fastORFile);
+ spdFOEff = (TH1F*) gFile->Get("spdFOEff");
+ physicsSelectionTask->GetPhysicsSelection()->Initialize(104867);
+ physicsSelectionTask->GetPhysicsSelection()->GetTriggerAnalysis()->SetSPDGFOEfficiency(spdFOEff);
+ }
+
+ // V0 syst. study
+ if (0)
+ {
+ Printf("NOTE: Systematic study for VZERO enabled!");
+ physicsSelectionTask->GetPhysicsSelection()->Initialize(104867);
+ for (Int_t i=0; i<1; i++)
+ {
+ // for MC and data
+ //physicsSelectionTask->GetPhysicsSelection()->GetTriggerAnalysis(i)->SetV0HwPars(15, 61.5, 86.5);
+ physicsSelectionTask->GetPhysicsSelection()->GetTriggerAnalysis(i)->SetV0AdcThr(6);
+ // only for MC
+ //physicsSelectionTask->GetPhysicsSelection()->GetTriggerAnalysis(i)->SetV0HwPars(0, 0, 125);
+ //physicsSelectionTask->GetPhysicsSelection()->GetTriggerAnalysis(i)->SetV0AdcThr(0);
+ }
+ }
+
+ // BG study
+ //physicsSelectionTask->GetPhysicsSelection()->AddCollisionTriggerClass("+CINT1A-ABCE-NOPF-ALL");
+ //physicsSelectionTask->GetPhysicsSelection()->AddCollisionTriggerClass("+CINT1C-ABCE-NOPF-ALL");
+
// Create, add task
if (runWhat == 0 || runWhat == 2)
{
//task->SetOnlyPrimaries();
//task->SetFillPhi();
//task->SetSymmetrize();
-
+
+ // INEL>0 definition
+ //task->SetMultAxisEta1();
+
task->SetTrigger(trigger);
task->SetAnalysisMode(analysisMode);
task->SetTrackCuts(esdTrackCuts);
- //task->SetDeltaPhiCut(0.05);
-
- if (requiredData == 2)
- task->SetCheckEventType();
- task->SetTriggerClasses(requireClass, rejectClass);
+ //task->SetDeltaPhiCut(0.064);
+ task->SetDiffTreatment(diffTreatment);
+
+ //if (requiredData == 2)
+ // task->SetCheckEventType();
+ //task->SetTriggerClasses(requireClass, rejectClass);
mgr->AddTask(task);
cOutput = mgr->CreateContainer("cOutput", TList::Class(), AliAnalysisManager::kOutputContainer);
mgr->ConnectOutput(task, 0, cOutput);
}
+
if (runWhat == 1 || runWhat == 2)
{
Load("AlidNdEtaCorrectionTask", aDebug);
//task2->SetOnlyPrimaries();
//task2->SetSymmetrize();
+ // to account for gaps in real life SPD geometry
+ task2->SetSkipParticles();
+
+ // INEL>0 definition
+ //task2->SetMultAxisEta1();
+
task2->SetTrigger(trigger);
task2->SetAnalysisMode(analysisMode);
task2->SetTrackCuts(esdTrackCuts);
- //task2->SetDeltaPhiCut(0.05);
+ //task2->SetDeltaPhiCut(0.064);
+ task2->SetDiffTreatment(diffTreatment);
mgr->AddTask(task2);
mgr->ConnectOutput(task2, 0, cOutput);
}
- if (requiredData == 1) {
+ if (requiredData == 1)
+ {
// Enable MC event handler
AliMCEventHandler* handler = new AliMCEventHandler;
handler->SetReadTR(kFALSE);
case AliTriggerAnalysis::kMB3: path += "/mb3"; break;
case AliTriggerAnalysis::kSPDGFO: path += "/spdgfo"; break;
case AliTriggerAnalysis::kSPDGFOBits: path += "/spdgfobits"; break;
+ case AliTriggerAnalysis::kAcceptAll: path += "/all"; break;
+ case AliTriggerAnalysis::kV0AND: path += "/v0and"; break;
+ case AliTriggerAnalysis::kV0OR: path += "/v0or"; break;
+ case AliTriggerAnalysis::kNSD1: path += "/nsd1"; break;
+ case AliTriggerAnalysis::kMB1Prime: path += "/mb1prime"; break;
default: Printf("ERROR: Trigger undefined for path to files"); return;
}
+ if (trigger & AliTriggerAnalysis::kOneParticle)
+ path += "-onepart";
+
if (strlen(requireClass) > 0 && strlen(rejectClass) == 0)
{
path += Form("/%s", requireClass);
if (analysisMode & AliPWG0Helper::kSPD)
path += "/spd";
+ if (analysisMode & AliPWG0Helper::kSPDOnlyL0)
+ path += "onlyL0";
+
if (analysisMode & AliPWG0Helper::kTPC)
path += "/tpc";
+ if (analysisMode & AliPWG0Helper::kTPCITS)
+ path += "/tpcits";
+
gSystem->mkdir(path, kTRUE);
if (runWhat == 0 || runWhat == 2)
{
}
if (runWhat == 1 || runWhat == 2)
{
- gSystem->Rename("correction_map.root", path + "/correction_map.root");
+ if (optStr.Contains("process-types"))
+ gSystem->Rename("correction_mapprocess-types.root", path + "/correction_mapprocess-types.root");
+ else
+ gSystem->Rename("correction_map.root", path + "/correction_map.root");
}
+ gSystem->Rename("event_stat.root", path + "/event_stat.root");
Printf(">>>>> Moved files to %s", path.Data());
}
{
gROOT->ProcessLine(".L CreateChainFromDataSet.C");
ds = gProof->GetDataSet(data)->GetStagedSubset();
- chain = CreateChainFromDataSet(ds);
- mgr->StartAnalysis("local", chain, nRuns, offset);
+ chain = CreateChainFromDataSet(ds, "esdTree", nRuns);
+ mgr->StartAnalysis("local", chain, 1234567890, offset);
}
else
{
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff