adding selector that creates histograms needed for systematic uncertainty
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 24 Jul 2006 17:10:19 +0000 (17:10 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 24 Jul 2006 17:10:19 +0000 (17:10 +0000)
  here: track2particle correction: influence of different particle compositions, influence of secondaries
some other small fixes

14 files changed:
PWG0/PWG0selectorsLinkDef.h
PWG0/dNdEta/AliMultiplicityESDSelector.cxx
PWG0/dNdEta/AlidNdEtaCorrection.cxx
PWG0/dNdEta/AlidNdEtaCorrectionSelector.cxx
PWG0/dNdEta/AlidNdEtaSystematicsSelector.cxx [new file with mode: 0644]
PWG0/dNdEta/AlidNdEtaSystematicsSelector.h [new file with mode: 0644]
PWG0/dNdEta/drawPlots.C
PWG0/dNdEta/drawSystematics.C [new file with mode: 0644]
PWG0/dNdEta/makeCorrection2.C
PWG0/dNdEta/makeSystematics.C [new file with mode: 0644]
PWG0/dNdEta/runMultiplicitySelector.C
PWG0/dNdEta/testAnalysis2.C
PWG0/esdTrackCuts/AliESDtrackCuts.cxx
PWG0/libPWG0selectors.pkg

index 6fea993ed7a2f0090feb28ec9d482de00501c737..ced196a0c72e66a2b40a86d7a3b24f153d1070dc 100644 (file)
@@ -13,5 +13,6 @@
 #pragma link C++ class AlidNdEtaAnalysisESDSelector+;
 #pragma link C++ class AlidNdEtaVertexRecEffSelector+;
 #pragma link C++ class AliMultiplicityESDSelector+;
+#pragma link C++ class AlidNdEtaSystematicsSelector+;
 
 #endif
index b08da1b69d7eea69822be1b69b782ddcb9f75ca2..2016aca6c8d7f810231ee1a7b338697604854980 100644 (file)
@@ -42,6 +42,8 @@ void AliMultiplicityESDSelector::Begin(TTree* tree)
 {
   // Begin function
 
+  AliSelector::Begin(tree);
+
   ReadUserObjects(tree);
 }
 
index 8509f472835d7c66c8df760903d8ee457a67b3eb..1fdd40bf5667b20aaa8d7cb018a34e6109745f79 100644 (file)
@@ -153,15 +153,19 @@ AlidNdEtaCorrection::Merge(TCollection* list) {
     collectionTriggerBias        ->Add(entry->GetTriggerBiasCorrection());
 
     count++;
+
+    fNEvents += entry->fNEvents;
+    fNTriggeredEvents += entry->fNTriggeredEvents;
   }
   fTrack2ParticleCorrection ->Merge(collectionNtrackToNparticle);
   fVertexRecoCorrection        ->Merge(collectionVertexReco);
   fTriggerBiasCorrection        ->Merge(collectionTriggerBias);
-  
+
   delete collectionNtrackToNparticle;
   delete collectionVertexReco;
   delete collectionTriggerBias;
 
+
   return count+1;
 }
 
index 84d96ad27cf27e83f9d4928d567930a212ff5458..7291a40363dddfdba6d15202c307470d66caf3a4 100644 (file)
@@ -213,9 +213,7 @@ Bool_t AlidNdEtaCorrectionSelector::Process(Long64_t entry)
       continue;
 
     if (SignOK(particle->GetPDG()) == kFALSE)
-        continue;
-
-    //if (TMath::Abs(particle->GetPdgCode()) == 2212)      continue;
+      continue;
 
     Float_t eta = particle->Eta();
     Float_t pt = particle->Pt();
@@ -268,8 +266,6 @@ Bool_t AlidNdEtaCorrectionSelector::Process(Long64_t entry)
     if (SignOK(particle->GetPDG()) == kFALSE)
         continue;
 
