PWG0/CorrectionMatrix2D.cxx

   1 // ------------------------------------------------------
   2 //
   3 // Class to handle 2d-corrections.
   4 //
   5 // ------------------------------------------------------
   6 //
   7
   8 /* $Id$ */
   9
  10 #include <TFile.h>
  11 #include <TCanvas.h>
  12
  13 #include <AliLog.h>
  14
  15 #include "CorrectionMatrix2D.h"
  16
  17 //____________________________________________________________________
  18 ClassImp(CorrectionMatrix2D)
  19
  20 //____________________________________________________________________
  21 CorrectionMatrix2D::CorrectionMatrix2D(const CorrectionMatrix2D& c)
  22   : TNamed(c)
  23 {
  24   // copy constructor
  25   ((CorrectionMatrix2D &)c).Copy(*this);
  26 }
  27
  28 //____________________________________________________________________
  29 CorrectionMatrix2D::CorrectionMatrix2D(Char_t* name, Char_t* title,
  30                                        Int_t nBinX, Float_t Xmin, Float_t Xmax,
  31                                        Int_t nBinY, Float_t Ymin, Float_t Ymax)
  32   : TNamed(name, title)
  33 {
  34   //
  35   // constructor
  36   //
  37
  38
  39   fhMeas  = new TH2F(Form("meas_%s",name), Form("meas_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
  40   fhGene  = new TH2F(Form("gene_%s",name), Form("gene_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
  41   fhCorr  = new TH2F(Form("corr_%s",name), Form("corr_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
  42
  43   fhMeas->Sumw2();
  44   fhGene->Sumw2();
  45   fhCorr->Sumw2();
  46 }
  47
  48 //____________________________________________________________________
  49 CorrectionMatrix2D::CorrectionMatrix2D(Char_t* name,Char_t* title,
  50                                        Int_t nBinX, Float_t *X, Int_t nBinY, Float_t *Y)
  51   : TNamed(name, title)
  52 {
  53   //
  54   // constructor
  55   //
  56
  57   fhMeas  = new TH2F(Form("meas_%s",name), Form("meas_%s",title),  nBinX, X, nBinY, Y);
  58   fhGene  = new TH2F(Form("gene_%s",name), Form("gene_%s",title),  nBinX, X, nBinY, Y);
  59   fhCorr  = new TH2F(Form("corr_%s",name), Form("corr_%s",title),  nBinX, X, nBinY, Y);
  60
  61   fhMeas->Sumw2();
  62   fhGene->Sumw2();
  63   fhCorr->Sumw2();
  64 }
  65
  66
  67 //____________________________________________________________________
  68 CorrectionMatrix2D::~CorrectionMatrix2D() {
  69   //
  70   // destructor
  71   //
  72   if (fhMeas)  delete fhMeas;
  73   if (fhGene)  delete fhGene;
  74   if (fhCorr)  delete fhCorr;
  75 }
  76
  77 //____________________________________________________________________
  78 CorrectionMatrix2D &CorrectionMatrix2D::operator=(const CorrectionMatrix2D &c)
  79 {
  80   // assigment operator
  81
  82   if (this != &c)
  83     ((CorrectionMatrix2D &) c).Copy(*this);
  84
  85   return *this;
  86 }
  87
  88 //____________________________________________________________________
  89 TH1F* CorrectionMatrix2D::Get1DCorrection(Char_t* opt) {
  90   //
  91   // integrate the correction over one variable
  92   //
  93
  94   TH1D* meas1D = 0;
  95   TH1D* gene1D = 0;
  96
  97   if (strcmp(opt,"x")==0) {
  98     meas1D = fhMeas->ProjectionX();
  99     gene1D = fhGene->ProjectionX();
 100   }
 101   if (strcmp(opt,"y")==0) {
 102     meas1D = fhMeas->ProjectionY();
 103     gene1D = fhGene->ProjectionY();
 104   }
 105   gene1D->Sumw2();
 106
 107   gene1D->SetName(Form("corr_1D_%s",fName.Data()));
 108   gene1D->SetTitle(Form("corr_1D_%s",fName.Data()));
 109
 110   gene1D->Divide(gene1D, meas1D, 1, 1, "B");
 111
 112   return (TH1F*)gene1D;
 113 }
 114
 115
 116 //____________________________________________________________________
 117 void
 118 CorrectionMatrix2D::Copy(TObject& c) const
 119 {
 120   // copy function
 121
 122   CorrectionMatrix2D& target = (CorrectionMatrix2D &) c;
 123
 124   target.fhMeas  = fhMeas;
 125   target.fhGene  = fhGene;
 126   target.fhCorr  = fhCorr;
 127 }
 128
 129
 130 //________________________________________________________________________
 131 void CorrectionMatrix2D::SetAxisTitles(Char_t* titleX, Char_t* titleY)
 132 {
 133   //
 134   // method for setting the axis titles of the histograms
 135   //
 136
 137   fhMeas ->SetXTitle(titleX);  fhMeas ->SetYTitle(titleY);
 138   fhGene ->SetXTitle(titleX);  fhGene ->SetYTitle(titleY);
 139   fhCorr ->SetXTitle(titleX);  fhCorr ->SetYTitle(titleY);
 140 }
 141
 142 //____________________________________________________________________
 143 Long64_t CorrectionMatrix2D::Merge(TCollection* list) {
 144   // Merge a list of CorrectionMatrix2D objects with this (needed for
 145   // PROOF).
 146   // Returns the number of merged objects (including this).
 147
 148   if (!list)
 149     return 0;
 150
 151   if (list->IsEmpty())
 152     return 1;
 153
 154   TIterator* iter = list->MakeIterator();
 155   TObject* obj;
 156
 157   // collections of measured and generated histograms
 158   TList* collectionMeas = new TList;
 159   TList* collectionGene = new TList;
 160
 161   Int_t count = 0;
 162   while ((obj = iter->Next())) {
 163
 164     CorrectionMatrix2D* entry = dynamic_cast<CorrectionMatrix2D*> (obj);
 165     if (entry == 0)
 166       continue;
 167
 168     collectionMeas->Add(entry->GetMeasuredHistogram());
 169     collectionGene->Add(entry->GetGeneratedHistogram());
 170
 171     count++;
 172   }
 173   fhMeas->Merge(collectionMeas);
 174   fhGene->Merge(collectionGene);
 175
 176   // is this really faster than just adding the histograms in the list???
 177   delete collectionMeas;
 178   delete collectionGene;
 179
 180
 181   return count+1;
 182 }
 183
 184
 185 //____________________________________________________________________
 186 void CorrectionMatrix2D::Divide() {
 187   //
 188   // divide the histograms to get the correction
 189   //
 190
 191   if (!fhMeas || !fhGene)  return;
 192
 193   fhCorr->Divide(fhGene, fhMeas, 1,1,"B");
 194
 195 }
 196
 197 //____________________________________________________________________
 198 void
 199 CorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge)
 200 {
 201   // remove edges of correction histogram by removing
 202   // - bins with content less than cut
 203   // - bins next to bins with zero bin content
 204
 205   Int_t nBinsX = fhCorr->GetNbinsX();
 206   Int_t nBinsY = fhCorr->GetNbinsY();
 207
 208   // set bin content to zero for bins with content smaller cut
 209   for (Int_t bx=0; bx<=nBinsX; bx++) {
 210     for (Int_t by=0; by<=nBinsY; by++) {
 211       if (fhCorr->GetBinContent(bx,by)>cut) {
 212           fhCorr->SetBinContent(bx,by,0);
 213           fhCorr->SetBinError(bx,by,0);
 214       }
 215     }
 216   }
 217
 218   // set bin content to zero for bins next to bins with zero
 219   TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
 220   tmp->Reset();
 221
 222   Bool_t done = kFALSE;
 223   Int_t nBinsXCount = 0;
 224   Int_t nBinsYCount = 0;
 225   while (!done) {
 226     if (nBinsXCount<nBinsXedge)
 227       for (Int_t bx=0; bx<=nBinsX; bx++) {
 228         for (Int_t by=0; by<=nBinsY; by++) {
 229           if ((fhCorr->GetBinContent(bx+1,by)==0)||
 230               (fhCorr->GetBinContent(bx-1,by)==0))
 231             tmp->SetBinContent(bx,by,1);
 232
 233         }
 234       }
 235     if (nBinsYCount<nBinsYedge)
 236       for (Int_t bx=0; bx<=nBinsX; bx++) {
 237         for (Int_t by=0; by<=nBinsY; by++) {
 238           if ((fhCorr->GetBinContent(bx,by+1)==0)||
 239               (fhCorr->GetBinContent(bx,by-1)==0))
 240             tmp->SetBinContent(bx,by,1);
 241         }
 242       }
 243     for (Int_t bx=0; bx<=nBinsX; bx++) {
 244       for (Int_t by=0; by<=nBinsY; by++) {
 245         if (tmp->GetBinContent(bx,by)==1) {
 246           fhCorr->SetBinContent(bx,by,0);
 247           fhCorr->SetBinError(bx,by,0);
 248         }
 249       }
 250     }
 251     nBinsXCount++;
 252     nBinsYCount++;
 253     if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
 254   }
 255   tmp->Delete();
 256
 257 }
 258
 259 //____________________________________________________________________
 260 Bool_t CorrectionMatrix2D::LoadHistograms(Char_t* fileName, Char_t* dir) {
 261   //
 262   // loads the histograms from a file
 263   //
 264
 265   TFile* fin = TFile::Open(fileName);
 266
 267   if(!fin) {
 268     //Info("LoadHistograms",Form(" %s file does not exist",fileName));
 269     return kFALSE;
 270   }
 271
 272   if(fhGene)  {delete fhGene;  fhGene=0;}
 273   if(fhCorr)  {delete fhCorr;  fhCorr=0;}
 274   if(fhMeas)  {delete fhMeas;  fhMeas=0;}
 275
 276   fhMeas  = (TH2F*)fin->Get(Form("%s/meas_%s", dir,GetName()));
 277   if(!fhMeas)  Info("LoadHistograms","No meas  hist available");
 278
 279   fhGene  = (TH2F*)fin->Get(Form("%s/gene_%s",dir, GetName()));
 280   if(!fhGene)  Info("LoadHistograms","No gene  hist available");
 281
 282   fhCorr  = (TH2F*)fin->Get(Form("%s/corr_%s",dir, GetName()));
 283   if(!fhCorr)
 284   {
 285     Info("LoadHistograms","No corr  hist available");
 286     return kFALSE;
 287   }
 288
 289   return kTRUE;
 290 }
 291
 292
 293 //____________________________________________________________________
 294 void
 295 CorrectionMatrix2D::SaveHistograms() {
 296   //
 297   // saves the histograms
 298   //
 299
 300   fhMeas ->Write();
 301   fhGene ->Write();
 302
 303   if (fhCorr)
 304     fhCorr->Write();
 305 }
 306
 307 //____________________________________________________________________
 308 void CorrectionMatrix2D::DrawHistograms()
 309 {
 310   //
 311   // draws all the four histograms on one TCanvas
 312   //
 313
 314   TCanvas* canvas = new TCanvas(Form("correction_%s",fName.Data()),
 315                                 Form("correction_%s",fName.Data()), 800, 800);
 316   canvas->Divide(2, 2);
 317
 318   canvas->cd(1);
 319   if (fhMeas)
 320     fhMeas->Draw("COLZ");
 321
 322   canvas->cd(2);
 323   if (fhGene)
 324     fhGene->Draw("COLZ");
 325
 326   canvas->cd(3);
 327   if (fhCorr)
 328     fhCorr->Draw("COLZ");
 329
 330   canvas->cd(4);
 331
 332   // add: draw here the stat. errors of the correction histogram
 333
 334 }
 335
 336