implemented analysis using 3d information (vtx_z, eta, pt)
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 13 Jun 2006 15:39:31 +0000 (15:39 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 13 Jun 2006 15:39:31 +0000 (15:39 +0000)
12 files changed:
PWG0/AliCorrectionMatrix3D.cxx
PWG0/AliCorrectionMatrix3D.h
PWG0/AliPWG0Helper.cxx
PWG0/AliPWG0Helper.h
PWG0/dNdEta/AlidNdEtaAnalysisESDSelector.cxx
PWG0/dNdEta/AlidNdEtaCorrection.cxx
PWG0/dNdEta/AlidNdEtaCorrection.h
PWG0/dNdEta/AlidNdEtaCorrectionSelector.cxx
PWG0/dNdEta/CreatedNdEta.C
PWG0/dNdEta/dNdEtaAnalysis.cxx
PWG0/dNdEta/dNdEtaAnalysis.h
PWG0/dNdEta/testAnalysis2.C

index 48d15b2..b856ae6 100644 (file)
@@ -12,6 +12,7 @@
 #include <AliLog.h>
 
 #include "AliCorrectionMatrix3D.h"
+#include "AliPWG0Helper.h"
 
 //____________________________________________________________________
 ClassImp(AliCorrectionMatrix3D)
@@ -119,27 +120,9 @@ void AliCorrectionMatrix3D::SaveHistograms()
 
   AliCorrectionMatrix::SaveHistograms();
 
-  WriteProjections(GetMeasuredHistogram());
-  WriteProjections(GetGeneratedHistogram());
+  AliPWG0Helper::CreateProjections(GetMeasuredHistogram());
+  AliPWG0Helper::CreateProjections(GetGeneratedHistogram());
 
   if (GetCorrectionHistogram())
-    WriteProjections(GetCorrectionHistogram());
-}
-
-//____________________________________________________________________
-void AliCorrectionMatrix3D::WriteProjections(TH3F* hist)
-{
-  // write some projections to disk
-
-  TH1* proj = hist->Project3D("yx");
-  proj->SetXTitle(hist->GetXaxis()->GetTitle());
-  proj->SetYTitle(hist->GetYaxis()->GetTitle());
-
-  proj = hist->Project3D("zx");
-  proj->SetXTitle(hist->GetXaxis()->GetTitle());
-  proj->SetYTitle(hist->GetZaxis()->GetTitle());
-
-  proj = hist->Project3D("zy");
-  proj->SetXTitle(hist->GetYaxis()->GetTitle());
-  proj->SetYTitle(hist->GetZaxis()->GetTitle());
+    AliPWG0Helper::CreateProjections(GetCorrectionHistogram());
 }
index 0ef3b57..f6619f6 100644 (file)
@@ -37,7 +37,6 @@ public:
   void RemoveEdges(Float_t cut=2, Int_t nBinsXedge = 0, Int_t nBinsYedge = 0, Int_t nBinsZedge = 0);
 
   virtual void SaveHistograms();
-  void WriteProjections(TH3F* hist);
 
 protected:
   ClassDef(AliCorrectionMatrix3D,1)
index 98420bc..9f036a0 100644 (file)
@@ -4,6 +4,8 @@
 
 #include <TParticle.h>
 #include <TParticlePDG.h>
+#include <TH1.h>
+#include <TH3F.h>
 
 #include <AliLog.h>
 #include <AliESD.h>
@@ -98,3 +100,22 @@ Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimari
 
   return kTRUE;
 }
+
+//____________________________________________________________________
+void AliPWG0Helper::CreateProjections(TH3F* hist)
+{
+  // create projections of 3d hists to all 2d combinations
+  // the histograms are not returned, just use them from memory or use this to create them in a file
+
+  TH1* proj = hist->Project3D("yx");
+  proj->SetXTitle(hist->GetXaxis()->GetTitle());
+  proj->SetYTitle(hist->GetYaxis()->GetTitle());
+
+  proj = hist->Project3D("zx");
+  proj->SetXTitle(hist->GetXaxis()->GetTitle());
+  proj->SetYTitle(hist->GetZaxis()->GetTitle());
+
+  proj = hist->Project3D("zy");
+  proj->SetXTitle(hist->GetYaxis()->GetTitle());
+  proj->SetYTitle(hist->GetZaxis()->GetTitle());
+}
index 404d454..3e79132 100644 (file)
@@ -9,6 +9,7 @@
 
 class AliESD;
 class TParticle;
