updated dndeta analysis
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 13 Apr 2010 13:59:28 +0000 (13:59 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 13 Apr 2010 13:59:28 +0000 (13:59 +0000)
16 files changed:
PWG0/AliCorrection.cxx
PWG0/AliPWG0Helper.cxx
PWG0/AliPWG0Helper.h
PWG0/dNdEta/AlidNdEtaCorrection.cxx
PWG0/dNdEta/AlidNdEtaCorrection.h
PWG0/dNdEta/AlidNdEtaCorrectionTask.cxx
PWG0/dNdEta/AlidNdEtaCorrectionTask.h
PWG0/dNdEta/AlidNdEtaTask.cxx
PWG0/dNdEta/AlidNdEtaTask.h
PWG0/dNdEta/correct.C
PWG0/dNdEta/dNdEtaAnalysis.cxx
PWG0/dNdEta/dNdEtaAnalysis.h
PWG0/dNdEta/drawPlots.C
PWG0/dNdEta/drawSystematics.C
PWG0/dNdEta/drawSystematicsNew.C
PWG0/dNdEta/run.C

index ba411f00a128a479f0bee2dd8157709b892ddd58..bca47670ba126eb641c76dc846f6e3230b8156d7 100644 (file)
@@ -77,12 +77,15 @@ AliCorrection::AliCorrection(const Char_t* name, const Char_t* title, AliPWG0Hel
   // 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,
@@ -96,25 +99,25 @@ AliCorrection::AliCorrection(const Char_t* name, const Char_t* title, AliPWG0Hel
                             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;
@@ -439,7 +442,10 @@ void AliCorrection::PrintStats(Float_t zRange, Float_t etaRange, Float_t ptCut)
   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);
index a242ae4cd1daa4fd9f3574e7d7cca1952aba0414..4ce88d4b4ca35a26887c2513942880d2cb34a70f 100644 (file)
@@ -19,6 +19,7 @@
 #include <AliESD.h>
 #include <AliESDEvent.h>
 #include <AliESDVertex.h>
+#include <AliESDRun.h>
 #include <AliVertexerTracks.h>
 #include <AliMultiplicity.h>
 
@@ -70,7 +71,7 @@ Bool_t AliPWG0Helper::TestVertex(const AliESDVertex* vertex, AnalysisMode analys
 }
 
 //____________________________________________________________________
-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
   //
@@ -78,17 +79,29 @@ const AliESDVertex* AliPWG0Helper::GetVertex(AliESDEvent* aEsd, AnalysisMode ana
   // 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);
@@ -356,13 +369,13 @@ AliPWG0Helper::MCProcessType AliPWG0Helper::GetDPMjetEventProcessType(AliGenEven
   }
 
 
-  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;
@@ -500,7 +513,7 @@ void AliPWG0Helper::NormalizeToBinWidth(TH2* hist)
 }
 
 //____________________________________________________________________
-void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, AliTriggerAnalysis::Trigger trigger)
+void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, AliTriggerAnalysis::Trigger trigger, AliPWG0Helper::DiffTreatment diffTreatment)
 {
   //
   // Prints the given configuration
@@ -529,8 +542,189 @@ void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, AliTriggerAnalysis::Tri
   
   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;
+}
index 763bba3d486810db8cdd5408b8720b4bb442bda1..de29b5bf26dd68b8b22b23b651f7c3dd3bae7d4a 100644 (file)
@@ -19,22 +19,25 @@ class AliGenEventHeader;
 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);
@@ -47,7 +50,10 @@ class AliPWG0Helper : public TObject
     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
index 0ff81baf751f34daf99b3cc7ae894a9803994c0d..65c6c88318786ebb17a6c639c537c166715cadd7 100644 (file)
@@ -37,7 +37,8 @@ AlidNdEtaCorrection::AlidNdEtaCorrection()
   fVertexRecoCorrection(0),
   fTriggerBiasCorrectionMBToINEL(0),
   fTriggerBiasCorrectionMBToNSD(0),
-  fTriggerBiasCorrectionMBToND(0)
+  fTriggerBiasCorrectionMBToND(0),
+  fTriggerBiasCorrectionMBToOnePart(0)
 {
   // default constructor
 }
@@ -49,7 +50,8 @@ AlidNdEtaCorrection::AlidNdEtaCorrection(const Char_t* name, const Char_t* title
   fVertexRecoCorrection(0),
   fTriggerBiasCorrectionMBToINEL(0),
   fTriggerBiasCorrectionMBToNSD(0),
-  fTriggerBiasCorrectionMBToND(0)
+  fTriggerBiasCorrectionMBToND(0),
+  fTriggerBiasCorrectionMBToOnePart(0)
 {
   //
   // constructor
@@ -61,6 +63,7 @@ AlidNdEtaCorrection::AlidNdEtaCorrection(const Char_t* name, const Char_t* title
   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);
 }
 
 //____________________________________________________________________
@@ -92,6 +95,11 @@ AlidNdEtaCorrection::~AlidNdEtaCorrection()
     delete fTriggerBiasCorrectionMBToND;
     fTriggerBiasCorrectionMBToND = 0;
   }
