PWG0/dNdEta/dNdEtaAnalysis.cxx

   1 /* $Id$ */
   2
   3 #include "dNdEtaAnalysis.h"
   4
   5 #include <TFile.h>
   6 #include <TH2F.h>
   7 #include <TH1D.h>
   8 #include <TMath.h>
   9 #include <TCanvas.h>
  10 #include <TCollection.h>
  11 #include <TIterator.h>
  12 #include <TList.h>
  13 #include <TLegend.h>
  14
  15 #include "dNdEtaCorrection.h"
  16
  17 //____________________________________________________________________
  18 ClassImp(dNdEtaAnalysis)
  19
  20 //____________________________________________________________________
  21 dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title) :
  22 TNamed(name, title)
  23 {
  24   hEtaVsVtx  = new TH2F(Form("%s_eta_vs_vtx", name),"",80,-20,20,120,-6,6);
  25   hEtaVsVtx->SetXTitle("vtx z [cm]");
  26   hEtaVsVtx->SetYTitle("#eta");
  27
  28   hEtaVsVtxUncorrected = dynamic_cast<TH2F*> (hEtaVsVtx->Clone(Form("%s_eta_vs_vtx_uncorrected", name)));
  29   hEtaVsVtxCheck = dynamic_cast<TH2F*> (hEtaVsVtx->Clone(Form("%s_eta_vs_vtx_check", name)));
  30   hVtx       = hEtaVsVtx->ProjectionX(Form("%s_vtx", name));
  31   for (Int_t i=0; i<kVertexBinning; ++i)
  32   {
  33     hdNdEta[i]    = hEtaVsVtx->ProjectionY(Form("%s_dNdEta_%d", name, i));
  34     hdNdEta[i]->SetYTitle("dN/d#eta");
  35   }
  36
  37   hEtaVsVtx->Sumw2();
  38   hVtx->Sumw2();
  39 }
  40
  41 //____________________________________________________________________
  42 void
  43 dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t c) {
  44   hEtaVsVtxUncorrected->Fill(vtx,eta);
  45   hEtaVsVtxCheck->Fill(vtx, eta, c);
  46 }
  47
  48 //____________________________________________________________________
  49 void
  50 dNdEtaAnalysis::FillEvent(Float_t vtx) {
  51   hVtx->Fill(vtx);
  52 }
  53
  54 //____________________________________________________________________
  55 void dNdEtaAnalysis::Finish(dNdEtaCorrection* correction)
  56 {
  57   // correct with correction values if available
  58   // TODO what do we do with the error?
  59   if (!correction)
  60     printf("INFO: No correction applied\n");
  61
  62   // this can be replaced by TH2F::Divide if we agree that the binning will be always the same
  63   for (Int_t iVtx=0; iVtx<=hEtaVsVtxUncorrected->GetNbinsX(); iVtx++)
  64   {
  65     for (Int_t iEta=0; iEta<=hEtaVsVtxUncorrected->GetNbinsY(); iEta++)
  66     {
  67       Float_t correctionValue = 1;
  68       if (correction)
  69         correctionValue = correction->GetCorrection(hEtaVsVtxUncorrected->GetXaxis()->GetBinCenter(iVtx), hEtaVsVtxUncorrected->GetYaxis()->GetBinCenter(iEta));
  70
  71       Float_t value = hEtaVsVtxUncorrected->GetBinContent(iVtx, iEta);
  72       Float_t error = hEtaVsVtxUncorrected->GetBinError(iVtx, iEta);
  73
  74       Float_t correctedValue = value * correctionValue;
  75       Float_t correctedError = error * correctionValue;
  76
  77       if (correctedValue != 0)
  78       {
  79         hEtaVsVtx->SetBinContent(iVtx, iEta, correctedValue);
  80         hEtaVsVtx->SetBinError(iVtx, iEta, correctedError);
  81       }
  82     }
  83   }
  84
  85   // normalize with n events (per vtx)
  86   for (Int_t iVtx=0; iVtx<=hVtx->GetNbinsX(); iVtx++) {
  87     Float_t nEvents      = hVtx->GetBinContent(iVtx);
  88     Float_t nEventsError = hVtx->GetBinError(iVtx);
  89
  90     if (nEvents==0) continue;
  91
  92     for (Int_t iEta=0; iEta<=hEtaVsVtx->GetNbinsY(); iEta++) {
  93       Float_t value = hEtaVsVtx->GetBinContent(iVtx, iEta) / nEvents;
  94       if (value==0) continue;
  95       Float_t error = hEtaVsVtx->GetBinError(iVtx, iEta)/nEvents;
  96       error = TMath::Sqrt(TMath::Power(hEtaVsVtx->GetBinError(iVtx, iEta)/
  97                                        hEtaVsVtx->GetBinContent(iVtx, iEta),2) +
  98                           TMath::Power(nEventsError/nEvents,2));
  99       hEtaVsVtx->SetBinContent(iVtx, iEta, value);
 100       hEtaVsVtx->SetBinError(iVtx,   iEta, error);
 101     }
 102
 103     //debug
 104     for (Int_t iEta=0; iEta<=hEtaVsVtxCheck->GetNbinsY(); iEta++) {
 105       Float_t value = hEtaVsVtxCheck->GetBinContent(iVtx, iEta) / nEvents;
 106       if (value==0) continue;
 107       Float_t error = hEtaVsVtxCheck->GetBinError(iVtx, iEta)/nEvents;
 108       error = TMath::Sqrt(TMath::Power(hEtaVsVtxCheck->GetBinError(iVtx, iEta)/
 109                                        hEtaVsVtxCheck->GetBinContent(iVtx, iEta),2) +
 110                           TMath::Power(nEventsError/nEvents,2));
 111       hEtaVsVtxCheck->SetBinContent(iVtx, iEta, value);
 112       hEtaVsVtxCheck->SetBinError(iVtx,   iEta, error);
 113     }
 114   }
 115
 116   // then take the wieghted average for each eta
 117   // is this the right way to do it???
 118   for (Int_t iEta=0; iEta<=hEtaVsVtx->GetNbinsY(); iEta++) {
 119     Float_t sum           = 0;
 120     Float_t sumError2     = 0;
 121     Int_t   nMeasurements = 0;
 122
 123     Float_t sumXw = 0;
 124     Float_t sumW  = 0;
 125
 126     // do we have several histograms for different vertex positions?
 127     Int_t vertexBinWidth = hVtx->GetNbinsX() / (kVertexBinning-1);
 128     for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
 129     {
 130       Int_t vertexBinBegin = 0;
 131       Int_t vertexBinEnd = hVtx->GetNbinsX();
 132
 133       // the first histogram is always for the whole vertex range
 134       if (vertexPos > 0)
 135       {
 136         vertexBinBegin = vertexBinWidth * (vertexPos-1);
 137         vertexBinEnd = vertexBinBegin + vertexBinWidth;
 138       }
 139
 140       for (Int_t iVtx=vertexBinBegin; iVtx<=vertexBinEnd; iVtx++) {
 141         if (hVtx->GetBinContent(iVtx)==0)             continue;
 142         if (hEtaVsVtxCheck->GetBinContent(iVtx, iEta)==0) continue;
 143
 144         Float_t w = 1/TMath::Power(hEtaVsVtx->GetBinError(iVtx, iEta),2);
 145         sumXw = sumXw + hEtaVsVtxCheck->GetBinContent(iVtx, iEta)*w;
 146         sumW  = sumW + w;
 147
 148         sum = sum + hEtaVsVtxCheck->GetBinContent(iVtx, iEta);
 149         sumError2 = sumError2 + TMath::Power(hEtaVsVtxCheck->GetBinError(iVtx, iEta),2);
 150         nMeasurements++;
 151       }
 152       Float_t dndeta = 0;
 153       Float_t error  = 0;
 154
 155       if (nMeasurements!=0) {
 156         dndeta = sum/Float_t(nMeasurements);
 157         error  = TMath::Sqrt(sumError2)/Float_t(nMeasurements);
 158
 159         dndeta = sumXw/sumW;
 160         error  = 1/TMath::Sqrt(sumW);
 161
 162         dndeta = dndeta/hdNdEta[vertexPos]->GetBinWidth(iEta);
 163         error  = error/hdNdEta[vertexPos]->GetBinWidth(iEta);
 164
 165         hdNdEta[vertexPos]->SetBinContent(iEta, dndeta);
 166         hdNdEta[vertexPos]->SetBinError(iEta, error);
 167       }
 168     }
 169   }
 170 }
 171
 172 //____________________________________________________________________
 173 void
 174 dNdEtaAnalysis::SaveHistograms() {
 175
 176   gDirectory->mkdir(GetName());
 177   gDirectory->cd(GetName());
 178
 179   hEtaVsVtx  ->Write();
 180   hEtaVsVtxUncorrected->Write();
 181   hVtx       ->Write();
 182   for (Int_t i=0; i<kVertexBinning; ++i)
 183     hdNdEta[i]    ->Write();
 184
 185   gDirectory->cd("../");
 186 }
 187
 188 //____________________________________________________________________
 189 void dNdEtaAnalysis::DrawHistograms()
 190 {
 191   TCanvas* canvas = new TCanvas("dNdEtaAnalysis", "dNdEtaAnalysis", 1200, 800);
 192   canvas->Divide(3, 2);
 193
 194   canvas->cd(1);
 195   if (hEtaVsVtx)
 196     hEtaVsVtx->Draw("COLZ");
 197
 198   canvas->cd(2);
 199   if (hEtaVsVtxCheck)
 200     hEtaVsVtxCheck->Draw("COLZ");
 201
 202   canvas->cd(3);
 203   if (hEtaVsVtxUncorrected)
 204     hEtaVsVtxUncorrected->Draw("COLZ");
 205
 206   canvas->cd(4);
 207   TH2F* clone = (TH2F*) hEtaVsVtxCheck->Clone("clone");
 208   clone->Divide(hEtaVsVtx);
 209   clone->Draw("COLZ");
 210
 211   canvas->cd(5);
 212   if (hVtx)
 213     hVtx->Draw();
 214
 215   canvas->cd(6);
 216   if (hdNdEta[0])
 217     hdNdEta[0]->Draw();
 218
 219     // histograms for different vertices?
 220   if (kVertexBinning > 0)
 221   {
 222       // doesnt work, but i dont get it, giving up...
 223     /*TCanvas* canvas2 =*/ new TCanvas("dNdEtaAnalysisVtx", "dNdEtaAnalysisVtx", 450, 450);
 224     //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
 225     //printf("%d\n", yPads);
 226     //canvas2->Divide(2, yPads);
 227
 228     TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9);
 229
 230     for (Int_t i=1; i<kVertexBinning; ++i)
 231     {
 232       //canvas2->cd(i-1);
 233       //printf("%d\n", i);
 234       if (hdNdEta[i])
 235       {
 236         hdNdEta[i]->SetLineColor(i);
 237         hdNdEta[i]->Draw((i == 1) ? "" : "SAME");
 238         legend->AddEntry(hdNdEta[i], Form("Vtx Bin %d", i-1));
 239       }
 240     }
 241
 242     legend->Draw();
 243   }
 244 }
 245
 246 Long64_t dNdEtaAnalysis::Merge(TCollection* list)
 247 {
 248   // Merges a list of dNdEtaAnalysis objects with this one.
 249   // This is needed for PROOF.
 250   // Returns the number of merged objects (including this)
 251
 252   if (!list)
 253     return 0;
 254
 255   if (list->IsEmpty())
 256     return 1;
 257
 258   TIterator* iter = list->MakeIterator();
 259   TObject* obj;
 260
 261   // sub collections
 262   const Int_t nCollections = kVertexBinning + 3;
 263   TList* collections[nCollections];
 264   for (Int_t i=0; i<nCollections; ++i)
 265     collections[i] = new TList;
 266
 267   Int_t count = 0;
 268   while ((obj = iter->Next()))
 269   {
 270     dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj);
 271     if (entry == 0)
 272       continue;
 273
 274     collections[0]->Add(entry->hEtaVsVtx);
 275     collections[1]->Add(entry->hEtaVsVtxUncorrected);
 276     collections[2]->Add(entry->hVtx);
 277
 278     for (Int_t i=0; i<kVertexBinning; ++i)
 279       collections[3+i]->Add(entry->hdNdEta[i]);
 280
 281     ++count;
 282   }
 283
 284   hEtaVsVtx->Merge(collections[0]);
 285   hEtaVsVtxUncorrected->Merge(collections[1]);
 286   hVtx->Merge(collections[2]);
 287   for (Int_t i=0; i<kVertexBinning; ++i)
 288     hdNdEta[i]->Merge(collections[3+i]);
 289
 290   for (Int_t i=0; i<nCollections; ++i)
 291     delete collections[i];
 292
 293   return count+1;
 294 }