+class TH3F;
 
 class AliPWG0Helper : public TObject
 {
@@ -17,6 +18,8 @@ class AliPWG0Helper : public TObject
     static Bool_t IsVertexReconstructed(AliESD* aEsd);
     static Bool_t IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimaries);
 
+    static void CreateProjections(TH3F* hist);
+
   protected:
     ClassDef(AliPWG0Helper, 0)
 };
index bca1f7f..e424e3f 100644 (file)
@@ -29,7 +29,7 @@ AlidNdEtaAnalysisESDSelector::AlidNdEtaAnalysisESDSelector() :
   // Constructor. Initialization of pointers
   //
 
-  AliLog::SetClassDebugLevel("AlidNdEtaAnalysisESDSelector", AliLog::kDebug);
+  //AliLog::SetClassDebugLevel("AlidNdEtaAnalysisESDSelector", AliLog::kDebug);
 }
 
 AlidNdEtaAnalysisESDSelector::~AlidNdEtaAnalysisESDSelector()
@@ -123,6 +123,9 @@ Bool_t AlidNdEtaAnalysisESDSelector::Process(Long64_t entry)
     return kFALSE;
   }
 
+  if (AliPWG0Helper::IsEventTriggered(fESD) == kFALSE)
+    return kTRUE;
+
   if (AliPWG0Helper::IsVertexReconstructed(fESD) == kFALSE)
     return kTRUE;
 
@@ -137,6 +140,12 @@ Bool_t AlidNdEtaAnalysisESDSelector::Process(Long64_t entry)
   Int_t nGoodTracks = list->GetEntries();
 
   Float_t vertexRecoCorr = fdNdEtaCorrection->GetVertexRecoCorrection(vtx[2], nGoodTracks);
+  if (vertexRecoCorr <= 0)
+  {
+    AliDebug(AliLog::kError, Form("INFO: Skipping event because vertexRecoCorr is <= 0 (%f)", vertexRecoCorr));
+    delete list;
+    return kTRUE;
+  }
 
   // loop over esd tracks
   for (Int_t t=0; t<nGoodTracks; t++)
@@ -156,12 +165,16 @@ Bool_t AlidNdEtaAnalysisESDSelector::Process(Long64_t entry)
     Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
     Float_t pt = vector.Pt();
 
-    // TODO pt cut
-
     Float_t track2particleCorr = fdNdEtaCorrection->GetTrack2ParticleCorrection(vtx[2], eta, pt);
 
-    fdNdEtaAnalysis->FillTrack(vtx[2], eta, pt, vertexRecoCorr * track2particleCorr);
+    Float_t weight = vertexRecoCorr * track2particleCorr;
+    if (weight <= 0)
+    {
+      AliDebug(AliLog::kError, Form("INFO: Skipping track because weight is <= 0 (track %d, weight %f)", t, weight));
+      continue;
+    }
 
+    fdNdEtaAnalysis->FillTrack(vtx[2], eta, pt, weight);
   } // end of track loop
 
   delete list;
index a69b244..6256bcc 100644 (file)
@@ -3,7 +3,8 @@
 #include "AlidNdEtaCorrection.h"
 
 #include <TCanvas.h>
-#include <TH2F.h>
+#include <TH3F.h>
+#include <TH1D.h>
 
 //____________________________________________________________________
 ClassImp(AlidNdEtaCorrection)
@@ -140,3 +141,49 @@ void AlidNdEtaCorrection::DrawHistograms()
   fTriggerBiasCorrection       ->DrawHistograms();
 
 }