+
+  if (fTriggerBiasCorrectionMBToOnePart) {
+    delete fTriggerBiasCorrectionMBToOnePart;
+    fTriggerBiasCorrectionMBToOnePart = 0;
+  }
 }
 
 //____________________________________________________________________
@@ -107,6 +115,7 @@ AlidNdEtaCorrection::Finish() {
   fTriggerBiasCorrectionMBToINEL->Divide();
   fTriggerBiasCorrectionMBToNSD->Divide();
   fTriggerBiasCorrectionMBToND->Divide();
+  fTriggerBiasCorrectionMBToOnePart->Divide();
 }
 
 //____________________________________________________________________
@@ -131,6 +140,7 @@ Long64_t AlidNdEtaCorrection::Merge(TCollection* list)
   TList* collectionTriggerBiasMBToINEL  = new TList;
   TList* collectionTriggerBiasMBToNSD   = new TList;
   TList* collectionTriggerBiasMBToND    = new TList;
+  TList* collectionTriggerBiasMBToOnePart    = new TList;
 
   Int_t count = 0;
   while ((obj = iter->Next())) {
@@ -144,6 +154,7 @@ Long64_t AlidNdEtaCorrection::Merge(TCollection* list)
     collectionTriggerBiasMBToINEL->Add(entry->fTriggerBiasCorrectionMBToINEL);
     collectionTriggerBiasMBToNSD ->Add(entry->fTriggerBiasCorrectionMBToNSD);
     collectionTriggerBiasMBToND  ->Add(entry->fTriggerBiasCorrectionMBToND);
+    collectionTriggerBiasMBToOnePart  ->Add(entry->fTriggerBiasCorrectionMBToOnePart);
 
     count++;
   }
@@ -152,12 +163,14 @@ Long64_t AlidNdEtaCorrection::Merge(TCollection* list)
   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;
 }
@@ -173,7 +186,22 @@ void AlidNdEtaCorrection::Add(AlidNdEtaCorrection* aCorrectionsToAdd, Float_t c)
   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);
 }
 
 //____________________________________________________________________
@@ -187,7 +215,7 @@ void AlidNdEtaCorrection::Reset(void) {
   fTriggerBiasCorrectionMBToINEL ->Reset();
   fTriggerBiasCorrectionMBToNSD  ->Reset();
   fTriggerBiasCorrectionMBToND   ->Reset();
-    
+  fTriggerBiasCorrectionMBToOnePart   ->Reset();
 }
 
 
@@ -211,6 +239,7 @@ Bool_t AlidNdEtaCorrection::LoadHistograms(const Char_t* dir)
   fTriggerBiasCorrectionMBToINEL ->LoadHistograms();
   fTriggerBiasCorrectionMBToNSD  ->LoadHistograms();
   fTriggerBiasCorrectionMBToND   ->LoadHistograms();
+  fTriggerBiasCorrectionMBToOnePart   ->LoadHistograms();
 
   gDirectory->cd("..");
 
@@ -232,6 +261,7 @@ void AlidNdEtaCorrection::SaveHistograms()
   fTriggerBiasCorrectionMBToINEL->SaveHistograms();
   fTriggerBiasCorrectionMBToNSD ->SaveHistograms();
   fTriggerBiasCorrectionMBToND  ->SaveHistograms();
+  fTriggerBiasCorrectionMBToOnePart->SaveHistograms();
 
   gDirectory->cd("..");
 }
@@ -248,6 +278,7 @@ void AlidNdEtaCorrection::DrawHistograms()
   fTriggerBiasCorrectionMBToINEL->DrawHistograms();
   fTriggerBiasCorrectionMBToNSD ->DrawHistograms();
   fTriggerBiasCorrectionMBToND  ->DrawHistograms();
+  fTriggerBiasCorrectionMBToOnePart  ->DrawHistograms();
 }
 
 //____________________________________________________________________
@@ -262,6 +293,7 @@ void AlidNdEtaCorrection::DrawOverview(const char* canvasName)
   fTriggerBiasCorrectionMBToINEL->DrawOverview(canvasName);
   fTriggerBiasCorrectionMBToNSD ->DrawOverview(canvasName);
   fTriggerBiasCorrectionMBToND  ->DrawOverview(canvasName);
+  fTriggerBiasCorrectionMBToOnePart  ->DrawOverview(canvasName);
 }
 
 //____________________________________________________________________
