CORRFW/test/testCFContainers.C

   1 #include <Riostream.h>
   2
   3 extern TRandom *gRandom;
   4 extern TBenchmark *gBenchmark;
   5 extern TSystem *gSystem;
   6
   7 void testCFContainers(){
   8
   9   // simple example macros for usage of a N-dim container (AliCFContainer)
  10   // handling a set of grids to accumulate data at different
  11   // selection steps, & derive efficiency (this is stored in AliCFEffGrid)
  12   // book, fill and draw some histos
  13   // The efficiency is then used to correct the data (trivially self-correct,
  14   // in this example)
  15
  16   gROOT->SetStyle("Plain");
  17   gStyle->SetPalette(1);
  18   gStyle->SetOptStat(111110);
  19   gStyle->SetPalette(1);
  20   gStyle->SetCanvasColor(0);
  21   gStyle->SetFrameFillColor(0);
  22
  23   gSystem->SetIncludePath("-I. -I$ALICE_ROOT/include  -I$ROOTSYS/include");
  24   gSystem->Load("libANALYSIS.so");
  25   gSystem->Load("$ALICE_ROOT/CORRFW/libCORRFW.so") ;
  26
  27   //Setting up the container grid...
  28
  29   const Int_t nstep=2; //number of selection steps  (just 2 in this ex)
  30
  31   const Int_t nvar=4; //number of variables on the grid:pt,vtx
  32
  33   const Int_t nbin1=6; //bins in pt
  34   const Int_t nbin2=10; //bins in eta
  35   const Int_t nbin3=18; //bins in phi
  36   const Int_t nbin4=10; //bins in vertex
  37
  38
  39   //Flag the sel steps. In this example, we have two, may be any nstep
  40   Int_t stepGen=0;
  41   Int_t stepRec=1;
  42
  43   //the sensitive variables, their indeces
  44   Int_t ipt =0;
  45   Int_t ieta=1;
  46   Int_t iphi=2;
  47   Int_t ivtx=3;
  48
  49   //arrays for the number of bins in each dimension
  50   const Int_t iBin[nvar] ={nbin1,nbin2,nbin3,nbin4};
  51
  52   //arrays for bin limits
  53   Double_t binLim1[nbin1+1];
  54   Double_t binLim2[nbin2+1];
  55   Double_t binLim3[nbin3+1];
  56   Double_t binLim4[nbin4+1];
  57
  58   for(Int_t i=0;i<=nbin1;i++){
  59     // pt [0,3] GeV/c
  60     binLim1[i]=i*0.5;
  61   }
  62
  63   for(Int_t i=0;i<=nbin2;i++){
  64     // eta [-1,1]
  65     binLim2[i]=i*0.2-1.;
  66   }
  67   for(Int_t i=0;i<=nbin3;i++){
  68     //phi [0,360]
  69     binLim3[i]=i*20.;
  70   }
  71   for(Int_t i=0;i<=nbin4;i++){
  72     //vertex [-20,20] cm
  73     binLim4[i]=i*4.-20.;
  74   }
  75
  76
  77   //the nstep grids "container"
  78   AliCFContainer *cont = new AliCFContainer("cont","example of  container",nstep,nvar,iBin);
  79
  80   //setting the bin limits
  81   cont->SetBinLimits(ipt,binLim1);
  82   cont->SetBinLimits(ieta,binLim2);
  83   cont->SetBinLimits(iphi,binLim3);
  84   cont->SetBinLimits(ivtx,binLim4);
  85
  86   //Start filling the mc and the data
  87
  88   //data sample (1M tracks)
  89   Int_t nev=1000000;
  90   Int_t seed =1234;
  91   gRandom->SetSeed(seed);
  92   Double_t Value[nvar];
  93   for(Int_t iev =0;iev<nev;iev++){
  94     Float_t y=gRandom->Rndm();
  95     Float_t pt=-TMath::Log(y)/0.5; //pt, exponential
  96     Double_t eta=2.*gRandom->Rndm()-1.;//flat in eta
  97     Double_t phi=360.*gRandom->Rndm(); //flat in phi
  98     Float_t vtx=gRandom->Gaus( 0,5.);//gaussian in vertex
  99     Value[ipt]=pt;
 100     Value[ieta]=eta;
 101     Value[iphi]=phi;
 102     Value[ivtx]=vtx;
 103     cont->Fill(Value, stepGen); //fill the efficiency denominator, sel step=0
 104     Float_t rndm=gRandom->Rndm();
 105     //simulate 80% constant efficiency everywhere
 106     if(rndm<0.8){
 107       cont->Fill(Value,stepRec); //fill the efficiency denominator, sel step =1
 108     }
 109   }
 110
 111 // Save it to a file
 112    cont->Save("container.root");
 113   //delete it
 114    delete cont;
 115
 116 // Read the  container from file
 117    TFile *file = new TFile("container.root");
 118    AliCFContainer *data = (AliCFContainer*) (file->Get("cont"));
 119
 120   // Make some 1 & 2-D projections..
 121   // pt and vertex, generator and reconstructed level
 122   TCanvas *cmc =new TCanvas("cmc","The  distributions",0,300,900,900);
 123   cmc->Divide(2,2);
 124   cmc->cd(1);
 125   TH1D *hpt1a = data->ShowProjection(ipt, stepGen);
 126   hpt1a->SetMinimum(0.01);
 127   hpt1a->Draw();
 128   cmc->cd(2);
 129   TH1D *hpt1b = data->ShowProjection(ipt, stepRec);
 130   hpt1b->SetMinimum(0.01);
 131   hpt1b->Draw();
 132   cmc->cd(3);
 133   TH2D *hptvtx1a = data->ShowProjection(ipt,ivtx, stepGen);
 134   hptvtx1a->SetMinimum(0.01);
 135   hptvtx1a->Draw("lego");
 136   cmc->cd(4);
 137   TH2D *hptvtx1b = data->ShowProjection(ipt,ivtx, stepRec);
 138   hptvtx1b->SetMinimum(0.01);
 139   hptvtx1b->Draw("lego");
 140   cmc->Print("data.gif");
 141
 142
 143   //construct the efficiency grid from the data container
 144   AliCFEffGrid *eff = new AliCFEffGrid("eff"," The efficiency",*data);
 145   eff->CalculateEfficiency(stepRec,stepGen); //eff= step1/step0
 146
 147   //The efficiency along pt and vertex, and 2-D projection
 148   TCanvas *ceff =new TCanvas("ceff"," Efficiency",0,300,900,300);
 149   ceff->Divide(3,1);
 150   ceff->cd(1);
 151   TH1D *hpt2a = eff->Project(ipt); //the efficiency vs pt
 152   hpt2a->SetMinimum(0.01);
 153   hpt2a->Draw();
 154   ceff->cd(2);
 155   TH1D *hvtx2a = eff->Project(ivtx); //the efficiency vs vtx
 156   hvtx2a->SetMinimum(0.01);
 157   hvtx2a->Draw();
 158   ceff->cd(3);
 159   TH2D *hptvtx2a = eff->Project(ipt,ivtx); //look at the numerator
 160   hptvtx2a->SetMinimum(0.01);
 161   hptvtx2a->SetMinimum(0.01);
 162   hptvtx2a->Draw("lego");
 163
 164   ceff->Print("eff.gif");
 165
 166   //get the corrected data grid
 167   AliCFDataGrid *corrdata = new AliCFDataGrid("corrdata","corrected data",*data);
 168   //correct selection step "reconstructed"
 169   corrdata->SetMeasured(stepRec); //set data to be corrected
 170   corrdata->ApplyEffCorrection(*eff);//apply the correction for efficiency
 171
 172   //The observed data, the corrected ones and the "MC truth" distributions
 173   //vs pt and vtx
 174   TCanvas *ccorrdata =new TCanvas("ccorrdata"," corrected data",0,300,900,900);
 175   ccorrdata->Divide(2,2);
 176   ccorrdata->cd(1);
 177   TH1D *hpt3a = corrdata->GetData()->Project(ipt); //uncorrected data
 178   hpt3a->SetMinimum(0.01);
 179   hpt3a->Draw();
 180   ccorrdata->cd(2);
 181   TH1D *hpt3b = corrdata->Project(ipt); //corrected data
 182   hpt3b->SetMinimum(0.01);
 183   hpt3b->Draw();
 184   ccorrdata->cd(3);
 185   TH1D *hvtx3a = corrdata->GetData()->Project(ivtx); //uncorrected data
 186   hvtx3a->SetMinimum(0.01);
 187   hvtx3a->Draw();
 188   ccorrdata->cd(4);
 189   TH1D *hvtx3b = corrdata->Project(ivtx); //corrected data
 190   hvtx3b->SetMinimum(0.01);
 191   hvtx3b->Draw();
 192   ccorrdata->Print("corrdata.gif");
 193
 194 }