]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG0/dNdEta/dNdEtaAnalysis.cxx
o) proper splitting of STEER dependent classes (PWG0dep) and ESD dependent classes
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / dNdEtaAnalysis.cxx
index ed2e1c184c34f8cdd43d96584cf64c4791c20a9a..4e21461b867b53204f493f368d660c79ab4af755 100644 (file)
@@ -19,168 +19,218 @@ ClassImp(dNdEtaAnalysis)
 
 //____________________________________________________________________
 dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title) :
-TNamed(name, title)
+  TNamed(name, title),
+  fEtaVsVtx(0),
+  fEtaVsVtxUncorrected(0),
+  fVtx(0)
 {
-  hEtaVsVtx  = new TH2F(Form("%s_eta_vs_vtx", name),"",80,-20,20,120,-6,6);
-  hEtaVsVtx->SetXTitle("vtx z [cm]");
-  hEtaVsVtx->SetYTitle("#eta");
+  // constructor
 
-  hEtaVsVtxUncorrected = dynamic_cast<TH2F*> (hEtaVsVtx->Clone(Form("%s_eta_vs_vtx_uncorrected", name)));
-  hEtaVsVtxCheck = dynamic_cast<TH2F*> (hEtaVsVtx->Clone(Form("%s_eta_vs_vtx_check", name)));
-  hVtx       = hEtaVsVtx->ProjectionX(Form("%s_vtx", name));
+  fEtaVsVtx  = new TH2F(Form("%s_eta_vs_vtx", name),"",80,-20,20,120,-6,6);
+  fEtaVsVtx->SetXTitle("vtx z [cm]");
+  fEtaVsVtx->SetYTitle("#eta");
+
+  fEtaVsVtxUncorrected = dynamic_cast<TH2F*> (fEtaVsVtx->Clone(Form("%s_eta_vs_vtx_uncorrected", name)));
+  fVtx       = fEtaVsVtx->ProjectionX(Form("%s_vtx", name));
   for (Int_t i=0; i<kVertexBinning; ++i)
   {
-    hdNdEta[i]    = hEtaVsVtx->ProjectionY(Form("%s_dNdEta_%d", name, i));
-    hdNdEta[i]->SetYTitle("dN/d#eta");
+    fdNdEta[i]    = fEtaVsVtx->ProjectionY(Form("%s_dNdEta_%d", name, i));
+    fdNdEta[i]->SetYTitle("dN/d#eta");
   }
 
-  hEtaVsVtx->Sumw2();
-  hVtx->Sumw2();
+  fEtaVsVtx->Sumw2();
+  fVtx->Sumw2();
 }
 
 //____________________________________________________________________