-    //if (TMath::Abs(particle->GetPdgCode()) == 2212)      continue;
-
     if (vertexReconstructed)
     {
       fdNdEtaCorrection->FillParticleWhenMeasuredTrack(vtxMC[2], particle->Eta(), particle->Pt());
diff --git a/PWG0/dNdEta/AlidNdEtaSystematicsSelector.cxx b/PWG0/dNdEta/AlidNdEtaSystematicsSelector.cxx
new file mode 100644 (file)
index 0000000..719a986
--- /dev/null
@@ -0,0 +1,374 @@
+/* $Id$ */
+
+#include "AlidNdEtaSystematicsSelector.h"
+
+#include <TStyle.h>
+#include <TSystem.h>
+#include <TCanvas.h>
+#include <TVector3.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TH3F.h>
+#include <TH1F.h>
+#include <TParticle.h>
+
+#include <AliLog.h>
+#include <AliESD.h>
+#include <AliGenEventHeader.h>
+#include <AliStack.h>
+#include <AliHeader.h>
+
+#include "esdTrackCuts/AliESDtrackCuts.h"
+#include "AliPWG0Helper.h"
+#include "AlidNdEtaCorrection.h"
+
+ClassImp(AlidNdEtaSystematicsSelector)
+
+AlidNdEtaSystematicsSelector::AlidNdEtaSystematicsSelector() :
+  AliSelectorRL(),
+  fSecondaries(0),
+  fEsdTrackCuts(0)
+{
+  //
+  // Constructor. Initialization of pointers
+  //
+
+  for (Int_t i=0; i<4; ++i)
+    fdNdEtaCorrection[i] = 0;
+}
+
+AlidNdEtaSystematicsSelector::~AlidNdEtaSystematicsSelector()
+{
+  //
+  // Destructor
+  //
+}
+
+void AlidNdEtaSystematicsSelector::Begin(TTree* tree)
+{
+  // Begin function
+
+  AliSelectorRL::Begin(tree);
+
+  ReadUserObjects(tree);
+}
+
+void AlidNdEtaSystematicsSelector::ReadUserObjects(TTree* tree)
+{
+  // read the user objects, called from slavebegin and begin
+
+  if (!fEsdTrackCuts && fInput)
+    fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("AliESDtrackCuts"));
+
+  if (!fEsdTrackCuts && tree)
+    fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (tree->GetUserInfo()->FindObject("AliESDtrackCuts"));
+
+  if (!fEsdTrackCuts)
+     AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list.");
+}
+
+void AlidNdEtaSystematicsSelector::SlaveBegin(TTree* tree)
+{
+  // The SlaveBegin() function is called after the Begin() function.
+  // When running with PROOF SlaveBegin() is called on each slave server.
+  // The tree argument is deprecated (on PROOF 0 is passed).
+
+  AliSelector::SlaveBegin(tree);
+
+  ReadUserObjects(tree);
+
+  TString option(GetOption());
+
+  printf("Running AlidNdEtaSystematicsSelector with options %s\n", option.Data());
+
+  if (option.Contains("secondaries"))
+  {
+    fSecondaries = new TH3F("fSecondaries", "fSecondaries;NacceptedTracks;CratioSecondaries;p_{T}", 201, -0.5, 205.5, 16, 0.45, 2.05, 10, 0, 10);
+  }
+
+  if (option.Contains("particle-composition"))
+  {
+    for (Int_t i=0; i<4; ++i)
+    {
+      TString name;
+      name.Form("correction_%d", i);
+      fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
+    }
+  }
+}
+
+Bool_t AlidNdEtaSystematicsSelector::Process(Long64_t entry)
+{
+  // The Process() function is called for each entry in the tree (or possibly
+  // keyed object in the case of PROOF) to be processed. The entry argument
+  // specifies which entry in the currently loaded tree is to be processed.
+  // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
+  // to read either all or the required parts of the data. When processing
+  // keyed objects with PROOF, the object is already loaded and is available
+  // via the fObject pointer.
+  //
+  // This function should contain the "body" of the analysis. It can contain
+  // simple or elaborate selection criteria, run algorithms on the data
+  // of the event and typically fill histograms.
+
+  // WARNING when a selector is used with a TChain, you must use
+  //  the pointer to the current TTree to call GetEntry(entry).
+  //  The entry is always the local entry number in the current tree.
+  //  Assuming that fTree is the pointer to the TChain being processed,
+  //  use fTree->GetTree()->GetEntry(entry).
+
+  if (AliSelectorRL::Process(entry) == kFALSE)
+    return kFALSE;
+
+  // Check prerequisites
+  if (!fESD)
+  {
+    AliDebug(AliLog::kError, "ESD branch not available");
+    return kFALSE;
+  }
+
+  if (!fEsdTrackCuts)
+  {
+    AliDebug(AliLog::kError, "fESDTrackCuts not available");
+    return kFALSE;
+  }
+
+  AliStack* stack = GetStack();
+  if (!stack)
+  {
+    AliDebug(AliLog::kError, "Stack not available");
+    return kFALSE;
+  }
+
+  TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
+
+  if (fdNdEtaCorrection[0])
+    FillCorrectionMaps(list);
+
+  if (fSecondaries)
+    FillSecondaries(list);
+
+  delete list;
+  list = 0;
+
+  return kTRUE;
+}
+
+void AlidNdEtaSystematicsSelector::FillCorrectionMaps(TObjArray* listOfTracks)
+{
+  // fills the correction maps for different particle species
+
+  AliStack* stack = GetStack();
+  AliHeader* header = GetHeader();
+
+  Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
+  Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
+
+  for (Int_t i=0; i<4; ++i)
+  {
+    fdNdEtaCorrection[i]->IncreaseEventCount();
+    if (eventTriggered)
+      fdNdEtaCorrection[i]->IncreaseTriggeredEventCount();
+  }
+
+  // get the MC vertex
+  AliGenEventHeader* genHeader = header->GenEventHeader();
+
+  TArrayF vtxMC(3);
+  genHeader->PrimaryVertex(vtxMC);
+
+  // loop over mc particles
+  Int_t nPrim  = stack->GetNprimary();
+
+  for (Int_t iMc = 0; iMc < nPrim; ++iMc)
+  {
+    TParticle* particle = stack->Particle(iMc);
+
+    if (!particle)
+    {
+      AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
+      continue;
+    }
+
+    if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
+      continue;
+
+    Float_t eta = particle->Eta();
+    Float_t pt = particle->Pt();
+
+    Int_t id = -1;
+    switch (TMath::Abs(particle->GetPdgCode()))
+    {
+      case 211: id = 0; break;
+      case 321: id = 1; break;
+      case 2212: id = 2; break;
+      default: id = 3; break;
+    }
+
+    if (vertexReconstructed)
+      fdNdEtaCorrection[id]->FillParticle(vtxMC[2], eta, pt);
+
+    fdNdEtaCorrection[id]->FillParticleAllEvents(eta, pt);
+    if (eventTriggered)
+      fdNdEtaCorrection[id]->FillParticleWhenEventTriggered(eta, pt);
+  }// end of mc particle
+
+  // loop over esd tracks
+  TIterator* iter = listOfTracks->MakeIterator();
+  TObject* obj = 0;
+  while ((obj = iter->Next()) != 0)
+  {
+    AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
+    if (!esdTrack)
+      continue;
+
+    // using the properties of the mc particle
+    Int_t label = TMath::Abs(esdTrack->GetLabel());
+    if (label == 0)
+    {
+      AliDebug(AliLog::kWarning, Form("WARNING: cannot find corresponding mc part for track %d.", label));
+      continue;
+    }
+
+    TParticle* particle = stack->Particle(label);
+    if (!particle)
+    {
+      AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
+      continue;
+    }
+
+    Int_t id = -1;
+    switch (TMath::Abs(particle->GetPdgCode()))
+    {
+      case 211: id = 0; break;
+      case 321: id = 1; break;
+      case 2212: id = 2; break;
+      default: id = 3; break;
+    }
+
+    if (vertexReconstructed)
+      fdNdEtaCorrection[id]->FillParticleWhenMeasuredTrack(vtxMC[2], particle->Eta(), particle->Pt());
+  } // end of track loop
+
+  delete iter;
+  iter = 0;
+}
+
+void AlidNdEtaSystematicsSelector::FillSecondaries(TObjArray* listOfTracks)
+{
+  // 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());
+
+  TIterator* iter = listOfTracks->MakeIterator();
+  TObject* obj = 0;
+  while ((obj = iter->Next()) != 0)
+  {
+    AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
+    if (!esdTrack)
+      continue;
+
+    Int_t label = TMath::Abs(esdTrack->GetLabel());
+    if (label == 0)
+    {
+      AliDebug(AliLog::kWarning, Form("WARNING: cannot find corresponding mc part for track %d.", label));
+      continue;
+    }
+
+    TParticle* particle = stack->Particle(label);
+    if (!particle)
+    {
+      AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
+      continue;
+    }
+
+    Int_t nPrim  = stack->GetNprimary();
+    if (label < nPrim)
+      nPrimaries->Fill(particle->Pt());
+    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);
+
+    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)));
+
+      for (Double_t factor = 0.5; factor < 2.01; factor += 0.1)
+        fSecondaries->Fill((Double_t) primaries + (Double_t) secondaries * factor, factor, nPrimaries->GetBinCenter(i));
+    }
+  }
+
+  delete nPrimaries;
+  nPrimaries = 0;
+
+  delete nSecondaries;
+  nSecondaries = 0;
+
+  delete iter;
+  iter = 0;
+}
+
+void AlidNdEtaSystematicsSelector::SlaveTerminate()
+{
+  // The SlaveTerminate() function is called after all entries or objects
+  // have been processed. When running with PROOF SlaveTerminate() is called
+  // on each slave server.
+
+  AliSelector::SlaveTerminate();
+
+  // Add the histograms to the output on each slave server
+  if (!fOutput)
+  {
+    AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
+    return;
+  }
+
+  if (fSecondaries)
+    fOutput->Add(fSecondaries);
+
+  for (Int_t i=0; i<4; ++i)
+    if (fdNdEtaCorrection[i])
+      fOutput->Add(fdNdEtaCorrection[i]);
+}
+
+void AlidNdEtaSystematicsSelector::Terminate()
+{
+  // The Terminate() function is the last function to be called during
+  // a query. It always runs on the client, it can be used to present
+  // the results graphically or save the results to file.
+
+  AliSelector::Terminate();
+
+  fSecondaries = dynamic_cast<TH3F*> (fOutput->FindObject("fSecondaries"));
+  for (Int_t i=0; i<4; ++i)
+    fdNdEtaCorrection[i] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
+
+  TFile* fout = TFile::Open("systematics.root", "RECREATE");
+
+  if (fEsdTrackCuts)
+    fEsdTrackCuts->SaveHistograms("esd_track_cuts");
+
+  if (fSecondaries)
+    fSecondaries->Write();
+
+  for (Int_t i=0; i<4; ++i)
+    if (fdNdEtaCorrection[i])
+      fdNdEtaCorrection[i]->SaveHistograms();
+
+  fout->Write();
+  fout->Close();
+
+  if (fSecondaries)
+  {
+    new TCanvas;
+    fSecondaries->Draw();
+  }
+}
diff --git a/PWG0/dNdEta/AlidNdEtaSystematicsSelector.h b/PWG0/dNdEta/AlidNdEtaSystematicsSelector.h
new file mode 100644 (file)
index 0000000..29aef2f
--- /dev/null
@@ -0,0 +1,39 @@
+/* $Id$ */
+
+#ifndef ALIDNDETASYSTEMATICSSELECTOR_H
+#define ALIDNDETASYSTEMATICSSELECTOR_H
+
+#include "AliSelectorRL.h"
+
+class AliESDtrackCuts;
+class AlidNdEtaCorrection;
+class TH3F;
+
+class AlidNdEtaSystematicsSelector : public AliSelectorRL {
+  public:
+    AlidNdEtaSystematicsSelector();
+    virtual ~AlidNdEtaSystematicsSelector();
+
+    virtual void    Begin(TTree* tree);
+    virtual void    SlaveBegin(TTree *tree);
+    virtual Bool_t  Process(Long64_t entry);
+    virtual void    SlaveTerminate();
+    virtual void    Terminate();
+
+ protected:
+    void ReadUserObjects(TTree* tree);
+
+    void FillCorrectionMaps(TObjArray* listOfTracks);
+    void FillSecondaries(TObjArray* listOfTracks);
+
+    TH3F* fSecondaries; // (accepted tracks) vs (tracks from sec)/(n * tracks from sec) vs pT
+    AlidNdEtaCorrection* fdNdEtaCorrection[4];      // correction for different particle species: here pi, K, p, others
+
+    AliESDtrackCuts* fEsdTrackCuts;     // Object containing the parameters of the esd track cuts
+
+ private:
+
+  ClassDef(AlidNdEtaSystematicsSelector, 0);
+};
+
+#endif
index 1aa110b9a0ecd90f3333d6f03793a601049ad70f..29e9bbb2cf639e6b67150885629895d4232ab476 100644 (file)
@@ -434,7 +434,7 @@ void Track2Particle1D(const char* fileName = "correction_map.root", Float_t uppe
 
   corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
 
-  TCanvas* canvas = new TCanvas("Track2Particle2D", "Track2Particle2D", 1200, 400);
+  TCanvas* canvas = new TCanvas("Track2Particle1D", "Track2Particle1D", 1200, 400);
   canvas->Divide(3, 1);
 
   canvas->cd(1);
diff --git a/PWG0/dNdEta/drawSystematics.C b/PWG0/dNdEta/drawSystematics.C
new file mode 100644 (file)
index 0000000..e140a12
--- /dev/null
@@ -0,0 +1,176 @@
+/* $Id$ */
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include "AliPWG0Helper.h"
+#include "dNdEtaAnalysis.h"
+#include "AlidNdEtaCorrection.h"
+
+#include <TCanvas.h>
+#include <TFile.h>
+#include <TH1.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <TLine.h>
+#include <TSystem.h>
+
+#endif
+
+void SetRanges(TAxis* axis)
+{
+  if (strcmp(axis->GetTitle(), "#eta") == 0)
+    axis->SetRangeUser(-1.7999, 1.7999);
+  if (strcmp(axis->GetTitle(), "p_{T} [GeV/c]") == 0)
+    axis->SetRangeUser(0, 9.9999);
+  if (strcmp(axis->GetTitle(), "vtx z [cm]") == 0)
+    axis->SetRangeUser(-15, 14.9999);
+  if (strcmp(axis->GetTitle(), "Ntracks") == 0)
+    axis->SetRangeUser(0, 99.9999);
+}
+
+void SetRanges(TH1* hist)
+{
+  SetRanges(hist->GetXaxis());
+  SetRanges(hist->GetYaxis());
+  SetRanges(hist->GetZaxis());
+}
+
+void Prepare3DPlot(TH3* hist)
+{
+  hist->GetXaxis()->SetTitleOffset(1.5);
+  hist->GetYaxis()->SetTitleOffset(1.5);
+  hist->GetZaxis()->SetTitleOffset(1.5);
+
+  hist->SetStats(kFALSE);
+}
+
+void Prepare2DPlot(TH2* hist)
+{
+  hist->SetStats(kFALSE);
+  hist->GetYaxis()->SetTitleOffset(1.4);
+
+  SetRanges(hist);
+}
+
+void Prepare1DPlot(TH1* hist)
+{
+  hist->SetLineWidth(2);
+  hist->SetStats(kFALSE);
+
+  SetRanges(hist);
+}
+
+void InitPad()
+{
+  gPad->Range(0, 0, 1, 1);
+  gPad->SetLeftMargin(0.15);
+  //gPad->SetRightMargin(0.05);
+  //gPad->SetTopMargin(0.13);
+  //gPad->SetBottomMargin(0.1);
+
+  //gPad->SetGridx();
+  //gPad->SetGridy();
+}
+
+void InitPadCOLZ()
+{
+  gPad->Range(0, 0, 1, 1);
+  gPad->SetRightMargin(0.15);
+  gPad->SetLeftMargin(0.12);
+
+  gPad->SetGridx();
+  gPad->SetGridy();
+}
+
+void Secondaries()
+{
+  TFile* file = TFile::Open("systematics.root");
+
+  TH3F* secondaries = dynamic_cast<TH3F*> (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++)
+  {
+    TGraph* graph = new TGraph;
+    graph->Clear();
+    graph->SetTitle(Form("%f < p_{T} < %f", secondaries->GetZaxis()->GetBinLowEdge(ptBin), secondaries->GetZaxis()->GetBinUpEdge(ptBin)));
+
+    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)
+        graph->SetPoint(graph->GetN(), secondaries->GetYaxis()->GetBinCenter(cBin), sum / count);
+    }
+
+    new TCanvas;
+    graph->SetMarkerStyle(21);
+    graph->Draw("AP");
+    graph->Print();
+  }
+}
+
+void Composition()
+{
+  gSystem->Load("libPWG0base");
+
+  AlidNdEtaCorrection* fdNdEtaCorrection[4];
+
+  TFile::Open("systematics.root");
+
+  for (Int_t i=0; i<4; ++i)
+  {
+    TString name;
+    name.Form("correction_%d", i);
+    fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
+    fdNdEtaCorrection[i]->LoadHistograms("systematics.root", name);
+  }
+
+  //fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetMeasuredHistogram()->Scale(2);
+  //fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetGeneratedHistogram()->Scale(2);
+
+  AlidNdEtaCorrection* finalCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
+
+  TList* collection = new TList;
+
+  for (Int_t i=0; i<4; ++i)
+    collection->Add(fdNdEtaCorrection[i]);
+
+  finalCorrection->Merge(collection);
+
+  delete collection;
+
+  finalCorrection->Finish();
+
+  TFile* file = TFile::Open("temp.root", "RECREATE");
+  finalCorrection->SaveHistograms();
+  file->Write();
+  file->Close();
+
+  gROOT->ProcessLine(".L drawPlots.C");
+  Track2Particle1D("temp.root");
+}
+
+void drawSystematics()
+{
+  Composition();
+}
index d175c438cba3a903f1195084afa97c804bf6db12..1913981aa45e8d31035e0a20d0e76ca70a4db37d 100644 (file)
@@ -14,25 +14,19 @@ void makeCorrection2(Char_t* dataDir, Int_t nRuns=20, Int_t offset = 0, Bool_t d
   gSystem->Load("libPWG0base");
   gSystem->Load("libPWG0dep");
 
-  TChain* chain = CreateESDChain(dataDir, nRuns, offset);
-
-  fEsdTrackCuts = new AliESDtrackCuts();
-  fEsdTrackCuts->DefineHistograms(1);
-
-  fEsdTrackCuts->SetMinNClustersTPC(50);
-  fEsdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
-  fEsdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
-  fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
+  gROOT->ProcessLine(".L CreateCuts.C");
 
-  fEsdTrackCuts->SetMinNsigmaToVertex(3);
-  fEsdTrackCuts->SetAcceptKingDaughters(kFALSE);
+  AliESDtrackCuts* esdTrackCuts = CreateTrackCuts();
+  if (!esdTrackCuts)
+  {
+    printf("ERROR: esdTrackCuts could not be created\n");
+    return;
+  }
 
-  chain->GetUserInfo()->Add(fEsdTrackCuts);
-
-  AliLog::SetClassDebugLevel("AlidNdEtaCorrectionSelector", AliLog::kInfo);
-  AliLog::SetClassDebugLevel("AliSelectorRL", AliLog::kInfo);
+  TChain* chain = CreateESDChain(dataDir, nRuns, offset);
+  chain->GetUserInfo()->Add(esdTrackCuts);
 
-  TString selector("AlidNdEtaCorrectionSelector.cxx+");
+  TString selector("AlidNdEtaCorrectionSelector.cxx++");
   if (debug != kFALSE)
     selector += "g";
 
diff --git a/PWG0/dNdEta/makeSystematics.C b/PWG0/dNdEta/makeSystematics.C
new file mode 100644 (file)
index 0000000..04a4b73
--- /dev/null
@@ -0,0 +1,31 @@
+/* $Id$ */
+
+//
+// Script to run the creation of input for systematics
+//
+
+#include "../CreateESDChain.C"
+
+void makeSystematics(Char_t* dataDir, Int_t nRuns=20, Int_t offset = 0, Bool_t debug = kFALSE, const Char_t* option = "")
+{
+  gSystem->Load("libPWG0base");
+  gSystem->Load("libPWG0dep");
+
+  gROOT->ProcessLine(".L CreateCuts.C");
+
+  AliESDtrackCuts* esdTrackCuts = CreateTrackCuts();
+  if (!esdTrackCuts)
+  {
+    printf("ERROR: esdTrackCuts could not be created\n");
+    return;
+  }
+
+  TChain* chain = CreateESDChain(dataDir, nRuns, offset);
+  chain->GetUserInfo()->Add(esdTrackCuts);
+
+  TString selector("AlidNdEtaSystematicsSelector.cxx+");
+  if (debug != kFALSE)
+    selector += "g";
+
+  chain->Process(selector, option);
+}
index b271a12c25826a30ade10e5116a4b74ab14aa222..98701acc326e13aa09d64a77c29005ff9328cede 100644 (file)
@@ -8,16 +8,13 @@
 
 void runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kFALSE, Bool_t aDebug = kFALSE, Bool_t aProof = kFALSE)
 {
-  /*gSystem->Load("libEG");
+  TStopwatch timer;
+  timer.Start();
+
+  gSystem->Load("libEG");
   gSystem->Load("libGeom");
-  gSystem->Load("libESD");*/
-  gSystem->cd("ESD");
-  gROOT->ProcessLine(".x PROOF-INF/SETUP.C");
-  gSystem->cd("..");
-  gSystem->cd("PWG0base");
-  gROOT->ProcessLine(".x PROOF-INF/SETUP.C");
-  gSystem->cd("..");
-  //gSystem->Load("libPWG0base");
+  gSystem->Load("libESD");
+  gSystem->Load("libPWG0base");
   if (aMC != kFALSE)
     gSystem->Load("libPWG0dep");
 
@@ -30,7 +27,7 @@ void runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_
 
   if (aProof != kFALSE)
   {
-    proof = TProof::Open("jgrosseo@lxb6043");
+    proof = TProof::Open("jgrosseo@lxb6046");
 
     if (!proof)
     {
@@ -83,9 +80,6 @@ void runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_
   if (aDebug != kFALSE)
     selectorName += "g";
 
-  TStopwatch timer;
-  timer.Start();
-
   Long64_t result = -1;
 
   if (proof != kFALSE)
index 711090205aa7e5b541490bd4878bc591d155c634..fe55cb3959a76767627adc597bada5170a00659f 100644 (file)
@@ -21,14 +21,13 @@ void testAnalysis2(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kF
 
   gROOT->ProcessLine(".L CreatedNdEta.C");
   gROOT->ProcessLine(".L CreateCuts.C");
+  gROOT->ProcessLine(".L drawPlots.C");
 
-  TChain* chain = 0;
+  TChain* chain = CreateESDChain(data, nRuns, offset);;
   TVirtualProof* proof = 0;
-  if (aProof == kFALSE)
-    chain = CreateESDChainFromDir(data, nRuns, offset);
-  else
+
+  if (aProof != kFALSE)
   {
-    chain = CreateESDChainFromList(data, nRuns, offset);
     proof = TProof::Open("jgrosseo@lxb6046");
 
     if (!proof)
@@ -70,30 +69,26 @@ void testAnalysis2(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kF
     return;
   }
 
-  //chain->GetUserInfo()->Add(esdTrackCuts);
-  //if (proof)
-  //  proof->AddInput(esdTrackCuts);
+  chain->GetUserInfo()->Add(esdTrackCuts);
+  if (proof)
+    proof->AddInput(esdTrackCuts);
 
   if (aMC == kFALSE)
   {
     AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
     dNdEtaCorrection->LoadHistograms("correction_map.root","dndeta_correction");
+    dNdEtaCorrection->ReduceInformation();
     //dNdEtaCorrection->RemoveEdges(2, 0, 2);
 
-    //chain->GetUserInfo()->Add(dNdEtaCorrection);
-    //if (proof)
-    //  proof->AddInput(dNdEtaCorrection);
+    chain->GetUserInfo()->Add(dNdEtaCorrection);
+    if (proof)
+      proof->AddInput(dNdEtaCorrection);
   }
 
   TString selectorName = ((aMC == kFALSE) ? "AlidNdEtaAnalysisESDSelector" : "AlidNdEtaAnalysisMCSelector");
   AliLog::SetClassDebugLevel(selectorName, AliLog::kInfo);
 
-  // workaround for a bug in PROOF that only allows header files for .C files
-  // please create symlink from <selector>.cxx to <selector>.C
-  if (proof != kFALSE)
-    selectorName += ".C++";
-  else
-    selectorName += ".cxx++";
+  selectorName += ".cxx++";
 
   if (aDebug != kFALSE)
     selectorName += "g";
@@ -117,6 +112,17 @@ void testAnalysis2(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kF
   timer.Stop();
   timer.Print();
 
-  CreatedNdEta(aMC ? kFALSE : kTRUE, aMC ? "analysis_mc.root" : "analysis_esd.root");
+  dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
+
+  TFile* file = TFile::Open(aMC ? "analysis_mc.root" : "analysis_esd.root");
+  if (!file)
+  {
+    cout << "Error. File not found" << endl;
+    return;
+  }
+  fdNdEtaAnalysis->LoadHistograms();
+  fdNdEtaAnalysis->DrawHistograms();
+
+  dNdEta(kTRUE);
 }
 
index 98ead8e1349c800487d0af42dabdb4851fbb2958..de124fd8209c17473b4975a7249d62ffb2669273 100644 (file)
@@ -303,8 +303,12 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
   Float_t nSigmaToVertex = -1;
   if (bRes[0]!=0 && bRes[1]!=0) {
     Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
-    nSigmaToVertex = d;//TMath::Sqrt(2)*(TMath::ErfInverse(1 - TMath::Exp(0.5*(-d*d))));
-    // JF solution: nSigmaToVertex = TMath::ErfInverse(TMath::Sqrt(2.0/TMath::Pi()) * TMath::Erf(d / TMath::Sqrt(2))) * TMath::Sqrt(2);
+    nSigmaToVertex = d;
+
+    // JF solution:
+    //nSigmaToVertex = TMath::ErfInverse(TMath::Sqrt(2.0/TMath::Pi()) * TMath::Erf(d / TMath::Sqrt(2))) * TMath::Sqrt(2);
+    // Claus solution:
+    //nSigmaToVertex = TMath::Sqrt(2)*(TMath::ErfInverse(1 - TMath::Exp(0.5*(-d*d))));
   }
 
   // getting the kinematic variables of the track 
@@ -357,7 +361,7 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
   if (nSigmaToVertex > fCutNsigmaToVertex) 
     cuts[11] = kTRUE;
   // if n sigma could not be calculated
-  if (nSigmaToVertex<0 && fCutSigmaToVertexRequired)   
+  if (nSigmaToVertex<0 && fCutSigmaToVertexRequired)
     cuts[12]=kTRUE;
   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0) 
     cuts[13]=kTRUE;
@@ -419,7 +423,7 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
     fhDXYvsDZ[0]->Fill(b[1],b[0]);
 
     if (bRes[0]!=0 && bRes[1]!=0) {
-      fhDZNormalized[0]->Fill(b[1]/bRes[1]);     
+      fhDZNormalized[0]->Fill(b[1]/bRes[1]);
       fhDXYNormalized[0]->Fill(b[0]/bRes[0]);    
       fhDXYvsDZNormalized[0]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
     }
@@ -439,7 +443,7 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
     
     fhC11[1]->Fill(extCov[0]);                 
     fhC22[1]->Fill(extCov[2]);                 
-    fhC33[1]->Fill(extCov[5]);                 
+    fhC33[1]->Fill(extCov[5]);
     fhC44[1]->Fill(extCov[9]);                                  
     fhC55[1]->Fill(extCov[14]);                                  
     
@@ -447,9 +451,12 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
     fhDXY[1]->Fill(b[0]);    
     fhDXYvsDZ[1]->Fill(b[1],b[0]);
 
-    fhDZNormalized[1]->Fill(b[1]/bRes[1]);     
-    fhDXYNormalized[1]->Fill(b[0]/bRes[0]);    
-    fhDXYvsDZNormalized[1]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
+    if (bRes[0]!=0 && bRes[1]!=0)
+    {
+      fhDZNormalized[1]->Fill(b[1]/bRes[1]);
+      fhDXYNormalized[1]->Fill(b[0]/bRes[0]);
+      fhDXYvsDZNormalized[1]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
+    }
   }
   
   return kTRUE;
index 21cfb4664527bd2342b9231ac179d25f443284a8..099367c3c06a172ebbcce4903c964d47fe06f353 100644 (file)
@@ -7,7 +7,8 @@ HDRS = dNdEta/AlidNdEtaCorrectionSelector.h \
        dNdEta/AlidNdEtaAnalysisMCSelector.h \
        dNdEta/AlidNdEtaAnalysisESDSelector.h \
        dNdEta/AlidNdEtaVertexRecEffSelector.h \
-       dNdEta/AliMultiplicityESDSelector.h
+       dNdEta/AliMultiplicityESDSelector.h \
+       dNdEta/AlidNdEtaSystematicsSelector.h
 
 SRCS = $(HDRS:.h=.cxx)