@@ -272,18 +304,22 @@ void AlidNdEtaCorrection::FillMCParticle(Float_t vtx, Float_t eta, Float_t pt, B
 
   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)
@@ -309,18 +345,22 @@ void AlidNdEtaCorrection::FillEvent(Float_t vtx, Float_t n, Bool_t trigger, Bool
 
   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)
@@ -422,6 +462,7 @@ void AlidNdEtaCorrection::ReduceInformation()
   fTriggerBiasCorrectionMBToINEL ->ReduceInformation();
   fTriggerBiasCorrectionMBToNSD  ->ReduceInformation();
   fTriggerBiasCorrectionMBToND   ->ReduceInformation();
+  fTriggerBiasCorrectionMBToOnePart   ->ReduceInformation();
 }
 
 //____________________________________________________________________
@@ -437,6 +478,7 @@ AliCorrection* AlidNdEtaCorrection::GetCorrection(CorrectionType correctionType)
     case kINEL :           return fTriggerBiasCorrectionMBToINEL;
     case kNSD :            return fTriggerBiasCorrectionMBToNSD;
     case kND :             return fTriggerBiasCorrectionMBToND;
+    case kOnePart :             return fTriggerBiasCorrectionMBToOnePart;
   }
 
   return 0;
index dbefa84ed5275e5b64bc41bdfbad00a3902b6ef4..c94978bbcd659845643f26a7b5d0f7a9c4464cd4 100644 (file)
@@ -31,7 +31,8 @@ public:
     kVertexReco,      // MB sample
     kINEL,
     kNSD,
-    kND
+    kND,
+    kOnePart
   };
 
   AlidNdEtaCorrection();
@@ -52,10 +53,12 @@ public:
   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);
@@ -73,12 +76,13 @@ protected:
   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
index 2dde4d98e56fad128b981a3e948d2c0324407209..0593bec3b4d544ec4301f9d16349406a69bdcbda 100644 (file)
 #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>
@@ -30,6 +32,7 @@
 #include "dNdEta/dNdEtaAnalysis.h"
 #include "dNdEta/AlidNdEtaCorrection.h"
 #include "AliTriggerAnalysis.h"
+#include "AliPhysicsSelection.h"
 
 ClassImp(AlidNdEtaCorrectionTask)
 
@@ -43,9 +46,12 @@ AlidNdEtaCorrectionTask::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),
@@ -81,6 +87,8 @@ AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask() :
   
   for (Int_t i=0; i<8; i++)
     fDeltaPhi[i] = 0;
+
+  AliLog::SetClassDebugLevel("AlidNdEtaCorrectionTask", AliLog::kWarning);
 }
 
 AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
@@ -93,9 +101,12 @@ 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),
@@ -135,6 +146,8 @@ AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
   
   for (Int_t i=0; i<8; i++)
     fDeltaPhi[i] = 0;
+
+  AliLog::SetClassDebugLevel("AlidNdEtaCorrectionTask", AliLog::kWarning);
 }
 
 AlidNdEtaCorrectionTask::~AlidNdEtaCorrectionTask()
@@ -255,12 +268,12 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
   }
 
   
-  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);
@@ -270,7 +283,7 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
   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);
@@ -298,15 +311,18 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
     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");
@@ -317,7 +333,8 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
   if (fEsdTrackCuts)
   {
     fEsdTrackCuts->SetName("fEsdTrackCuts");
-    fOutput->Add(fEsdTrackCuts);
+    // TODO like this we send an empty object back...
+    fOutput->Add(fEsdTrackCuts->Clone());
   }
 }
 
@@ -338,11 +355,26 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
   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");
 
@@ -376,12 +408,16 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
     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();
@@ -394,27 +430,19 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
   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;
@@ -424,48 +452,81 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
   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)
@@ -476,6 +537,8 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       return;
     }
     
+    Bool_t foundInEta10 = kFALSE;
+    
     if (vtxESD)
     {
       // get multiplicity from ESD tracks
@@ -500,7 +563,8 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
           continue;
         }
         
-        //Printf("status is: %u", esdTrack->GetStatus());
+        if (esdTrack->Pt() < 0.15)
+          continue;
         
         if (fOnlyPrimaries)
         {
@@ -510,7 +574,11 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
           
           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)
@@ -524,7 +592,8 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
 
       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++)
@@ -584,14 +653,59 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
         }
       }
     }
+    
+    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);
@@ -609,7 +723,7 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
 
     if (SignOK(particle->GetPDG()) == kFALSE)
       continue;
-
+    
     if (fPIDParticles && TMath::Abs(particle->Eta()) < 1.0)
       fPIDParticles->Fill(particle->GetPdgCode());
 
@@ -635,21 +749,24 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
     //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"))