+
+//____________________________________________________________________
+Float_t AlidNdEtaCorrection::GetMeasuredFraction(Float_t ptCutOff, Float_t eta, Bool_t debug)
+{
+  // calculates the fraction of particles measured (some are missed due to the pt cut off)
+  // uses the generated particle histogram from fTrack2ParticleCorrection
+
+  TH3F* generated = fTrack2ParticleCorrection->GetGeneratedHistogram();
+
+  // find eta borders, if eta is negative assume -0.8 ... 0.8
+  Int_t etaBegin = 0;
+  Int_t etaEnd = 0;
+  if (eta < 0)
+  {
+    etaBegin = generated->GetYaxis()->FindBin(-0.8);
+    etaEnd = generated->GetYaxis()->FindBin(0.8);
+  }
+  else
+  {
+    etaBegin = generated->GetYaxis()->FindBin(eta);
+    etaEnd = etaBegin;
+  }
+
+  Int_t vertexBegin = generated->GetXaxis()->FindBin(-10);
+  Int_t vertexEnd = generated->GetXaxis()->FindBin(10);
+
+  TH1D* ptProj = dynamic_cast<TH1D*> (generated->ProjectionZ(Form("%s_pt", GetName()), vertexBegin, vertexEnd, etaBegin, etaEnd));
+
+  Int_t ptBin = ptProj->FindBin(ptCutOff);
+  Float_t abovePtCut = ptProj->Integral(ptBin, ptProj->GetNbinsX());
+  Float_t all = ptProj->Integral();
+
+  if (all == 0)
+    return -1;
+
+  Float_t fraction = abovePtCut / all;
+
+  if (debug)
+  {
+    new TCanvas;
+    ptProj->Draw();
+  }
+
+  return fraction;
+}
+
index a581e80..8c754b2 100644 (file)
@@ -61,11 +61,13 @@ public:
   
   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);}
 
   Float_t GetTriggerBiasCorrection(Float_t eta, Float_t pt=0) {return fTriggerBiasCorrection->GetCorrection(eta, pt);}
 
+  Float_t GetMeasuredFraction(Float_t ptCutOff, Float_t eta = -1, Bool_t debug = kFALSE);
+
 protected:
   AliCorrectionMatrix3D* fTrack2ParticleCorrection; // handles the track-to-particle correction, function of vtx_z, eta, pt
   AliCorrectionMatrix2D* fVertexRecoCorrection;     // handles the vertex reconstruction efficiency, function of n_clustersITS and vtx_z
index 8abc3c6..64b46cb 100644 (file)
@@ -134,7 +134,6 @@ Bool_t AlidNdEtaCorrectionSelector::Process(Long64_t entry)
 
   Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
 
-  // check if the event was triggered
   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
 
   fdNdEtaCorrection->IncreaseEventCount();
@@ -214,9 +213,13 @@ Bool_t AlidNdEtaCorrectionSelector::Process(Long64_t entry)
       fdNdEtaCorrection->FillParticleWhenMeasuredTrack(vtxMC[2], particle->Eta(), particle->Pt());
   } // end of track loop
 
-  fdNdEtaCorrection->FillEvent(vtxMC[2], nGoodTracks);
-  if (vertexReconstructed)
-    fdNdEtaCorrection->FillEventWithReconstructedVertex(vtxMC[2], nGoodTracks);
+
+  if (eventTriggered)
+  {
+    fdNdEtaCorrection->FillEvent(vtxMC[2], nGoodTracks);
+    if (vertexReconstructed)
+      fdNdEtaCorrection->FillEventWithReconstructedVertex(vtxMC[2], nGoodTracks);
+  }
 
   return kTRUE;
 }
index c7b6fd2..ce97337 100644 (file)
@@ -4,12 +4,12 @@ void CreatedNdEta(Bool_t correct = kTRUE, const Char_t* filename = "analysis_esd
 {
   gSystem->Load("libPWG0base");
 
-  dNdEtaCorrection* dNdEtaCorrection = 0;
+  AlidNdEtaCorrection* dNdEtaCorrection = 0;
   if (correct)
   {
-    dNdEtaCorrection = new dNdEtaCorrection();
+    dNdEtaCorrection = new AlidNdEtaCorrection();
     dNdEtaCorrection->LoadHistograms("correction_map.root","dndeta_correction");
-    dNdEtaCorrection->RemoveEdges(2, 0, 2);
+    //dNdEtaCorrection->RemoveEdges(2, 0, 2);
   }
 
   fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
@@ -22,7 +22,7 @@ void CreatedNdEta(Bool_t correct = kTRUE, const Char_t* filename = "analysis_esd
   }
   fdNdEtaAnalysis->LoadHistograms();
 
-  fdNdEtaAnalysis->Finish(dNdEtaCorrection);
+  fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3);
 
   fdNdEtaAnalysis->DrawHistograms();
 }
