.so cleanup: removed from gSystem->Load()
[u/mrichter/AliRoot.git] / CORRFW / test / testCFContainers.C
CommitLineData
310aed5c 1#include <Riostream.h>
2
3extern TRandom *gRandom;
4extern TBenchmark *gBenchmark;
5extern TSystem *gSystem;
6
7void 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
4070f709 23 gSystem->Load("libANALYSIS");
24 gSystem->Load("libANALYSISalice");
25 gSystem->Load("libCORRFW") ;
310aed5c 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"
82711910 78 AliCFContainer *cont = new AliCFContainer("cont","example of container",nstep,nvar,iBin);
310aed5c 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)
fb494025 89 Int_t nev=100000;
310aed5c 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
fb494025 111 // Save it to a file
310aed5c 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
fb494025 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");
310aed5c 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
db6722a5 164 ceff->Print("eff.gif");
310aed5c 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();
310aed5c 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();
310aed5c 192 ccorrdata->Print("corrdata.gif");
193
194}