introducing monalisa monitoring
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 Aug 2006 17:41:51 +0000 (17:41 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 Aug 2006 17:41:51 +0000 (17:41 +0000)
systematic uncertainty: particle composition

PWG0/dNdEta/AliMultiplicityESDSelector.cxx
PWG0/dNdEta/AliMultiplicityESDSelector.h
PWG0/dNdEta/AlidNdEtaAnalysisESDSelector.cxx
PWG0/dNdEta/AlidNdEtaCorrection.cxx
PWG0/dNdEta/AlidNdEtaCorrection.h
PWG0/dNdEta/AlidNdEtaSystematicsSelector.cxx
PWG0/dNdEta/AlidNdEtaSystematicsSelector.h
PWG0/dNdEta/drawSystematics.C
PWG0/dNdEta/makeSystematics.C

index 6046133f48ab47939d8478297d7c85294d1cebed..3b874ccdd08a23322145fea7f4c8f8c8a87d3c5b 100644 (file)
 #include "esdTrackCuts/AliESDtrackCuts.h"
 #include "AliPWG0Helper.h"
 
+// uncomment this to enable mona lisa monitoring
+//#define ALISELECTOR_USEMONALISA
+
+#ifdef ALISELECTOR_USEMONALISA
+  #include <TMonaLisaWriter.h>
+#endif
+
 ClassImp(AliMultiplicityESDSelector)
 
 AliMultiplicityESDSelector::AliMultiplicityESDSelector() :
   AliSelector(),
   fMultiplicity(0),
-  fEsdTrackCuts(0)
+  fEsdTrackCuts(0),
+  fMonaLisaWriter(0)
 {
   //
   // Constructor. Initialization of pointers
@@ -72,6 +80,24 @@ void AliMultiplicityESDSelector::SlaveBegin(TTree* tree)
   ReadUserObjects(tree);
 
   fMultiplicity = new TH1F("fMultiplicity", "multiplicity", 201, 0.5, 200.5);
+
+  #ifdef ALISELECTOR_USEMONALISA
+    TNamed *nm = 0;
+    if (fInput)
+      nm = dynamic_cast<TNamed*> (fInput->FindObject("PROOF_QueryTag"));
+    if (!nm)
+    {
+      AliDebug(AliLog::kError, "Query tag not found. Cannot enable monitoring");
+      return;
+    }
+
+    TString option = GetOption();
+    option.ReplaceAll("#+", "");
+
+    TString id;
+    id.Form("%s_%s%d", gSystem->HostName(), nm->GetTitle(), gSystem->GetPid());
+    fMonaLisaWriter = new TMonaLisaWriter(option, id, "CAF", "aliendb6.cern.ch");
+  #endif
 }
 
 Bool_t AliMultiplicityESDSelector::Process(Long64_t entry)
@@ -132,6 +158,14 @@ void AliMultiplicityESDSelector::SlaveTerminate()
 
   AliSelector::SlaveTerminate();
 
+  #ifdef ALISELECTOR_USEMONALISA
+    if (fMonaLisaWriter)
+    {
+      delete fMonaLisaWriter;
+      fMonaLisaWriter = 0;
+    }
+  #endif
+
   // Add the histograms to the output on each slave server
   if (!fOutput)
   {
index 3bfc48640c390ab02c9b16e1b5fd1dc69c998406..cb6f9d983873018f76930c8b6b6c5b278d0e7a99 100644 (file)
@@ -8,6 +8,8 @@
 class AliESDtrackCuts;
 class TH1F;
 
+class TMonaLisaWriter;
+
 class AliMultiplicityESDSelector : public AliSelector {
   public:
     AliMultiplicityESDSelector();
@@ -30,6 +32,8 @@ class AliMultiplicityESDSelector : public AliSelector {
     AliMultiplicityESDSelector(const AliMultiplicityESDSelector&);
     AliMultiplicityESDSelector& operator=(const AliMultiplicityESDSelector&);
 
+    TMonaLisaWriter* fMonaLisaWriter; //! ML instance for monitoring
+
   ClassDef(AliMultiplicityESDSelector, 0);
 };
 
index 1e01297b9b57f604a1a17fb0da7b27b02da8ede0..017853bdd3dca1eeeb9ad57ec9a9f5a00d8ede4c 100644 (file)
@@ -27,7 +27,7 @@ AlidNdEtaAnalysisESDSelector::AlidNdEtaAnalysisESDSelector() :
   fdNdEtaAnalysis(0),
   fEsdTrackCuts(0),
   fdNdEtaCorrection(0)
-{
+  {
   //
   // Constructor. Initialization of pointers
   //
@@ -62,7 +62,7 @@ AlidNdEtaAnalysisESDSelector::~AlidNdEtaAnalysisESDSelector()
     fdNdEtaAnalysisMBVtx = 0;
   }
 
-  if (fEsdTrackCuts)
+  /*if (fEsdTrackCuts)
   {
     delete fEsdTrackCuts;
     fEsdTrackCuts = 0;
@@ -72,7 +72,7 @@ AlidNdEtaAnalysisESDSelector::~AlidNdEtaAnalysisESDSelector()
   {
     delete fdNdEtaCorrection;
     fdNdEtaCorrection = 0;
-  }
+  }*/
 }
 
 void AlidNdEtaAnalysisESDSelector::Begin(TTree* tree)
index add6030ce54d0018c379caab849e23ce828b93a9..b3476c7d630eb83797582187a04d319f647c97ab 100644 (file)
@@ -15,9 +15,9 @@ AlidNdEtaCorrection::AlidNdEtaCorrection()
   fTrack2ParticleCorrection(0),
   fVertexRecoCorrection(0),
   fTriggerCorrection(0),
-  fTriggerBiasCorrection(0),
-  fNEvents(0),
-  fNTriggeredEvents(0)
+  fTriggerBiasCorrection(0)
+  //fNEvents(0),
+  //fNTriggeredEvents(0)
 {
   // default constructor
 }
@@ -28,9 +28,9 @@ AlidNdEtaCorrection::AlidNdEtaCorrection(const Char_t* name, const Char_t* title
   fTrack2ParticleCorrection(0),
   fVertexRecoCorrection(0),
   fTriggerCorrection(0),
-  fTriggerBiasCorrection(0),
-  fNEvents(0),
-  fNTriggeredEvents(0)
+  fTriggerBiasCorrection(0)
+  //fNEvents(0),
+  //fNTriggeredEvents(0)
 {
   // constructor
   //
@@ -118,25 +118,25 @@ AlidNdEtaCorrection::Finish() {
   fVertexRecoCorrection->Divide();
   fTriggerCorrection->Divide();
 
-  if (fNEvents == 0)
+  if (GetNevents() <= 0)
   {
-    printf("ERROR: fNEvents is empty. Cannot scale histogram. Skipping processing of fTriggerBiasCorrection\n");
+    printf("ERROR: GetNevents() returns 0. Cannot scale histogram. Skipping processing of fTriggerBiasCorrection\n");
     return;
   }
-  fTriggerBiasCorrection->GetMeasuredHistogram()->Scale(Double_t(fNTriggeredEvents)/Double_t(fNEvents));
+  fTriggerBiasCorrection->GetMeasuredHistogram()->Scale(Double_t(GetNtriggeredEvents())/Double_t(GetNevents()));
   fTriggerBiasCorrection->Divide();
 }
 
 //____________________________________________________________________
-Long64_t 
+Long64_t
 AlidNdEtaCorrection::Merge(TCollection* list) {
   // Merge a list of dNdEtaCorrection objects with this (needed for
-  // PROOF). 
+  // PROOF).
   // Returns the number of merged objects (including this).
 
   if (!list)
     return 0;
-  
+
   if (list->IsEmpty())
     return 1;
 
@@ -147,36 +147,38 @@ AlidNdEtaCorrection::Merge(TCollection* list) {
   TList* collectionNtrackToNparticle = new TList;
   TList* collectionVertexReco        = new TList;
   TList* collectionTriggerBias       = new TList;
+  TList* collectionTrigger           = new TList;
 
   Int_t count = 0;
   while ((obj = iter->Next())) {
-    
+
     AlidNdEtaCorrection* entry = dynamic_cast<AlidNdEtaCorrection*> (obj);
-    if (entry == 0) 
+    if (entry == 0)
       continue;
 
     collectionNtrackToNparticle ->Add(entry->GetTrack2ParticleCorrection());
     collectionVertexReco        ->Add(entry->GetVertexRecoCorrection());
     collectionTriggerBias        ->Add(entry->GetTriggerBiasCorrection());
+    collectionTrigger        ->Add(entry->GetTriggerCorrection());
 
     count++;
 
-    fNEvents += entry->fNEvents;
-    fNTriggeredEvents += entry->fNTriggeredEvents;
+    //fNEvents += entry->fNEvents;
+    //fNTriggeredEvents += entry->fNTriggeredEvents;
   }
   fTrack2ParticleCorrection ->Merge(collectionNtrackToNparticle);
   fVertexRecoCorrection        ->Merge(collectionVertexReco);
   fTriggerBiasCorrection        ->Merge(collectionTriggerBias);
+  fTriggerCorrection        ->Merge(collectionTrigger);
 
   delete collectionNtrackToNparticle;
   delete collectionVertexReco;
   delete collectionTriggerBias;
-
+  delete collectionTrigger;
 
   return count+1;
 }
 
-
 //____________________________________________________________________
 Bool_t
 AlidNdEtaCorrection::LoadHistograms(const Char_t* fileName, const Char_t* dir) {
@@ -192,7 +194,6 @@ AlidNdEtaCorrection::LoadHistograms(const Char_t* fileName, const Char_t* dir) {
   return kTRUE;
 }
 
-
 //____________________________________________________________________
 void
 AlidNdEtaCorrection::SaveHistograms() {
@@ -282,3 +283,23 @@ void AlidNdEtaCorrection::ReduceInformation()
   fTriggerCorrection->ReduceInformation();
   fTriggerBiasCorrection->ReduceInformation();
 }
+
+Int_t AlidNdEtaCorrection::GetNevents()
+{
+  // Return total number of events used for this correction map (taken from fTriggerCorrection)
+
+  if (!fTriggerCorrection)
+    return -1;
+
+  return (Int_t) fTriggerCorrection->GetGeneratedHistogram()->Integral(0, fTriggerCorrection->GetGeneratedHistogram()->GetNbinsX()+1);
+}
+
+Int_t AlidNdEtaCorrection::GetNtriggeredEvents()
+{
+  // Return total number of triggered events used for this correction map (taken from fTriggerCorrection)
+
+  if (!fTriggerCorrection)
+    return -1;
+
+  return (Int_t) fTriggerCorrection->GetMeasuredHistogram()->Integral(0, fTriggerCorrection->GetMeasuredHistogram()->GetNbinsX()+1);
+}
index 7b249df26a9b935e6349e5a7a17dcd5042d95d1c..eb1703106e50ccd57fc4ebce771bf3eb2216ed51 100644 (file)
@@ -44,14 +44,15 @@ public:
   void FillParticleAllEvents(Float_t eta, Float_t pt)          {fTriggerBiasCorrection->FillGene(eta, pt);}
   void FillParticleWhenEventTriggered(Float_t eta, Float_t pt) {fTriggerBiasCorrection->FillMeas(eta, pt);}
 
-  void IncreaseEventCount() { fNEvents++; }
-  void IncreaseTriggeredEventCount() { fNTriggeredEvents++; }
+  //void IncreaseEventCount() { fNEvents++; }
+  //void IncreaseTriggeredEventCount() { fNTriggeredEvents++; }
 
   void Finish();
 
   AliCorrectionMatrix3D* GetTrack2ParticleCorrection()    {return fTrack2ParticleCorrection;}
   AliCorrectionMatrix2D* GetVertexRecoCorrection()        {return fVertexRecoCorrection;}
   AliCorrectionMatrix2D* GetTriggerBiasCorrection()       {return fTriggerBiasCorrection;}
+  AliCorrectionMatrix2D* GetTriggerCorrection()       {return fTriggerCorrection;}
 
   virtual Long64_t Merge(TCollection* list);
 
@@ -64,7 +65,7 @@ public:
   
   //  void RemoveEdges(Float_t cut=2, Int_t nBinsVtx=0, Int_t nBinsEta=0);
   
-  Float_t GetTrack2ParticleCorrection(Float_t vtx, Float_t eta, Float_t pt) 
+  Float_t GetTrack2ParticleCorrection(Float_t vtx, Float_t eta, Float_t pt)
     {return fTrack2ParticleCorrection->GetCorrection(vtx, eta, pt);}
 
   Float_t GetVertexRecoCorrection(Float_t vtx, Float_t n) {return fVertexRecoCorrection->GetCorrection(vtx, n);}
@@ -75,7 +76,10 @@ public:
 
   Float_t GetMeasuredFraction(Float_t ptCutOff, Float_t eta = -1, Bool_t debug = kFALSE);
 
-  void SetNEvents(Long64_t events) { fNEvents = events; }
+  //void SetNEvents(Long64_t events) { fNEvents = events; }
+
+  Int_t GetNevents();
+  Int_t GetNtriggeredEvents();
 
   void ReduceInformation();
 
@@ -86,8 +90,8 @@ protected:
 
   AliCorrectionMatrix2D* fTriggerBiasCorrection;          //-> MB to desired sample
 
-  Long64_t fNEvents;
-  Long64_t fNTriggeredEvents;
+  //Long64_t fNEvents;
+  //Long64_t fNTriggeredEvents;
 
 private:
   AlidNdEtaCorrection(const AlidNdEtaCorrection&);
index 7750b22d1669ba1419ee59e3d4b47a2d9a6adaa7..7807fa8d4ad33bb75e14e90e73ffcbdebd150a87 100644 (file)
@@ -8,7 +8,7 @@
 #include <TVector3.h>
 #include <TChain.h>
 #include <TFile.h>
-#include <TH3F.h>
+#include <TH2F.h>
 #include <TH1F.h>
 #include <TParticle.h>
 
@@ -29,8 +29,6 @@ AlidNdEtaSystematicsSelector::AlidNdEtaSystematicsSelector() :
   fSecondaries(0),
   fSigmaVertex(0),
   fEsdTrackCuts(0),
-  fOverallPrimaries(0),
-  fOverallSecondaries(0),
   fPIDParticles(0),
   fPIDTracks(0)
 {
@@ -88,7 +86,15 @@ void AlidNdEtaSystematicsSelector::SlaveBegin(TTree* tree)
 
   if (option.Contains("secondaries"))
   {
-    fSecondaries = new TH3F("fSecondaries", "fSecondaries;NacceptedTracks;CratioSecondaries;p_{T}", 2000, -0.5, 205.5, 16, 0.45, 2.05, 10, 0, 10);
+    fSecondaries = new TH2F("fSecondaries", "fSecondaries;Case;N", 8, 0.5, 8.5, 2001, -0.5, 2000.5);
+    fSecondaries->GetXaxis()->SetBinLabel(1, "All Primaries");
+    fSecondaries->GetXaxis()->SetBinLabel(2, "All Secondaries");
+    fSecondaries->GetXaxis()->SetBinLabel(3, "Primaries pT > 0.3 GeV/c");
+    fSecondaries->GetXaxis()->SetBinLabel(4, "Secondaries pT > 0.3 GeV/c");
+    fSecondaries->GetXaxis()->SetBinLabel(5, "Tracks from Primaries");
+    fSecondaries->GetXaxis()->SetBinLabel(6, "Tracks from Secondaries");
+    fSecondaries->GetXaxis()->SetBinLabel(7, "Accepted Tracks from Primaries");
+    fSecondaries->GetXaxis()->SetBinLabel(8, "Accepted Tracks from Secondaries");
   }
 
   if (option.Contains("particle-composition"))
@@ -160,7 +166,7 @@ Bool_t AlidNdEtaSystematicsSelector::Process(Long64_t entry)
     FillCorrectionMaps(list);
 
   if (fSecondaries)
-    FillSecondaries(list);
+    FillSecondaries();
 
   if (fSigmaVertex)
     FillSigmaVertex();
@@ -219,17 +225,17 @@ void AlidNdEtaSystematicsSelector::FillCorrectionMaps(TObjArray* listOfTracks)
       case 211: id = 0; break;
       case 321: id = 1; break;
       case 2212: id = 2; break;
-      case 11: 
+      case 11:
         {
-          if (pt < 0.1 && particle->GetMother(0) > -1)
+          /*if (pt < 0.1 && particle->GetMother(0) > -1)
           {
             TParticle* mother = stack->Particle(particle->GetMother(0));
             printf("Mother pdg code is %d\n", mother->GetPdgCode());
-          }   
+          } */
           //particle->Dump();
           //if (particle->GetUniqueID() == 1)
 
-          //printf("Found illegal particle mc id = %d file = %s event = %d\n", iMc,  fTree->GetCurrentFile()->GetName(), fTree->GetTree()->GetReadEntry()); 
+          //printf("Found illegal particle mc id = %d file = %s event = %d\n", iMc,  fTree->GetCurrentFile()->GetName(), fTree->GetTree()->GetReadEntry());
         }
       default: id = 3; break;
     }
@@ -237,7 +243,7 @@ void AlidNdEtaSystematicsSelector::FillCorrectionMaps(TObjArray* listOfTracks)
     if (vertexReconstructed)
     {
       fdNdEtaCorrection[id]->FillParticle(vtxMC[2], eta, pt);
-      if (pt < 0.1)
+      //if (pt < 0.1)
         fPIDParticles->Fill(particle->GetPdgCode());
     }
 
@@ -264,8 +270,36 @@ void AlidNdEtaSystematicsSelector::FillCorrectionMaps(TObjArray* listOfTracks)
       continue;
     }
 
+    TParticle* mother = particle;
+    // find primary particle that created this particle
+    while (label >= nPrim)
+    {
+      //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
+
+      if (mother->GetMother(0) == -1)
+      {
+        AliDebug(AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
+        mother = 0;
+        break;
+      }
+
+      label = mother->GetMother(0);
+
+      mother = stack->Particle(label);
+      if (!mother)
+      {
+        AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (find mother loop).", label));
+        break;
+      }
+    }
+
+    if (!mother)
+      continue;
+
+    //printf("We continue with particle %d (pdg %d)\n", label, mother->GetPdgCode());
+
     Int_t id = -1;
-    switch (TMath::Abs(particle->GetPdgCode()))
+    switch (TMath::Abs(mother->GetPdgCode()))
     {
       case 211:  id = 0; break;
       case 321:  id = 1; break;
@@ -276,31 +310,92 @@ void AlidNdEtaSystematicsSelector::FillCorrectionMaps(TObjArray* listOfTracks)
     if (vertexReconstructed)
     {
       fdNdEtaCorrection[id]->FillParticleWhenMeasuredTrack(vtxMC[2], particle->Eta(), particle->Pt());
-      if (particle->Pt() < 0.1)
+      //if (particle->Pt() < 0.1)
         fPIDTracks->Fill(particle->GetPdgCode());
     }
   } // end of track loop
 
+  Int_t nGoodTracks = listOfTracks->GetEntries();
+
+  for (Int_t i=0; i<4; ++i)
+  {
+    fdNdEtaCorrection[i]->FillEvent(vtxMC[2], nGoodTracks);
+    if (eventTriggered)
+    {
+      fdNdEtaCorrection[i]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
+      if (vertexReconstructed)
+        fdNdEtaCorrection[i]->FillEventWithTriggerWithReconstructedVertex(vtxMC[2], nGoodTracks);
+    }
+  }
+
   delete iter;
   iter = 0;
 }
 
-void AlidNdEtaSystematicsSelector::FillSecondaries(TObjArray* listOfTracks)
+void AlidNdEtaSystematicsSelector::FillSecondaries()
 {
   // fills the secondary histograms
 
   AliStack* stack = GetStack();
 
-  TH1* nPrimaries = new TH1F("secondaries_primaries", "secondaries_primaries", fSecondaries->GetZaxis()->GetNbins(), fSecondaries->GetZaxis()->GetXmin(), fSecondaries->GetZaxis()->GetXmax());
-  TH1* nSecondaries = new TH1F("secondaries_secondaries", "secondaries_secondaries", fSecondaries->GetZaxis()->GetNbins(), fSecondaries->GetZaxis()->GetXmin(), fSecondaries->GetZaxis()->GetXmax());
+  Int_t particlesPrimaries = 0;
+  Int_t particlesSecondaries = 0;
+  Int_t particlesPrimariesPtCut = 0;
+  Int_t particlesSecondariesPtCut = 0;
 
-  TIterator* iter = listOfTracks->MakeIterator();
-  TObject* obj = 0;
-  while ((obj = iter->Next()) != 0)
+  for (Int_t iParticle = 0; iParticle < stack->GetNtrack(); iParticle++)
   {
-    AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
+    TParticle* particle = stack->Particle(iParticle);
+    if (!particle)
+    {
+      AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (particle loop).", iParticle));
+      continue;
+    }
+
+    if (TMath::Abs(particle->Eta()) > 0.9)
+      continue;
+
+    Bool_t isPrimary = kFALSE;
+    if (iParticle < stack->GetNprimary())
+    {
+      if (AliPWG0Helper::IsPrimaryCharged(particle, stack->GetNprimary()) == kFALSE)
+        continue;
+
+      isPrimary = kTRUE;
+    }
+
+    if (isPrimary)
+      particlesPrimaries++;
+    else
+      particlesSecondaries++;
+
+    if (particle->Pt() < 0.3)
+      continue;
+
+    if (isPrimary)
+      particlesPrimariesPtCut++;
+    else
+      particlesSecondariesPtCut++;
+  }
+
+  fSecondaries->Fill(1, particlesPrimaries);
+  fSecondaries->Fill(2, particlesSecondaries);
+  fSecondaries->Fill(3, particlesPrimariesPtCut);
+  fSecondaries->Fill(4, particlesSecondariesPtCut);
+
+  Int_t allTracksPrimaries = 0;
+  Int_t allTracksSecondaries = 0;
+  Int_t acceptedTracksPrimaries = 0;
+  Int_t acceptedTracksSecondaries = 0;
+
+  for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++)
+  {
+    AliESDtrack* esdTrack = fESD->GetTrack(iTrack);
     if (!esdTrack)
+    {
+      AliDebug(AliLog::kError, Form("UNEXPECTED: Could not get track %d.", iTrack));
       continue;
+    }
 
     Int_t label = TMath::Abs(esdTrack->GetLabel());
     TParticle* particle = stack->Particle(label);
@@ -310,42 +405,26 @@ void AlidNdEtaSystematicsSelector::FillSecondaries(TObjArray* listOfTracks)
       continue;
     }
 
-    Int_t nPrim  = stack->GetNprimary();
-    if (label < nPrim)
-      nPrimaries->Fill(particle->Pt());
+   if (label < stack->GetNprimary())
+      allTracksPrimaries++;
     else
-      nSecondaries->Fill(particle->Pt());
-  }
-
-  for (Int_t i=1; i<=nPrimaries->GetNbinsX(); ++i)
-  {
-    Int_t primaries = (Int_t) nPrimaries->GetBinContent(i);
-    Int_t secondaries = (Int_t) nSecondaries->GetBinContent(i);
+      allTracksSecondaries++;
 
-    if (primaries + secondaries > 0)
-    {
-      AliDebug(AliLog::kDebug, Form("The ratio between primaries and secondaries is %d:%d = %f", primaries, secondaries, ((secondaries > 0) ? (Double_t) primaries / secondaries : -1)));
+    if (!fEsdTrackCuts->AcceptTrack(esdTrack))
+      continue;
 
-      for (Double_t factor = 0.5; factor < 2.01; factor += 0.1)
-      {
-        Double_t nTracks = (Double_t) primaries + (Double_t) secondaries * factor;
-        fSecondaries->Fill(nTracks, factor, nPrimaries->GetBinCenter(i));
-        //if (secondaries > 0) printf("We fill: %f %f %f\n", nTracks, factor, nPrimaries->GetBinCenter(i));
-      }
-    }
+    if (label < stack->GetNprimary())
+      acceptedTracksPrimaries++;
+    else
+      acceptedTracksSecondaries++;
   }
 
-  fOverallPrimaries += (Int_t) nPrimaries->Integral();
-  fOverallSecondaries += (Int_t) nSecondaries->Integral();
-
-  delete nPrimaries;
-  nPrimaries = 0;
-
-  delete nSecondaries;
-  nSecondaries = 0;
+  fSecondaries->Fill(5, allTracksPrimaries);
+  fSecondaries->Fill(6, allTracksSecondaries);
+  fSecondaries->Fill(7, acceptedTracksPrimaries);
+  fSecondaries->Fill(8, acceptedTracksSecondaries);
 
-  delete iter;
-  iter = 0;
+  //printf("P = %d, S = %d, P_pt = %d, S_pt = %d, T_P = %d, T_S = %d, T_P_acc = %d, T_S_acc = %d\n", particlesPrimaries, particlesSecondaries, particlesPrimariesPtCut, particlesSecondariesPtCut, allTracksPrimaries, allTracksSecondaries, acceptedTracksPrimaries, acceptedTracksSecondaries);
 }
 
 void AlidNdEtaSystematicsSelector::FillSigmaVertex()
@@ -421,7 +500,7 @@ void AlidNdEtaSystematicsSelector::Terminate()
 
   AliSelector::Terminate();
 
-  fSecondaries = dynamic_cast<TH3F*> (fOutput->FindObject("fSecondaries"));
+  fSecondaries = dynamic_cast<TH2F*> (fOutput->FindObject("fSecondaries"));
   for (Int_t i=0; i<4; ++i)
     fdNdEtaCorrection[i] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
   fSigmaVertex = dynamic_cast<TH1F*> (fOutput->FindObject("fSigmaVertex"));
@@ -441,10 +520,7 @@ void AlidNdEtaSystematicsSelector::Terminate()
     fEsdTrackCuts->SaveHistograms("esd_track_cuts");
 
   if (fSecondaries)
-  {
     fSecondaries->Write();
-    printf("We had %d primaries and %d secondaries.\n", (Int_t) fOverallPrimaries, (Int_t) fOverallSecondaries);
-  }
 
   if (fSigmaVertex)
     fSigmaVertex->Write();
@@ -459,7 +535,7 @@ void AlidNdEtaSystematicsSelector::Terminate()
   if (fSecondaries)
   {
     new TCanvas;
-    fSecondaries->Draw();
+    fSecondaries->Draw("COLZ");
   }
 
   if (fSigmaVertex)
index 727e9dcc5fa2ca63cf580fe57d47c6aea572b23c..e805381e85c85421d5b7d6daefafb9286a2d136d 100644 (file)
@@ -7,7 +7,8 @@
 
 class AliESDtrackCuts;
 class AlidNdEtaCorrection;
-class TH3F;
+
+class TH2F;
 class TH1F;
 
 class AlidNdEtaSystematicsSelector : public AliSelectorRL {
@@ -25,18 +26,17 @@ class AlidNdEtaSystematicsSelector : public AliSelectorRL {
     void ReadUserObjects(TTree* tree);
 
     void FillCorrectionMaps(TObjArray* listOfTracks);
-    void FillSecondaries(TObjArray* listOfTracks);
+    void FillSecondaries();
     void FillSigmaVertex();
 
-    TH3F* fSecondaries; // (accepted tracks) vs (tracks from sec)/(n * tracks from sec) vs pT
+    TH2F* fSecondaries; // (Nprim/Nsec for the cases: all/above3GeV/reconstructed tracks/accepted tracks) vs (particle count)
+
     AlidNdEtaCorrection* fdNdEtaCorrection[4];      // correction for different particle species: here pi, K, p, others
+
     TH1F* fSigmaVertex; // (accepted tracks) vs (n of sigma to vertex cut)
 
     AliESDtrackCuts* fEsdTrackCuts;     // Object containing the parameters of the esd track cuts
 
-    Long64_t fOverallPrimaries; // count of all primaries
-    Long64_t fOverallSecondaries; // count of all secondaries
-
     TH1F* fPIDParticles; // pid of primary particles
     TH1F* fPIDTracks; // pid of reconstructed tracks
 
index f24cf5f624ace53ca13a1c3e8bd3f9464fac59a2..6f834ce8ac0618baf3caa72e6971313399df39a7 100644 (file)
@@ -60,12 +60,16 @@ void Prepare2DPlot(TH2* hist)
   SetRanges(hist);
 }
 
-void Prepare1DPlot(TH1* hist)
+void Prepare1DPlot(TH1* hist, Bool_t setRanges = kTRUE)
 {
   hist->SetLineWidth(2);
   hist->SetStats(kFALSE);
 
-  SetRanges(hist);
+  hist->GetXaxis()->SetTitleOffset(1.2);
+  hist->GetYaxis()->SetTitleOffset(1.2);
+
+  if (setRanges)
+    SetRanges(hist);
 }
 
 void InitPad()
@@ -97,51 +101,22 @@ void Secondaries()
 {
   TFile* file = TFile::Open("systematics.root");
 
-  TH3F* secondaries = dynamic_cast<TH3F*> (file->Get("fSecondaries"));
+  TH2F* secondaries = dynamic_cast<TH2F*> (file->Get("fSecondaries"));
   if (!secondaries)
   {
     printf("Could not read histogram\n");
     return;
   }
 
-  for (Int_t ptBin=1; ptBin<=secondaries->GetNbinsZ(); ptBin++)
-  //for (Int_t ptBin = 1; ptBin<=2; ptBin++)
+  TCanvas* canvas = new TCanvas("Secondaries", "Secondaries", 1000, 1000);
+  canvas->Divide(3, 3);
+  for (Int_t i=1; i<=8; i++)
   {
-    TH1F* hist = new TH1F(Form("secondaries_%d", ptBin), Form("secondaries_%d", ptBin), secondaries->GetNbinsY(), secondaries->GetYaxis()->GetXmin(), secondaries->GetYaxis()->GetXmax());
-
-    hist->SetTitle(Form("%f < p_{T} < %f", secondaries->GetZaxis()->GetBinLowEdge(ptBin), secondaries->GetZaxis()->GetBinUpEdge(ptBin)));
+    TH1D* hist = secondaries->ProjectionY(Form("proj_%d", i), i, i);
+    hist->SetTitle(secondaries->GetXaxis()->GetBinLabel(i));
 
-    for (Int_t cBin=1; cBin<=secondaries->GetNbinsY(); ++cBin)
-    {
-      if (secondaries->GetBinContent(0, cBin, ptBin) > 0)
-        printf("WARNING: Underflow bin not empty!");
-      if (secondaries->GetBinContent(secondaries->GetNbinsX()+1, cBin, ptBin) > 0)
-        printf("WARNING: Overflow bin not empty!");
-
-      Double_t sum = 0;
-      Double_t count = 0;
-      for (Int_t nBin=1; nBin<=secondaries->GetNbinsX(); ++nBin)
-      {
-        //printf("%f %f\n", secondaries->GetXaxis()->GetBinCenter(nBin), secondaries->GetBinContent(nBin, cBin, ptBin));
-        sum += secondaries->GetXaxis()->GetBinCenter(nBin) * secondaries->GetBinContent(nBin, cBin, ptBin);
-        count += secondaries->GetBinContent(nBin, cBin, ptBin);
-      }
-
-      printf("%f %f\n", sum, count);
-
-      if (count > 0)
-      {
-        hist->SetBinContent(cBin, sum / count);
-        hist->SetBinError(cBin, TMath::Sqrt(sum) / count);
-      }
-    }
-
-    hist->Scale(1.0 / hist->GetBinContent(hist->GetXaxis()->FindBin(1)));
-    hist->Add(new TF1("flat", "-1", 0, 2));
-
-    new TCanvas;
-    hist->SetMarkerStyle(21);
-    hist->Draw("");
+    canvas->cd(i);
+    hist->Draw();
   }
 }
 
@@ -214,7 +189,7 @@ TH1** DrawRatios(const char* fileName = "systematics.root")
 
   TH1** ptDists = new TH1*[4];
 
-  TLegend* legend = new TLegend(0.7, 0.7, 0.95, 0.95);
+  TLegend* legend = new TLegend(0.73, 0.73, 0.98, 0.98);
 
   const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3" };
   const char* particleNames[] = { "#pi", "K", "p", "other" };
@@ -255,19 +230,43 @@ TH1** DrawRatios(const char* fileName = "systematics.root")
   for (Int_t i=0; i<4; ++i)
   {
     ptDists[i]->Divide(total);
-    ptDists[i]->SetTitle(";p_{T};Ratio");
+    ptDists[i]->SetStats(kFALSE);
+    ptDists[i]->SetTitle(";p_{T};Fraction of total");
     ptDists[i]->GetYaxis()->SetRangeUser(0, 1);
     ptDists[i]->Draw((i == 0) ? "" : "SAME");
   }
+  legend->SetFillColor(0);
   legend->Draw();
 
   canvas->SaveAs("DrawRatios.gif");
 
+
+  canvas = new TCanvas("PythiaRatios", "PythiaRatios", 500, 500);
+  for (Int_t i=0; i<4; ++i)
+  {
+    TH1* hist = ptDists[i]->Clone();
+    hist->GetXaxis()->SetRangeUser(0, 1.9);
+    hist->Draw((i == 0) ? "" : "SAME");
+  }
+  legend->Draw();
+
+  canvas->SaveAs("pythiaratios.eps");
+
   file->Close();
 
   return ptDists;
 }
 
+void DrawCompareToReal()
+{
+  gROOT->ProcessLine(".L drawPlots.C");
+
+  const char* fileNames[] = { "correction_map.root", "new_compositions.root" };
+  const char* folderNames[] = { "dndeta_correction", "PythiaRatios" };
+
+  Track2Particle1DComposition(fileNames, 2, folderNames);
+}
+
 void DrawDifferentSpecies()
 {
   gROOT->ProcessLine(".L drawPlots.C");
@@ -308,7 +307,7 @@ void DrawBoosts()
   Track2Particle1DComposition(fileNames, 4, folderNames);
 }
 
-TH1F* HistogramDifferences(const char* filename, const char* folder1, const char* folder2)
+TH2F* HistogramDifferences(const char* filename, const char* folder1, const char* folder2)
 {
   gSystem->Load("libPWG0base");
 
@@ -322,15 +321,19 @@ TH1F* HistogramDifferences(const char* filename, const char* folder1, const char
   fdNdEtaCorrection[1] = new AlidNdEtaCorrection(folder2, folder2);
   fdNdEtaCorrection[1]->LoadHistograms(filename, folder2);
 
-  TH1F* difference = new TH1F("difference", Form(";#DeltaC_{pT, z, #eta} %s - %s;Count", folder2, folder1), 1000, -0.5, 0.5);
-
   TH3F* hist1 = fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
   TH3F* hist2 = fdNdEtaCorrection[1]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
 
+  //TH1F* difference = new TH1F("difference", Form(";#DeltaC_{pT, z, #eta} %s / %s;Count", folder2, folder1), 1000, 0.9, 1.1);
+  TH2F* difference = new TH2F(Form("difference_%s_%s", folder1, folder2), Form(";#Sigma (C_{pT, z} %s / C_{pT, z} %s);#eta;Count", folder2, folder1), 100, 0.9, 1.1, hist1->GetYaxis()->GetNbins(), hist1->GetYaxis()->GetXmin(), hist1->GetYaxis()->GetXmax());
+
   for (Int_t x=hist1->GetXaxis()->FindBin(-10); x<=hist1->GetXaxis()->FindBin(10); ++x)
     for (Int_t y=hist1->GetYaxis()->FindBin(-0.8); y<=hist1->GetYaxis()->FindBin(0.8); ++y)
       for (Int_t z=hist1->GetZaxis()->FindBin(0.3); z<=hist1->GetZaxis()->FindBin(9.9); ++z)
-        difference->Fill(hist2->GetBinContent(x, y, z) - hist1->GetBinContent(x, y, z));
+        if (hist1->GetBinContent(x, y, z) != 0)
+          difference->Fill(hist2->GetBinContent(x, y, z) / hist1->GetBinContent(x, y, z), hist1->GetYaxis()->GetBinCenter(y));
+
+  difference->GetYaxis()->SetRangeUser(-0.8, 0.8);
 
   printf("Over-/Underflow bins: %d %d\n", difference->GetBinContent(0), difference->GetBinContent(difference->GetNbinsX()+1));
 
@@ -339,26 +342,78 @@ TH1F* HistogramDifferences(const char* filename, const char* folder1, const char
 
 void HistogramDifferences()
 {
-  TH1F* PiBoosted = HistogramDifferences("new_compositions.root", "PythiaRatios", "PiBoosted");
-  TH1F* KBoosted = HistogramDifferences("new_compositions.root", "PythiaRatios", "KBoosted");
-  TH1F* pBoosted = HistogramDifferences("new_compositions.root", "PythiaRatios", "pBoosted");
+  TH2F* KBoosted = HistogramDifferences("new_compositions.root", "PythiaRatios", "KBoosted");
+  TH2F* pBoosted = HistogramDifferences("new_compositions.root", "PythiaRatios", "pBoosted");
+  TH2F* KReduced = HistogramDifferences("new_compositions.root", "PythiaRatios", "KReduced");
+  TH2F* pReduced = HistogramDifferences("new_compositions.root", "PythiaRatios", "pReduced");
 
-  TCanvas* canvas = new TCanvas("HistogramDifferences", "HistogramDifferences", 1200, 400);
-  canvas->Divide(3, 1);
+  TCanvas* canvas = new TCanvas("HistogramDifferences", "HistogramDifferences", 1000, 1000);
+  canvas->Divide(2, 2);
 
   canvas->cd(1);
-  PiBoosted->GetXaxis()->SetRangeUser(-0.01, 0.01);
-  PiBoosted->Draw();
+  KBoosted->GetXaxis()->SetRangeUser(-0.05, 0.05);
+  KBoosted->Draw("COLZ");
 
   canvas->cd(2);
-  KBoosted->GetXaxis()->SetRangeUser(-0.01, 0.01);
-  KBoosted->Draw();
+  KReduced->GetXaxis()->SetRangeUser(-0.05, 0.05);
+  KReduced->Draw("COLZ");
 
   canvas->cd(3);
-  pBoosted->GetXaxis()->SetRangeUser(-0.01, 0.01);
-  pBoosted->Draw();
+  pBoosted->GetXaxis()->SetRangeUser(-0.02, 0.02);
+  pBoosted->Draw("COLZ");
+
+  canvas->cd(4);
+  pReduced->GetXaxis()->SetRangeUser(-0.02, 0.02);
+  pReduced->Draw("COLZ");
 
   canvas->SaveAs("HistogramDifferences.gif");
+
+  canvas = new TCanvas("HistogramDifferencesProfile", "HistogramDifferencesProfile", 1000, 500);
+  canvas->Divide(2, 1);
+
+  canvas->cd(1);
+  gPad->SetBottomMargin(0.13);
+  KBoosted->SetTitle("a) Ratio of correction maps");
+  KBoosted->SetStats(kFALSE);
+  KBoosted->GetXaxis()->SetTitleOffset(1.4);
+  KBoosted->GetXaxis()->SetLabelOffset(0.02);
+  KBoosted->Draw("COLZ");
+
+  canvas->cd(2);
+  gPad->SetGridx();
+  gPad->SetGridy();
+
+  TLegend* legend = new TLegend(0.73, 0.73, 0.98, 0.98);
+
+  for (Int_t i=0; i<4; ++i)
+  {
+    TH2F* hist = 0;
+    TString name;
+    switch (i)
+    {
+      case 0: hist = KBoosted; name = "K enhanced"; break;
+      case 1: hist = KReduced; name = "K reduced"; break;
+      case 2: hist = pBoosted; name = "p enhanced"; break;
+      case 3: hist = pReduced; name = "p reduced"; break;
+    }
+
+    TProfile* profile = hist->ProfileY();
+    profile->SetTitle("b) Mean and RMS");
+    profile->GetXaxis()->SetRange(hist->GetYaxis()->GetFirst(), hist->GetYaxis()->GetLast());
+    profile->GetXaxis()->SetTitleOffset(1.2);
+    profile->GetXaxis()->SetLabelOffset(0.02);
+    profile->GetYaxis()->SetRangeUser(0.98, 1.02);
+    profile->SetStats(kFALSE);
+    profile->SetLineColor(i+1);
+    profile->SetMarkerColor(i+1);
+    profile->DrawCopy(((i > 0) ? "SAME" : ""));
+
+
+    legend->AddEntry(profile, name);
+  }
+
+  legend->Draw();
+  canvas->SaveAs("particlecomposition_result.eps");
 }
 
 
@@ -416,26 +471,44 @@ const char* ChangeComposition(void** correctionPointer, Int_t index)
       break;
 
     case 1: // each species rated with pythia ratios
-      TH1** ptDists = DrawRatios("pythiaratios.root");
+      /*TH1** ptDists = DrawRatios("pythiaratios.root");
 
       for (Int_t i=0; i<3; ++i)
       {
         ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram(), ptDists[i]);
         ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram(), ptDists[i]);
-      }
+      }*/
       return "PythiaRatios";
       break;
 
+      // one species enhanced / reduced
+    case 2: // + 50% pions
+    case 3: // - 50% pions
+    case 4: // + 50% kaons
+    case 5: // - 50% kaons
+    case 6: // + 50% protons
+    case 7: // - 50% protons
+      Int_t correctionIndex = (index - 2) / 2;
+      Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5;
+
+      fdNdEtaCorrection[correctionIndex]->GetTrack2ParticleCorrection()->GetMeasuredHistogram()->Scale(scaleFactor);
+      fdNdEtaCorrection[correctionIndex]->GetTrack2ParticleCorrection()->GetGeneratedHistogram()->Scale(scaleFactor);
+
+      TString* str = new TString;
+      str->Form("%s%s", (correctionIndex == 0) ? "Pi" : ((correctionIndex == 1) ? "K" : "p"), (index % 2 == 0) ? "Boosted" : "Reduced");
+      return str->Data();
+      break;
+
       // each species rated with pythia ratios
-    case 2: // + 10% pions
-    case 3: // - 10% pions
-    case 4: // + 10% kaons
-    case 5: // - 10% kaons
-    case 6: // + 10% protons
-    case 7: // - 10% protons
+    case 12: // + 50% pions
+    case 13: // - 50% pions
+    case 14: // + 50% kaons
+    case 15: // - 50% kaons
+    case 16: // + 50% protons
+    case 17: // - 50% protons
       TH1** ptDists = DrawRatios("pythiaratios.root");
       Int_t functionIndex = (index - 2) / 2;
-      Double_t scaleFactor = (index % 2 == 0) ? 1.1 : 0.9;
+      Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5;
       ptDists[functionIndex]->Scale(scaleFactor);
 
       for (Int_t i=0; i<3; ++i)
@@ -674,3 +747,43 @@ void drawSystematics()
 
   Sigma2VertexSimulation();
 }
+
+void DrawdNdEtaDifferences()
+{
+  TH1* hists[5];
+
+  TCanvas* canvas = new TCanvas("DrawdNdEtaDifferences", "DrawdNdEtaDifferences", 1000, 500);
+  canvas->Divide(2, 1);
+
+  canvas->cd(1);
+
+  for (Int_t i=0; i<5; ++i)
+  {
+    TFile* file = 0;
+
+    switch(i)
+    {
+      case 0 : file = TFile::Open("systematics_dndeta_reference.root"); break;
+      case 1 : file = TFile::Open("systematics_dndeta_KBoosted.root"); break;
+      case 2 : file = TFile::Open("systematics_dndeta_KReduced.root"); break;
+      case 3 : file = TFile::Open("systematics_dndeta_pBoosted.root"); break;
+      case 4 : file = TFile::Open("systematics_dndeta_pReduced.root"); break;
+      default: return;
+    }
+
+    hists[i] = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2");
+
+    hists[i]->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
+    hists[i]->SetLineColor(i+1);
+    hists[i]->DrawCopy(((i > 0) ? "SAME" : ""));
+  }
+
+  canvas->cd(2);
+
+  for (Int_t i=1; i<5; ++i)
+  {
+    hists[i]->Divide(hists[0]);
+    hists[i]->GetYaxis()->SetRangeUser(0.98, 1.02);
+    hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
+  }
+}
index 04a4b73228075cd59ea835350c22b63e081b3e41..a8d89a45406d3f5dba2d8fc27fcf55c926a9ebf3 100644 (file)
@@ -29,3 +29,52 @@ void makeSystematics(Char_t* dataDir, Int_t nRuns=20, Int_t offset = 0, Bool_t d
 
   chain->Process(selector, option);
 }
+
+void runAnalysisWithDifferentMaps(Char_t* dataDir, Int_t nRuns=20, Int_t offset = 0, Bool_t debug = kFALSE, Bool_t proof = kFALSE)
+{
+  // This functions runs the dN/dEta analysis with different correction maps to gather systematics
+  // It runs with the "normal map", and with 4 other different cases where particle species are enhanced
+  // or reduced.
+  // The normal map is expected in correction_map.root, created by AlidNdEtaCorrectionSelector
+  // The others in new_compositions.root in the folders (K|p)(Boosted|Reduced), created
+  //   by AlidNdEtaSystematicsSelector and Composition() out of drawSystematics.C
+
+  gROOT->ProcessLine(".L testAnalysis2.C");
+
+  gSystem->Exec("rm analysis_esd.root");
+
+  testAnalysis2(dataDir, nRuns, offset, kFALSE, debug, proof, "correction_map.root", "dndeta_correction");
+  if (gSystem->Exec("mv analysis_esd.root systematics_dndeta_reference.root") != 0)
+  {
+    printf("systematics_dndeta_reference.root failed\n");
+    return;
+  }
+
+  testAnalysis2(dataDir, nRuns, offset, kFALSE, debug, proof, "new_compositions.root", "KBoosted");
+  if (gSystem->Exec("mv analysis_esd.root systematics_dndeta_KBoosted.root") != 0)
+  {
+    printf("systematics_dndeta_KBoosted.root failed\n");
+    return;
+  }
+
+  testAnalysis2(dataDir, nRuns, offset, kFALSE, debug, proof, "new_compositions.root", "KReduced");
+  if (gSystem->Exec("mv analysis_esd.root systematics_dndeta_KReduced.root") != 0)
+  {
+    printf("systematics_dndeta_KReduced.root failed\n");
+    return;
+  }
+
+  testAnalysis2(dataDir, nRuns, offset, kFALSE, debug, proof, "new_compositions.root", "pBoosted");
+  if (gSystem->Exec("mv analysis_esd.root systematics_dndeta_pBoosted.root") != 0)
+  {
+    printf("systematics_dndeta_pBoosted.root failed\n");
+    return;
+  }
+
+  testAnalysis2(dataDir, nRuns, offset, kFALSE, debug, proof, "new_compositions.root", "pReduced");
+  if (gSystem->Exec("mv analysis_esd.root systematics_dndeta_pReduced.root") != 0)
+  {
+    printf("systematics_dndeta_pReduced.root failed\n");
+    return;
+  }
+}