index 3c86850..fad1fa0 100644 (file)
@@ -4,6 +4,7 @@
 
 #include <TFile.h>
 #include <TH3F.h>
+#include <TH2D.h>
 #include <TH1D.h>
 #include <TMath.h>
 #include <TCanvas.h>
@@ -13,6 +14,7 @@
 #include <TLegend.h>
 
 #include "AlidNdEtaCorrection.h"
+#include "AliPWG0Helper.h"
 
 //____________________________________________________________________
 ClassImp(dNdEtaAnalysis)
@@ -27,7 +29,7 @@ dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title) :
 {
   // constructor
 
-  fData  = new TH3F(Form("%s_analysis", name),"",80,-20,20,120,-6,6,100, 0, 10);
+  fData  = new TH3F(Form("%s_analysis", name),"dNdEtaAnalysis",80,-20,20,120,-6,6,100, 0, 10);
   fData->SetXTitle("vtx z [cm]");
   fData->SetYTitle("#eta");
   fData->SetZTitle("p_{T}");
@@ -128,13 +130,13 @@ void dNdEtaAnalysis::FillEvent(Float_t vtx, Float_t weight)
 {
   // fills an event into the histograms
 
-  fVtx->Fill(vtx, weight); // TODO vtx distribution with or without weight?
+  fVtx->Fill(vtx, weight);
 
   fNEvents += weight;
 }
 
 //____________________________________________________________________
-void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction)
+void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut)
 {
   // correct with correction values if available
 
@@ -142,30 +144,21 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction)
   if (!correction)
     printf("INFO: No correction applied\n");
 
-  // this can be replaced by TH2F::Divide if we agree that the binning will be always the same
- /* for (Int_t iVtx=0; iVtx<=fDataUncorrected->GetNbinsX(); iVtx++)
-  {
-    for (Int_t iEta=0; iEta<=fDataUncorrected->GetNbinsY(); iEta++)
-    {
-      Float_t correctionValue = 1;
-      if (correction)
-        correctionValue = correction->GetTrack2ParticleCorrection(fDataUncorrected->GetXaxis()->GetBinCenter(iVtx), fDataUncorrected->GetYaxis()->GetBinCenter(iEta), 1.0);
-
-      Float_t value = fDataUncorrected->GetBinContent(iVtx, iEta);
-      Float_t error = fDataUncorrected->GetBinError(iVtx, iEta);
-
-      Float_t correctedValue = value * correctionValue;
-      Float_t correctedError = error * correctionValue;
+  // In fData we have the track2particle and vertex reconstruction efficiency correction already applied
 
-      if (correctedValue != 0)
-      {
-        fData->SetBinContent(iVtx, iEta, correctedValue);
-        fData->SetBinError(iVtx, iEta, correctedError);
-      }
-    }
+  // integrate over pt (with pt cut)
+  fData->GetZaxis()->SetRange(fData->GetZaxis()->FindBin(ptCut), fData->GetZaxis()->GetNbins());
+  TH2D* vtxVsEta = dynamic_cast<TH2D*> (fData->Project3D("yx2"));
+  if (vtxVsEta == 0)
+  {
+    printf("ERROR: pt integration failed\n");
+    return;
   }
 
-  for (Int_t iEta=0; iEta<=fData->GetNbinsY(); iEta++)
+  new TCanvas;
+  vtxVsEta->Draw("COLZ");
+
+  for (Int_t iEta=0; iEta<=vtxVsEta->GetNbinsY(); iEta++)
   {
     // do we have several histograms for different vertex positions?
     Int_t vertexBinWidth = fVtx->GetNbinsX() / (kVertexBinning-1);
@@ -192,15 +185,23 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction)
       Float_t sumError2 = 0;
       for (Int_t iVtx = vertexBinBegin; iVtx < vertexBinEnd; iVtx++)
       {
-        if (fData->GetBinContent(iVtx, iEta) != 0)
+        if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
         {
-          sum = sum + fData->GetBinContent(iVtx, iEta);
-          sumError2 = sumError2 + TMath::Power(fData->GetBinError(iVtx, iEta),2);
+          sum = sum + vtxVsEta->GetBinContent(iVtx, iEta);
+          sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
         }
       }
 
-      Float_t dndeta = sum / totalEvents;
-      Float_t error  = TMath::Sqrt(sumError2) / totalEvents;
+      Float_t ptCutOffCorrection = correction->GetMeasuredFraction(ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta));
+      //ptCutOffCorrection = 1;
+      if (ptCutOffCorrection <= 0)
+      {
+        printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd);
+        continue;
+      }
+
+      Float_t dndeta = sum / totalEvents / ptCutOffCorrection;
+      Float_t error  = TMath::Sqrt(sumError2) / totalEvents / ptCutOffCorrection;
 
       dndeta = dndeta/fdNdEta[vertexPos]->GetBinWidth(iEta);
       error  = error/fdNdEta[vertexPos]->GetBinWidth(iEta);
@@ -208,7 +209,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction)
       fdNdEta[vertexPos]->SetBinContent(iEta, dndeta);
       fdNdEta[vertexPos]->SetBinError(iEta, error);
     }
-  }*/
+  }
 }
 
 //____________________________________________________________________
