New implementation of the forward multiplicity analysis.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / older / DrawResVtx.C
1 #include "DrawRes1D.C"
2
3 /** 
4  * @defgroup pwg2_forward_analysis_scripts PWG2 Forward analysis - scripts
5  *
6  * @ingroup pwg2_forward_analysis
7  */
8 /** 
9  * Example macro to loop over the event-by-event 2D histogram of 
10  * @f[
11  *   \frac{d^{2}N_{ch}}{d\eta\,d\phi}
12  * @f]
13  * stored in an AOD.  
14  * 
15  * The class needs the files &lt;<i>base</i>&gt;<tt>_hists.root</tt> 
16  * containing the histograms generated by AliForwardMultiplicity and 
17  * the file &lt;<i>base</i>&gt;<tt>_aods.root</tt> containing the tree 
18  * with AliAODEvent objects where AliAODForwardMult objects have been 
19  * added to in the branch <tt>Forward</tt>
20  * 
21  * @ingroup pwg2_forward_analysis_scripts
22  */
23 class DrawResVtx : public DrawRes1D
24 {
25 public:
26   /** 
27    * Constructor 
28    * 
29    * @param special If true, add to the list of 'specials'
30    */
31   DrawResVtx()
32     : fByVtx1D(0)
33   {}
34   //__________________________________________________________________
35   /** 
36    * Open the files &lt;<i>base</i>&gt;<tt>_hists.root</tt> 
37    * containing the histograms generated by AliForwardMultiplicity and 
38    * the file &lt;<i>base</i>&gt;<tt>_aods.root</tt> containing the tree 
39    * with AliAODEvent objects.
40    * 
41    * @param base  Base name of files 
42    * @param vzMin Minimum collision vertex z position to use
43    * @param vzMax Maximum collision vertex z position to use
44    * @param rebin Rebinning factor 
45    * 
46    * @return true on success, false otherwise 
47    */
48   Bool_t Open(const char* base, 
49               Double_t    vzMin=-10, 
50               Double_t    vzMax=10)
51   {
52     // Clear previously created data objects 
53     fByVtx1D.Delete();
54     return DrawRes1D::Open(base, vzMin, vzMax);
55   }
56
57   //__________________________________________________________________
58   /** 
59    * Make a 1D histogram of @f$ dN_{ch}/d\eta@f$ with the given name and 
60    * title 
61    * 
62    * @param name   Name of histogram 
63    * @param title  Title of histogram 
64    * @param a1     Eta axis 
65    * 
66    * @return Newly allocated 1D histogram object 
67    */
68   TH1D* Make1D(const char* name, const char* title, const TAxis* a1)
69   {
70     TH1D* ret = new TH1D(name, title,
71                          a1->GetNbins(), a1->GetXmin(), a1->GetXmax());
72     ret->SetXTitle("#eta");
73     ret->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
74     ret->Sumw2();
75     ret->SetMarkerColor(kRed+1);
76     ret->SetLineColor(kRed+1);
77     ret->SetMarkerStyle(20);
78     ret->SetMarkerSize(1);
79     ret->SetStats(0);
80     ret->SetDirectory(0);
81
82     return ret;
83   }
84   //__________________________________________________________________
85   /** 
86    * Utility function to set-up histograms based on the input 
87    * @f$ dd^{2}N_{ch}/d\eta\,d\phi@f$ histogram.   This member function 
88    * is called on the first event so that we have the proper binning 
89    * 
90    * @param templ Input histogram
91    * 
92    * @return true on succcess
93    */
94   Bool_t FirstEvent(const TH2D& templ) 
95   { 
96     if (!DrawRes1D::FirstEvent(templ)) return kFALSE;
97
98     const TAxis* etaAxis = templ.GetXaxis();
99     
100     // Generate the per-vertex bin histograms.  These will be the sum 
101     // of each of the per-event input histograms for a given vertex range. 
102     for (Int_t i = fBinVzMin; i <= fBinVzMax; i++) { 
103       TH1* h1 = Make1D(Form("dndeta_vz%02d", i), 
104                        Form("#frac{1}{N}#frac{dN_{ch}}{d#eta}, vtxbin=%d", i),
105                        etaAxis);
106       fByVtx1D.AddAtAndExpand(h1, i);
107     }
108     return kTRUE;
109   }
110
111
112   //__________________________________________________________________
113   /** 
114    * Process the events 
115    * 
116    * 
117    * @return true on success, false otherwise 
118    */
119   void AddContrib(const TH2D& hist, Int_t vtxBin)
120   {
121     // Get 'by-vertex' histograms 
122     TH1* tmp1 = static_cast<TH1D*>(fByVtx1D.At(vtxBin));
123     if (!tmp1) { 
124       Warning("Process", "No histogram for vertex bin %d)",vtxBin);
125       return;
126     }
127     // Make a projection on the X axis of the input histogram 
128     TH1D*    proj  = hist.ProjectionX("_px", 0, -1, "e");
129     
130     // Add to per-vertex histogram 
131     tmp1->Add(proj); 
132     
133     delete proj;
134   } 
135   //__________________________________________________________________
136   /** 
137    * Called at the end of the job. 
138    *
139    * It plots the result of the analysis in tree canvases.   
140    * - One shows the per-vertex accumalted histograms and compares them 
141    *   to the per-vertex histograms made in the analysis. 
142    * - Another shows the final @f$ dN_{ch}/d\eta@f$ calculated in 
143    *   different ways and compare those to the  @f$ dN_{ch}/d\eta@f$ made 
144    *   during the analysis. 
145    * - The last canvas shows the @f$ dd^{2}N_{ch}/d\eta\,d\phi@f$ histogram. 
146    * 
147    * @return true on succes, false otherwise 
148    */
149   TH1D* GetResult() 
150   {
151     // Get number of bins 
152     Int_t nBin = fBinVzMax-fBinVzMin+1;
153
154     // Make first canvas 
155     TCanvas* cVtx = new TCanvas("cVtx", "By Vertex", 1000, 700);
156     cVtx->SetBorderSize(0);
157     cVtx->SetBorderMode(0);
158     cVtx->Divide((nBin+.5)/2, 2, 0.0005, 0.0005);
159
160     // Loop over vertex histograms 
161     TDirectory* savdir = gDirectory;
162     for (Int_t i = fBinVzMin; i <= fBinVzMax; i++) { 
163       TDirectory* vtxDir = savdir->mkdir(Form("vtx%02d", i));
164       vtxDir->cd();
165
166       TH1D* vh1 = static_cast<TH1D*>(fByVtx1D.At(i));
167       if (!vh1) { 
168         Error("Finish", "No histogram at %d", i);
169         continue;
170       }
171       
172       fTotal1D->Add(vh1);
173
174       // Scale and add to output 
175       if (fVtx->GetBinContent(i) > 0)
176         vh1->Scale(1. / fVtx->GetBinContent(i), "width");
177
178       // Write to output file 
179       vh1->Write(); 
180
181       // Divide the pad
182       TVirtualPad* pad = cVtx->cd(i-fBinVzMin+1);
183       pad->Divide(1,2,0.005,0);
184       
185       // Draw the per-vertex histogram 
186       pad->cd(1);
187       vh1->DrawClone();
188       
189       // Get the same histogram from the analysis and draw it 
190       TProfile* p = 
191         static_cast<TProfile*>(fCollect
192                                ->FindObject(Form("dndeta_vtx%02d",i-1)));
193       p->SetMarkerColor(kBlue+1);
194       p->SetLineColor(kBlue+1);
195       p->SetMarkerSize(0.8);
196       p->SetMarkerStyle(20);
197       p->DrawClone("same");
198       p->Write();
199
200       // Make the ratio of the two histograms and draw it 
201       pad->cd(2);
202       TH1* ratio = static_cast<TH1*>(vh1->Clone(Form("ratio_vtx%02d", i-1)));
203       // ratio->Scale(1./fEvents->GetBinContent(i));
204       ratio->Divide(p);
205       ratio->SetMarkerSize(.8);
206       ratio->SetLineColor(kGray);
207       ratio->DrawClone("P");
208     }
209     cVtx->cd();
210     savdir->cd();
211
212     // Scale the histogram summed over the vertex bins.  This must be 
213     // divided by the number of bins we have summed over.  If we do 
214     // rebinning, we should scale it again. 
215     fTotal1D->Divide(fNorm);
216     fTotal1D->Scale(1, "width"); 
217     return fTotal1D;
218   }
219   /** 
220    * Browse this object 
221    * 
222    * @param b Browser to use 
223    */
224   void Browse(TBrowser* b) 
225   {
226     b->Add(&fByVtx1D);
227     DrawRes1D::Browse(b);
228   }
229   /** 
230    * Create a new object of this class, and add it to the list of 
231    * specials, and create a browse and browse to this object 
232    * 
233    * @return Newly created object
234    */
235   static DrawResVtx* Create() 
236   {
237     DrawResVtx* dr = new DrawResVtx;
238     gROOT->GetListOfSpecials()->Add(dr);
239     TBrowser* b = new TBrowser("b");
240     b->BrowseObject(gROOT->GetListOfSpecials());
241     return dr;
242   }
243
244 protected:
245   TObjArray          fByVtx1D;   // List of per-vertex 1D histograms
246
247   ClassDef(DrawResVtx,0)
248 };
249
250 //
251 // EOF
252 //