@@ -662,7 +779,7 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
         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)
@@ -674,7 +791,8 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       nAccepted++;
   }
 
-  fEventStats->Fill(AliPWG0Helper::GetLastProcessType(), biny);
+  if (AliPWG0Helper::GetLastProcessType() >= -1)
+    fEventStats->Fill(AliPWG0Helper::GetLastProcessType(), biny);
 
   fMultAll->Fill(nAccepted);
   if (eventTriggered) {
@@ -694,11 +812,17 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
   }
 
   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];
 
@@ -715,7 +839,8 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       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
@@ -784,7 +909,7 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       {
         if (fAnalysisMode & AliPWG0Helper::kSPD && !fFillPhi)
         {
-          thirdDim = inputCount;
+          thirdDim = (fMultAxisEta1) ? nEta10 : inputCount;
         }
         else
           thirdDim = thirdDimArr[i];
@@ -816,8 +941,11 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       }
 
       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)
       {
@@ -936,43 +1064,50 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
   {
     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;
index ba26c410b3e9636e3fdd0a54cc3262f7d02e03e8..8a4946b6b8721b3652f2752fa900fe77b55556e3 100644 (file)
@@ -36,6 +36,9 @@ class AlidNdEtaCorrectionTask : public AliAnalysisTask {
     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; }
 
@@ -51,10 +54,13 @@ class AlidNdEtaCorrectionTask : public AliAnalysisTask {
     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
 
index 08c56501a69871dbc188607c14ddb04b67ad42fb..ae33c25ff276fa7362060f46f340c0a4d03a2244 100644 (file)
@@ -69,8 +69,9 @@ AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
   fUseMCKine(kFALSE),
   fCheckEventType(kFALSE),
   fSymmetrize(kFALSE),
+  fMultAxisEta1(kFALSE),
+  fDiffTreatment(AliPWG0Helper::kMCFlags),
   fEsdTrackCuts(0),
-  fPhysicsSelection(0),
   fTriggerAnalysis(0),
   fdNdEtaAnalysisESD(0),
   fMult(0),
@@ -80,6 +81,7 @@ AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
   fdNdEtaAnalysis(0),
   fdNdEtaAnalysisND(0),
   fdNdEtaAnalysisNSD(0),
+  fdNdEtaAnalysisOnePart(0),
   fdNdEtaAnalysisTr(0),
   fdNdEtaAnalysisTrVtx(0),
   fdNdEtaAnalysisTracks(0),
@@ -95,7 +97,6 @@ AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
   fFiredChips(0),
   fTrackletsVsClusters(0),
   fTrackletsVsUnassigned(0),
-  fTriggerVsTime(0),
   fStats(0),
   fStats2(0)
 {
@@ -211,7 +212,7 @@ void AlidNdEtaTask::CreateOutputObjects()
 
   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]);
   }
@@ -219,22 +220,15 @@ void AlidNdEtaTask::CreateOutputObjects()
   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");
@@ -246,11 +240,11 @@ void AlidNdEtaTask::CreateOutputObjects()
   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");
@@ -270,9 +264,9 @@ void AlidNdEtaTask::CreateOutputObjects()
   
   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);
@@ -311,6 +305,9 @@ void AlidNdEtaTask::CreateOutputObjects()
     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);
 
@@ -331,17 +328,7 @@ void AlidNdEtaTask::CreateOutputObjects()
     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*)
@@ -365,37 +352,24 @@ 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);
@@ -413,134 +387,102 @@ void AlidNdEtaTask::Exec(Option_t*)
       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)
         {
@@ -548,6 +490,15 @@ void AlidNdEtaTask::Exec(Option_t*)
           filled = kTRUE;
         }
       }
+      
+      if (TMath::Abs(vtx[2]) > 10)
+      {
+        if (!filled)
+        {
+          fStats2->Fill(4, vZero);
+          filled = kTRUE;
+        }
+      }
         
       const AliMultiplicity* mult = fESD->GetMultiplicity();
       if (!mult)
@@ -557,26 +508,17 @@ void AlidNdEtaTask::Exec(Option_t*)
       {
         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);
@@ -586,14 +528,16 @@ void AlidNdEtaTask::Exec(Option_t*)
     // 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());
@@ -607,7 +551,7 @@ void AlidNdEtaTask::Exec(Option_t*)
           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());
             }
@@ -616,10 +560,13 @@ void AlidNdEtaTask::Exec(Option_t*)
               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;
@@ -673,7 +620,7 @@ void AlidNdEtaTask::Exec(Option_t*)
         return;
       }
     }
-
+    
     // create list of (label, eta, pt) tuples
     Int_t inputCount = 0;
     Int_t* labelArr = 0;