@@ -220,7 +221,10 @@ void dNdEtaAnalysis::SaveHistograms()
   gDirectory->cd(GetName());
 
   fData  ->Write();
+  AliPWG0Helper::CreateProjections(fData);
   fDataUncorrected->Write();
+  AliPWG0Helper::CreateProjections(fDataUncorrected);
+
   fVtx       ->Write();
   for (Int_t i=0; i<kVertexBinning; ++i)
     fdNdEta[i]    ->Write();
index 9db2a91..1d76b6c 100644 (file)
@@ -40,7 +40,7 @@ public:
   void FillTrack(Float_t vtx, Float_t eta, Float_t pt, Float_t weight);
   void FillEvent(Float_t vtx, Float_t weight);
 
-  void Finish(AlidNdEtaCorrection* correction);
+  void Finish(AlidNdEtaCorrection* correction, Float_t ptCut);
 
   void DrawHistograms();
   void LoadHistograms();
index 376d6e1..fb0f467 100644 (file)
@@ -58,7 +58,7 @@ void testAnalysis2(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kF
       }
     }
 
-    chain->SetProof(proof);
+    //chain->SetProof(proof);
   }
 
   // ########################################################
@@ -74,6 +74,17 @@ void testAnalysis2(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kF
   if (proof)
     proof->AddInput(esdTrackCuts);
 
+  if (aMC == kFALSE)
+  {
+    AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection();
+    dNdEtaCorrection->LoadHistograms("correction_map.root","dndeta_correction");
+    //dNdEtaCorrection->RemoveEdges(2, 0, 2);
+
+    chain->GetUserInfo()->Add(dNdEtaCorrection);
+    if (proof)
+      proof->AddInput(dNdEtaCorrection);
+  }
+
   TString selectorName = ((aMC == kFALSE) ? "AlidNdEtaAnalysisESDSelector" : "AlidNdEtaAnalysisMCSelector");
   AliLog::SetClassDebugLevel(selectorName, AliLog::kInfo);
 
@@ -90,7 +101,13 @@ void testAnalysis2(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kF
   TStopwatch timer;
   timer.Start();
 
-  Long64_t result = chain->Process(selectorName);
+  Long64_t result = -1;
+
+  if (proof != kFALSE)
+    result = chain->MakeTDSet()->Process(selectorName);
+  else
+    result = chain->Process(selectorName);
+
   if (result != 0)
   {
     printf("ERROR: Executing process failed with %d.\n", result);