Implemented ZDC time cut in phsyics selection and in trigger analysis for MC. Trackle...
[u/mrichter/AliRoot.git] / 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->Load("libANALYSIS.so");
24   gSystem->Load("libANALYSISalice.so");
25   gSystem->Load("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=100000;
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 }