@@ -681,99 +628,114 @@ void AlidNdEtaTask::Exec(Option_t*)
     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)
     {
@@ -783,6 +745,8 @@ void AlidNdEtaTask::Exec(Option_t*)
         return;
       }
 
+      Bool_t foundInEta10 = kFALSE;
+      
       if (vtxESD)
       {
         // get multiplicity from ESD tracks
@@ -821,6 +785,13 @@ void AlidNdEtaTask::Exec(Option_t*)
           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)
@@ -850,13 +821,16 @@ void AlidNdEtaTask::Exec(Option_t*)
           ++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;
@@ -868,8 +842,6 @@ void AlidNdEtaTask::Exec(Option_t*)
       fMult->Fill(inputCount);
       fdNdEtaAnalysisESD->FillTriggeredEvent(inputCount);
 
-      fTriggerVsTime->SetPoint(fTriggerVsTime->GetN(), fESD->GetTimeStamp(), 1);
-
       if (vtxESD)
       {
         // control hist
@@ -877,8 +849,8 @@ void AlidNdEtaTask::Exec(Option_t*)
         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];
@@ -895,13 +867,26 @@ void AlidNdEtaTask::Exec(Option_t*)
             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)
@@ -931,7 +916,7 @@ void AlidNdEtaTask::Exec(Option_t*)
                 thirdDim = particle->Phi();
               }
               else
-                thirdDim = inputCount;
+                thirdDim = multAxis;
             }
             else
               thirdDim = particle->Pt();
@@ -943,11 +928,10 @@ void AlidNdEtaTask::Exec(Option_t*)
           } // 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)
@@ -994,18 +978,19 @@ void AlidNdEtaTask::Exec(Option_t*)
 
     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)
@@ -1025,7 +1010,10 @@ void AlidNdEtaTask::Exec(Option_t*)
       //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++;
     }
@@ -1069,6 +1057,9 @@ void AlidNdEtaTask::Exec(Option_t*)
 
       if (processType == AliPWG0Helper::kND)
         fdNdEtaAnalysisND->FillTrack(vtxMC[2], eta, thirdDim);
+      
+      if (oneParticleEvent)
+        fdNdEtaAnalysisOnePart->FillTrack(vtxMC[2], eta, thirdDim);
 
       if (eventTriggered)
       {
@@ -1090,6 +1081,8 @@ void AlidNdEtaTask::Exec(Option_t*)
       fdNdEtaAnalysisNSD->FillEvent(vtxMC[2], nAcceptedParticles);
     if (processType == AliPWG0Helper::kND)
       fdNdEtaAnalysisND->FillEvent(vtxMC[2], nAcceptedParticles);
+    if (oneParticleEvent)
+      fdNdEtaAnalysisOnePart->FillEvent(vtxMC[2], nAcceptedParticles);
 
     if (eventTriggered)
     {
@@ -1118,7 +1111,7 @@ void AlidNdEtaTask::Terminate(Option_t *)
     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"));
@@ -1133,13 +1126,10 @@ void AlidNdEtaTask::Terminate(Option_t *)
     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)
@@ -1164,10 +1154,11 @@ void AlidNdEtaTask::Terminate(Option_t *)
 
     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)
       {
@@ -1207,15 +1198,6 @@ void AlidNdEtaTask::Terminate(Option_t *)
     if (fEsdTrackCuts)
       fEsdTrackCuts->SaveHistograms("esd_track_cuts");
 
-    if (fPhysicsSelection)
-    {
-      fPhysicsSelection->SaveHistograms("physics_selection");
-      fPhysicsSelection->Print();
-    }
-    
-    if (fTriggerAnalysis)
-      fTriggerAnalysis->SaveHistograms();
-
     if (fMult)
       fMult->Write();
 
@@ -1230,7 +1212,11 @@ void AlidNdEtaTask::Terminate(Option_t *)
       fEvents->Write();
 
     if (fVertexResolution)
+    {
       fVertexResolution->Write();
+      fVertexResolution->ProjectionX();
+      fVertexResolution->ProjectionY();
+    }
 
     if (fDeltaPhi)
       fDeltaPhi->Write();
@@ -1263,9 +1249,6 @@ void AlidNdEtaTask::Terminate(Option_t *)
     if (fTrackletsVsUnassigned)
       fTrackletsVsUnassigned->Write();
     
-    if (fTriggerVsTime)
-      fTriggerVsTime->Write();
-
     if (fStats)
       fStats->Write();
 
@@ -1291,6 +1274,7 @@ void AlidNdEtaTask::Terminate(Option_t *)
       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"));
@@ -1306,6 +1290,7 @@ void AlidNdEtaTask::Terminate(Option_t *)
     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);