-void
-dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t c) {
-  hEtaVsVtxUncorrected->Fill(vtx,eta);
-  hEtaVsVtxCheck->Fill(vtx, eta, c);
+dNdEtaAnalysis::~dNdEtaAnalysis()
+{
+  // destructor
+
+  delete fEtaVsVtx;
+  fEtaVsVtx = 0;
+
+  delete fEtaVsVtxUncorrected;
+  fEtaVsVtxUncorrected = 0;
+
+  delete fVtx;
+  fVtx = 0;
+
+  for (Int_t i=0; i<kVertexBinning; ++i)
+  {
+    delete fdNdEta[i];
+    fdNdEta[i] = 0;
+  }
+}
+
+//_____________________________________________________________________________
+dNdEtaAnalysis::dNdEtaAnalysis(const dNdEtaAnalysis &c) :
+  TNamed(c),
+  fEtaVsVtx(0),
+  fEtaVsVtxUncorrected(0),
+  fVtx(0)
+{
+  //
+  // dNdEtaAnalysis copy constructor
+  //
+
+  ((dNdEtaAnalysis &) c).Copy(*this);
+}
+
+//_____________________________________________________________________________
+dNdEtaAnalysis &dNdEtaAnalysis::operator=(const dNdEtaAnalysis &c)
+{
+  //
+  // Assignment operator
+  //
+
+  if (this != &c) ((dNdEtaAnalysis &) c).Copy(*this);
+  return *this;
+}
+
+//_____________________________________________________________________________
+void dNdEtaAnalysis::Copy(TObject &c) const
+{
+  //
+  // Copy function
+  //
+
+  dNdEtaAnalysis& target = (dNdEtaAnalysis &) c;
+
+  target.fEtaVsVtx = dynamic_cast<TH2F*> (fEtaVsVtx->Clone());
+  target.fEtaVsVtxUncorrected = dynamic_cast<TH2F*> (fEtaVsVtxUncorrected->Clone());
+  target.fVtx = dynamic_cast<TH1D*> (fVtx->Clone());
+
+  for (Int_t i=0; i<kVertexBinning; ++i)
+    target.fdNdEta[i] = dynamic_cast<TH1D*> (fdNdEta[i]->Clone());
+
+  TNamed::Copy((TNamed &) c);
 }
 
 //____________________________________________________________________
-void
-dNdEtaAnalysis::FillEvent(Float_t vtx) {
-  hVtx->Fill(vtx);
+void dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta)
+{
+  // fills a track into the histograms
+
+  fEtaVsVtxUncorrected->Fill(vtx,eta);
+}
+
+//____________________________________________________________________
+void dNdEtaAnalysis::FillEvent(Float_t vtx)
+{
+  // fills an event into the histograms
+
+  fVtx->Fill(vtx);
 }
 
 //____________________________________________________________________
 void dNdEtaAnalysis::Finish(dNdEtaCorrection* correction)
 {
   // correct with correction values if available
+
   // TODO what do we do with the error?
   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<=hEtaVsVtxUncorrected->GetNbinsX(); iVtx++)
+  for (Int_t iVtx=0; iVtx<=fEtaVsVtxUncorrected->GetNbinsX(); iVtx++)
   {
-    for (Int_t iEta=0; iEta<=hEtaVsVtxUncorrected->GetNbinsY(); iEta++)
+    for (Int_t iEta=0; iEta<=fEtaVsVtxUncorrected->GetNbinsY(); iEta++)
     {
       Float_t correctionValue = 1;
       if (correction)
-        correctionValue = correction->GetCorrection(hEtaVsVtxUncorrected->GetXaxis()->GetBinCenter(iVtx), hEtaVsVtxUncorrected->GetYaxis()->GetBinCenter(iEta));
+        correctionValue = correction->GetCorrection(fEtaVsVtxUncorrected->GetXaxis()->GetBinCenter(iVtx), fEtaVsVtxUncorrected->GetYaxis()->GetBinCenter(iEta));
 
-      Float_t value = hEtaVsVtxUncorrected->GetBinContent(iVtx, iEta);
-      Float_t error = hEtaVsVtxUncorrected->GetBinError(iVtx, iEta);
+      Float_t value = fEtaVsVtxUncorrected->GetBinContent(iVtx, iEta);
+      Float_t error = fEtaVsVtxUncorrected->GetBinError(iVtx, iEta);
 
       Float_t correctedValue = value * correctionValue;
       Float_t correctedError = error * correctionValue;
 
       if (correctedValue != 0)
       {
-        hEtaVsVtx->SetBinContent(iVtx, iEta, correctedValue);
-        hEtaVsVtx->SetBinError(iVtx, iEta, correctedError);
+        fEtaVsVtx->SetBinContent(iVtx, iEta, correctedValue);
+        fEtaVsVtx->SetBinError(iVtx, iEta, correctedError);
       }
     }
   }
 
-  // normalize with n events (per vtx)
-  for (Int_t iVtx=0; iVtx<=hVtx->GetNbinsX(); iVtx++) {
-    Float_t nEvents      = hVtx->GetBinContent(iVtx);
-    Float_t nEventsError = hVtx->GetBinError(iVtx);
-
-    if (nEvents==0) continue;
-
-    for (Int_t iEta=0; iEta<=hEtaVsVtx->GetNbinsY(); iEta++) {
-      Float_t value = hEtaVsVtx->GetBinContent(iVtx, iEta) / nEvents;
-      if (value==0) continue;
-      Float_t error = hEtaVsVtx->GetBinError(iVtx, iEta)/nEvents;
-      error = TMath::Sqrt(TMath::Power(hEtaVsVtx->GetBinError(iVtx, iEta)/
-                                      hEtaVsVtx->GetBinContent(iVtx, iEta),2) +
-                         TMath::Power(nEventsError/nEvents,2));
-      hEtaVsVtx->SetBinContent(iVtx, iEta, value);
-      hEtaVsVtx->SetBinError(iVtx,   iEta, error);
-    }
-
-    //debug
-    for (Int_t iEta=0; iEta<=hEtaVsVtxCheck->GetNbinsY(); iEta++) {
-      Float_t value = hEtaVsVtxCheck->GetBinContent(iVtx, iEta) / nEvents;
-      if (value==0) continue;
-      Float_t error = hEtaVsVtxCheck->GetBinError(iVtx, iEta)/nEvents;
-      error = TMath::Sqrt(TMath::Power(hEtaVsVtxCheck->GetBinError(iVtx, iEta)/
-                                      hEtaVsVtxCheck->GetBinContent(iVtx, iEta),2) +
-                         TMath::Power(nEventsError/nEvents,2));
-      hEtaVsVtxCheck->SetBinContent(iVtx, iEta, value);
-      hEtaVsVtxCheck->SetBinError(iVtx,   iEta, error);
-    }
-  }
-
-  // then take the wieghted average for each eta
-  // is this the right way to do it???
-  for (Int_t iEta=0; iEta<=hEtaVsVtx->GetNbinsY(); iEta++) {
-    Float_t sum           = 0;
-    Float_t sumError2     = 0;
-    Int_t   nMeasurements = 0;    
-
-    Float_t sumXw = 0;
-    Float_t sumW  = 0;
-    
+  for (Int_t iEta=0; iEta<=fEtaVsVtx->GetNbinsY(); iEta++)
+  {
     // do we have several histograms for different vertex positions?
-    Int_t vertexBinWidth = hVtx->GetNbinsX() / (kVertexBinning-1);
+    Int_t vertexBinWidth = fVtx->GetNbinsX() / (kVertexBinning-1);
     for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
     {
-      Int_t vertexBinBegin = 0;
-      Int_t vertexBinEnd = hVtx->GetNbinsX();
+      Int_t vertexBinBegin = 1;
+      Int_t vertexBinEnd = fVtx->GetNbinsX() + 1;
 
       // the first histogram is always for the whole vertex range
       if (vertexPos > 0)
       {
-        vertexBinBegin = vertexBinWidth * (vertexPos-1);
+        vertexBinBegin = 1 + vertexBinWidth * (vertexPos-1);
         vertexBinEnd = vertexBinBegin + vertexBinWidth;
       }
 
-      for (Int_t iVtx=vertexBinBegin; iVtx<=vertexBinEnd; iVtx++) {
-        if (hVtx->GetBinContent(iVtx)==0)             continue;
-        if (hEtaVsVtx->GetBinContent(iVtx, iEta)==0) continue;
-
-        Float_t w = 1/TMath::Power(hEtaVsVtx->GetBinError(iVtx, iEta),2);
-        sumXw = sumXw + hEtaVsVtx->GetBinContent(iVtx, iEta)*w;
-        sumW  = sumW + w;
-
-        sum = sum + hEtaVsVtx->GetBinContent(iVtx, iEta);
-        sumError2 = sumError2 + TMath::Power(hEtaVsVtx->GetBinError(iVtx, iEta),2);      
-        nMeasurements++;
+      Float_t totalEvents = fVtx->Integral(vertexBinBegin, vertexBinEnd - 1);
+      if (totalEvents == 0)
+      {
+        printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
+        continue;
       }
-      Float_t dndeta = 0;
-      Float_t error  = 0;
 
-      if (nMeasurements!=0) {
-        dndeta = sum/Float_t(nMeasurements);
-        error  = TMath::Sqrt(sumError2)/Float_t(nMeasurements);
+      Float_t sum = 0;
+      Float_t sumError2 = 0;
+      for (Int_t iVtx = vertexBinBegin; iVtx < vertexBinEnd; iVtx++)
+      {
+        if (fEtaVsVtx->GetBinContent(iVtx, iEta) != 0)
+        {
+          sum = sum + fEtaVsVtx->GetBinContent(iVtx, iEta);
+          sumError2 = sumError2 + TMath::Power(fEtaVsVtx->GetBinError(iVtx, iEta),2);
+        }
+      }
 
-        dndeta = sumXw/sumW;
-        error  = 1/TMath::Sqrt(sumW);
+      Float_t dndeta = sum / totalEvents;
+      Float_t error  = TMath::Sqrt(sumError2) / totalEvents;
 
-        dndeta = dndeta/hdNdEta[vertexPos]->GetBinWidth(iEta);
-        error  = error/hdNdEta[vertexPos]->GetBinWidth(iEta);
+      dndeta = dndeta/fdNdEta[vertexPos]->GetBinWidth(iEta);
+      error  = error/fdNdEta[vertexPos]->GetBinWidth(iEta);
 
-        hdNdEta[vertexPos]->SetBinContent(iEta, dndeta);
-        hdNdEta[vertexPos]->SetBinError(iEta, error);
-      }
+      fdNdEta[vertexPos]->SetBinContent(iEta, dndeta);
+      fdNdEta[vertexPos]->SetBinError(iEta, error);
     }
   }
 }
 
 //____________________________________________________________________
-void
-dNdEtaAnalysis::SaveHistograms() {
+void dNdEtaAnalysis::SaveHistograms()
+{
+  // save the histograms to a directory with the name of this class (retrieved from TNamed)
 
   gDirectory->mkdir(GetName());
   gDirectory->cd(GetName());
-  
-  hEtaVsVtx  ->Write();
-  hEtaVsVtxUncorrected->Write();
-  hVtx       ->Write();
+
+  fEtaVsVtx  ->Write();
+  fEtaVsVtxUncorrected->Write();
+  fVtx       ->Write();
   for (Int_t i=0; i<kVertexBinning; ++i)
-    hdNdEta[i]    ->Write();
+    fdNdEta[i]    ->Write();
+
+  gDirectory->cd("../");
+}
+
+void dNdEtaAnalysis::LoadHistograms()
+{
+  // loads the histograms from a directory with the name of this class (retrieved from TNamed)
+
+  gDirectory->cd(GetName());
+
+  fEtaVsVtx = dynamic_cast<TH2F*> (gDirectory->Get(fEtaVsVtx->GetName()));
+  fEtaVsVtxUncorrected = dynamic_cast<TH2F*> (gDirectory->Get(fEtaVsVtxUncorrected->GetName()));
+
+  fVtx = dynamic_cast<TH1D*> (gDirectory->Get(fVtx->GetName()));
+
+  for (Int_t i=0; i<kVertexBinning; ++i)
+    fdNdEta[i] = dynamic_cast<TH1D*> (gDirectory->Get(fdNdEta[i]->GetName()));
 
   gDirectory->cd("../");
 }
@@ -188,28 +238,26 @@ dNdEtaAnalysis::SaveHistograms() {
 //____________________________________________________________________
 void dNdEtaAnalysis::DrawHistograms()
 {
+  // draws the histograms
+
   TCanvas* canvas = new TCanvas("dNdEtaAnalysis", "dNdEtaAnalysis", 800, 800);
   canvas->Divide(2, 2);
 
   canvas->cd(1);
-  if (hEtaVsVtx)
-    hEtaVsVtx->Draw("COLZ");
+  if (fEtaVsVtx)
+    fEtaVsVtx->Draw("COLZ");
 
   canvas->cd(2);
-  if (hEtaVsVtxCheck)
-    hEtaVsVtxCheck->Draw("COLZ");
+  if (fEtaVsVtxUncorrected)
+    fEtaVsVtxUncorrected->Draw("COLZ");
 
   canvas->cd(3);
-  if (hEtaVsVtxUncorrected)
-    hEtaVsVtxUncorrected->Draw("COLZ");
-
-  /*canvas->cd(3);
-  if (hVtx)
-    hVtx->Draw();*/
+  if (fVtx)
+    fVtx->Draw();
 
   canvas->cd(4);
-  if (hdNdEta[0])
-    hdNdEta[0]->Draw();
+  if (fdNdEta[0])
+    fdNdEta[0]->Draw();
 
     // histograms for different vertices?
   if (kVertexBinning > 0)
@@ -222,15 +270,15 @@ void dNdEtaAnalysis::DrawHistograms()
 
     TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9);
 
-    for (Int_t i=1; i<kVertexBinning; ++i)
+    for (Int_t i=0; i<kVertexBinning; ++i)
     {
       //canvas2->cd(i-1);
       //printf("%d\n", i);
-      if (hdNdEta[i])
+      if (fdNdEta[i])
       {
-        hdNdEta[i]->SetLineColor(i);
-        hdNdEta[i]->Draw((i == 1) ? "" : "SAME");
-        legend->AddEntry(hdNdEta[i], Form("Vtx Bin %d", i-1));
+        fdNdEta[i]->SetLineColor(i+1);
+        fdNdEta[i]->Draw((i == 0) ? "" : "SAME");
+        legend->AddEntry(fdNdEta[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
       }
     }
 
@@ -266,21 +314,21 @@ Long64_t dNdEtaAnalysis::Merge(TCollection* list)
     if (entry == 0)
       continue;
 
-    collections[0]->Add(entry->hEtaVsVtx);
-    collections[1]->Add(entry->hEtaVsVtxUncorrected);
-    collections[2]->Add(entry->hVtx);
+    collections[0]->Add(entry->fEtaVsVtx);
+    collections[1]->Add(entry->fEtaVsVtxUncorrected);
+    collections[2]->Add(entry->fVtx);
 
     for (Int_t i=0; i<kVertexBinning; ++i)
-      collections[3+i]->Add(entry->hdNdEta[i]);
+      collections[3+i]->Add(entry->fdNdEta[i]);
 
     ++count;
   }
 
-  hEtaVsVtx->Merge(collections[0]);
-  hEtaVsVtxUncorrected->Merge(collections[1]);
-  hVtx->Merge(collections[2]);
+  fEtaVsVtx->Merge(collections[0]);
+  fEtaVsVtxUncorrected->Merge(collections[1]);
+  fVtx->Merge(collections[2]);
   for (Int_t i=0; i<kVertexBinning; ++i)
-    hdNdEta[i]->Merge(collections[3+i]);
+    fdNdEta[i]->Merge(collections[3+i]);
 
   for (Int_t i=0; i<nCollections; ++i)
     delete collections[i];