@@ -1319,6 +1304,7 @@ void AlidNdEtaTask::Terminate(Option_t *)
     fdNdEtaAnalysis->SaveHistograms();
     fdNdEtaAnalysisND->SaveHistograms();
     fdNdEtaAnalysisNSD->SaveHistograms();
+    fdNdEtaAnalysisOnePart->SaveHistograms();
     fdNdEtaAnalysisTr->SaveHistograms();
     fdNdEtaAnalysisTrVtx->SaveHistograms();
     fdNdEtaAnalysisTracks->SaveHistograms();
index 2a42b2e373f3c942b63a62427cf05dc9d950e21d..66298a55741f711b3fa42997d4e0a487baa76911 100644 (file)
@@ -14,8 +14,6 @@ class TH1F;
 class TH2F;
 class TH3F;
 class AliESDEvent;
-class TGraph;
-class AliPhysicsSelection;
 class AliTriggerAnalysis;
 
 class AlidNdEtaTask : public AliAnalysisTask {
@@ -41,6 +39,8 @@ 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; }
 
@@ -62,24 +62,26 @@ class AlidNdEtaTask : public AliAnalysisTask {
     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)
@@ -100,7 +102,6 @@ class AlidNdEtaTask : public AliAnalysisTask {
     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
 
index 6264bb5d7ba8ad47d43d6965b73c7a32a5749a4e..058cf5f8b808bd89b6d094215a129f5cd772904d 100644 (file)
@@ -12,11 +12,16 @@ void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const ch
 {
   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)
@@ -24,19 +29,47 @@ void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const ch
     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();
@@ -44,7 +77,7 @@ void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const ch
   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();
@@ -52,7 +85,7 @@ void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const ch
   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();
@@ -60,7 +93,7 @@ void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const ch
   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();
@@ -68,7 +101,15 @@ void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const ch
   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();
@@ -92,13 +133,15 @@ void* FinishAnalysis(const char* analysisFile = "analysis_esd_raw.root", const c
     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);
@@ -115,9 +158,27 @@ void* FinishAnalysis(const char* analysisFile = "analysis_esd_raw.root", const c
   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);
 }
index e4584cdb297090a76631f74f48f6b48481b7d93b..ebd0793ec74d6a151825cb9602ba648fd7e7577c 100644 (file)
@@ -194,7 +194,7 @@ void dNdEtaAnalysis::FillTriggeredEvent(Float_t n)
 }
 
 //____________________________________________________________________
-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
@@ -204,6 +204,9 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
 
   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();
@@ -246,6 +249,12 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
         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;
     }
   }
@@ -256,7 +265,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
 
   fData->ResetErrorsOnCorrections();
   fData->Multiply();
-
+  
   if (correctionType >= AlidNdEtaCorrection::kVertexReco)
   {
     // There are no events with vertex that have 0 multiplicity, therefore
@@ -269,9 +278,10 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
     //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");
 
@@ -283,10 +293,16 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
     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++)
     {
@@ -296,22 +312,35 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
       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);
     }
@@ -396,7 +425,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
   }
 
   //new TCanvas(tag, tag, 800, 600);  vtxVsEta->DrawCopy("COLZ");
-
+  
   // clear result histograms
   for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
   {
@@ -416,7 +445,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
       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));
@@ -425,20 +454,22 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
       // 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;
@@ -529,6 +560,15 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
         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);
 
@@ -541,7 +581,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
         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);
       }
     }
   }
index 6eea884152e6f653b4385e10359404c6826af8c0..0d317ecbdd4756b9700c74fc14cb6bf7609c8e5f 100644 (file)
@@ -44,7 +44,7 @@ public:
   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);
index e91d270cf9ae8a8e691715b8947408d8d66fd86f..d574d825f948e563d50e381e0a9fb1620b75e32e 100644 (file)
@@ -442,8 +442,11 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
   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");
@@ -452,7 +455,10 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
   TH1* histESDMBTracksNoPt = (TH1*) file->Get("dndetaTracks/dNdEta");
 
   Prepare1DPlot(histESD);
+  Prepare1DPlot(histESD1);
+  Prepare1DPlot(histESD2);
   Prepare1DPlot(histESDnsd);
+  Prepare1DPlot(histESDonePart);
   Prepare1DPlot(histESDMB);
   Prepare1DPlot(histESDMBVtx);
 
@@ -463,6 +469,7 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
 
   histESD->SetLineWidth(0);
   histESDnsd->SetLineWidth(0);
+  histESDonePart->SetLineWidth(0);
   histESDMB->SetLineWidth(0);
   histESDMBVtx->SetLineWidth(0);
 
@@ -472,11 +479,13 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
 
   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);
 
@@ -487,6 +496,7 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
 
   histESD->SetMarkerStyle(20);
   histESDnsd->SetMarkerStyle(29);
+  histESDonePart->SetMarkerStyle(24);
   histESDMB->SetMarkerStyle(21);
   histESDMBVtx->SetMarkerStyle(22);
 
@@ -495,7 +505,7 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
   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;
@@ -504,6 +514,7 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
   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);
@@ -520,6 +531,7 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
   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);
@@ -536,17 +548,91 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
   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;
 
@@ -554,18 +640,19 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
 
   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);
@@ -633,6 +720,7 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
   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);
@@ -713,7 +801,7 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
   ratioNSD->Draw("HIST EP SAME");
   legend7->Draw();
   
-  c7->SaveAs("ratios.eps");
+  //c7->SaveAs("ratios.eps");
 
   new TCanvas;
   dummy2->DrawCopy();
@@ -740,6 +828,7 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
 
   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)");
@@ -834,6 +923,165 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
   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);
@@ -1036,8 +1284,8 @@ void TriggerBias1D(const char* fileName = "correction_map.root", const char* fol
   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);
@@ -1065,10 +1313,21 @@ void TriggerBias1D(const char* fileName = "correction_map.root", const char* fol
 
   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()
@@ -1524,10 +1783,10 @@ void CompareTrack2Particle1D(const char* file1, const char* file2, Float_t upper
   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();
@@ -1556,7 +1815,7 @@ TCanvas* Track2Particle2D(const char* fileName = "correction_map.root", const ch
 {
   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"));
@@ -1873,12 +2132,12 @@ void CompareCorrection2Measured(Float_t ptMin = 0.301, const char* dataInput = "
   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");
@@ -2854,13 +3113,34 @@ void PrintEventStats(Int_t corrID = 3, const char* fileName = "correction_map.ro
   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 +++");
   
@@ -2874,10 +3154,6 @@ void PrintEventStats(Int_t corrID = 3, const char* fileName = "correction_map.ro
   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);
@@ -3150,6 +3426,34 @@ void FitDiamond()
   }
 }
 
+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");
@@ -3198,10 +3502,17 @@ void CompareQualityHists(const char* fileName1, const char* fileName2, const cha
   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);
@@ -3218,31 +3529,63 @@ void CompareQualityHists(TH1* hist1, TH1* hist2, Int_t rebin1 = 1, Int_t 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);
@@ -3250,6 +3593,7 @@ void CompareQualityHists(TH1* hist1, TH1* hist2, Int_t rebin1 = 1, Int_t rebin2
     new TCanvas;
     hist1->SetMarkerStyle(20);
     hist1->DrawCopy("P");
+    */
   }
 }
 
@@ -3274,7 +3618,7 @@ void DrawClustersVsTracklets()
   
   c->SaveAs("clusters_vs_tracklets.eps");
 }
-  
+
 void VertexPlotBackgroundNote()
 {
   TFile::Open("all.root");
@@ -3292,8 +3636,6 @@ void VertexPlotBackgroundNote()
   
   proj->SetLineColor(2);
   proj->Draw("SAME");
-  
-  
 }
 
 void BackgroundAnalysis(const char* signal, const char* background)
@@ -3373,3 +3715,112 @@ void DrawStats(Bool_t all = kFALSE)
     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
index 38de7a779d093596dc94d83b6af2d352211bc365..f2870e0c4efcecaf94e698a78cfd705cfa337961 100644 (file)
@@ -504,34 +504,42 @@ const char* ChangeComposition(void** correctionPointer, Int_t index)
       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;
@@ -579,7 +587,13 @@ void Composition()
   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)
   {
@@ -621,7 +635,7 @@ void Composition()
 
     // 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);
@@ -776,7 +790,7 @@ void DrawdNdEtaDifferences()
   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
@@ -789,6 +803,17 @@ void mergeCorrectionsWithDifferentCrosssections(Int_t correctionTarget = 3, Char
   // 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();
 
@@ -805,16 +830,133 @@ void mergeCorrectionsWithDifferentCrosssections(Int_t correctionTarget = 3, Char
   //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};
@@ -825,9 +967,9 @@ void mergeCorrectionsWithDifferentCrosssections(Int_t correctionTarget = 3, Char
   
   // 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);
@@ -840,6 +982,7 @@ void mergeCorrectionsWithDifferentCrosssections(Int_t correctionTarget = 3, Char
   correctionStandard->GetTriggerBiasCorrectionINEL()->Reset();
   correctionStandard->GetTriggerBiasCorrectionNSD()->Reset();
   correctionStandard->GetTriggerBiasCorrectionND()->Reset();
+  correctionStandard->GetTriggerBiasCorrectionOnePart()->Reset();
 
   AlidNdEtaCorrection* corrections[100];
   TH1F* hRatios[100];
@@ -902,8 +1045,8 @@ void mergeCorrectionsWithDifferentCrosssections(Int_t correctionTarget = 3, Char
       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);
@@ -912,6 +1055,10 @@ void mergeCorrectionsWithDifferentCrosssections(Int_t correctionTarget = 3, Char
         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
@@ -927,6 +1074,14 @@ void mergeCorrectionsWithDifferentCrosssections(Int_t correctionTarget = 3, Char
 
       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;
 
@@ -936,7 +1091,7 @@ void mergeCorrectionsWithDifferentCrosssections(Int_t correctionTarget = 3, Char
       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_");
@@ -949,7 +1104,12 @@ void mergeCorrectionsWithDifferentCrosssections(Int_t correctionTarget = 3, Char
       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);
 
@@ -974,12 +1134,8 @@ void mergeCorrectionsWithDifferentCrosssections(Int_t correctionTarget = 3, Char
   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
@@ -988,51 +1144,120 @@ void CreateCorrectionsWithUA5CrossSections(Int_t origin, const Char_t* correctio
   //   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
@@ -1043,8 +1268,6 @@ void CreateCorrectionsWithUA5CrossSections(Int_t origin, const Char_t* correctio
 //       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");
@@ -1056,6 +1279,7 @@ void CreateCorrectionsWithUA5CrossSections(Int_t origin, const Char_t* correctio
   correctionStandard->GetTriggerBiasCorrectionINEL()->Reset();
   correctionStandard->GetTriggerBiasCorrectionNSD()->Reset();
   correctionStandard->GetTriggerBiasCorrectionND()->Reset();
+  correctionStandard->GetTriggerBiasCorrectionOnePart()->Reset();
 
   AlidNdEtaCorrection* corrections[100];
   TH1F* hRatios[100];
@@ -1114,6 +1338,10 @@ void CreateCorrectionsWithUA5CrossSections(Int_t origin, const Char_t* correctio
   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();
@@ -1128,6 +1356,14 @@ void CreateCorrectionsWithUA5CrossSections(Int_t origin, const Char_t* correctio
   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();
 
@@ -1784,3 +2020,173 @@ void ChangePtInCorrection(const char* fileName = "correction_map.root", const ch
   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()));
+}
index 4592693c8e31ec02d33636b125e113a16d891ece..0d9939550c41cb2bcb3b05b7a72d4985635d0c15 100644 (file)
@@ -68,7 +68,7 @@ TPad* DrawChange(Bool_t spd, const char* basename, const char** changes, Int_t n
 {  
   Float_t etaMax = 1.05;
   if (spd)
-    etaMax = 1.79;
+    etaMax = 1.49;
 
   TH1F* hRatios[100];
   for(Int_t i=0; i<nChanges; i++) {
@@ -79,7 +79,7 @@ TPad* DrawChange(Bool_t spd, const char* basename, const char** changes, Int_t n
     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);
@@ -173,7 +173,7 @@ void DrawEffectOfChangeInCrossSection(Bool_t spd = kFALSE, const char* fileName
   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};
@@ -200,14 +200,14 @@ void DrawEffectOfChangeInComposition(Bool_t spd = kFALSE, const char* fileName =
 {
   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) {
index 77479248b835a75f1ecb365db9e859926d4c2b6d..8c54df8081a00cd9389836a4dbab9bea48d0c29c 100644 (file)
@@ -1,7 +1,7 @@
 void Load(const char* taskName, Bool_t debug)
 {
   TString compileTaskName;
-  compileTaskName.Form("%s.cxx++", taskName);
+  compileTaskName.Form("%s.cxx+", taskName);
   if (debug)
     compileTaskName += "g";
 
@@ -57,10 +57,7 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
 
   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
@@ -85,8 +82,8 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
     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);
       
@@ -119,15 +116,30 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
   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))
@@ -154,6 +166,73 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
     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)
   {
@@ -169,15 +248,19 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
     //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);
 
@@ -188,6 +271,7 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
     cOutput = mgr->CreateContainer("cOutput", TList::Class(), AliAnalysisManager::kOutputContainer);
     mgr->ConnectOutput(task, 0, cOutput);
   }
+
   if (runWhat == 1 || runWhat == 2)
   {
     Load("AlidNdEtaCorrectionTask", aDebug);
@@ -198,10 +282,17 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
     //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);
 
@@ -213,7 +304,8 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
     mgr->ConnectOutput(task2, 0, cOutput);
   }
 
-  if (requiredData == 1) {
+  if (requiredData == 1) 
+  {
     // Enable MC event handler
     AliMCEventHandler* handler = new AliMCEventHandler;
     handler->SetReadTR(kFALSE);
@@ -247,9 +339,17 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
         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);
@@ -260,9 +360,15 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
       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)
       {
@@ -272,8 +378,12 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
       }
       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());
     }
@@ -282,8 +392,8 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
   {
     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
   {