82fbe2a6aaf029cd2daf340398c0db59558e8631
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / scripts / CompareSecMaps.C
1 /** 
2  * Make ratio of two specific maps 
3  * 
4  * @param d        Detector
5  * @param r        Ring
6  * @param v        Vertex bin (1 based)
7  * @param first    First correction
8  * @param second   Second correction
9  * 
10  * @return Ratio of the two, or null
11  */
12 TH2*
13 Compare2Maps(UShort_t d, Char_t r, UShort_t v, 
14              const AliFMDCorrSecondaryMap& first, 
15              const AliFMDCorrSecondaryMap& second)
16 {
17   TH2* h1 = first.GetCorrection(d, r, v);
18   TH2* h2 = second.GetCorrection(d, r, v);
19   
20   if (!h1) { 
21     Error("Compare2Maps", "Map for FMD%d%c, vtxbin %3d not found in first", 
22           d, r, v);
23     return 0;
24   }
25   if (!h1) { 
26     Error("Compare2Maps", "Map for FMD%d%c, vtxbin %3d not found in second", 
27           d, r, v);
28     return 0;
29   }
30   
31   Double_t vl    = first.GetVertexAxis().GetBinLowEdge(v);
32   Double_t vh    = first.GetVertexAxis().GetBinUpEdge(v);
33   TH2*     ratio = static_cast<TH2*>(h1->Clone(Form("tmpFMD%d%c_%3d",d,r,v)));
34   ratio->SetName(Form("FMD%d%c_vtx%03d_ratio", d, r, v));
35   ratio->SetTitle(Form("%+5.1f<v_{z}<%-+5.1f", vl, vh));
36   ratio->Divide(h2);
37   ratio->SetStats(0);
38   ratio->SetDirectory(0);
39   ratio->SetZTitle("ratio");
40   // ratio->SetMinimum(0.9);
41   // ratio->SetMaximum(1.5);
42
43   return ratio;
44 }
45
46 //____________________________________________________________________
47 TPad* 
48 ClearCanvas(TCanvas* c, UShort_t nVtx, UShort_t d, Char_t r, 
49             const char* n1, const char* n2)
50 {
51   c->Clear();
52   TPad* p1 = new TPad("top", "Top", 0, .95, 1, 1, 0, 0);
53   p1->Draw();
54   p1->cd();
55
56   TLatex* l = new TLatex(.5, .5, Form("Ratio of secondary maps for "
57                                        "FMD%d%c (%s / %s)", d, r, n1, n2));
58   l->SetNDC();
59   l->SetTextAlign(22);
60   l->SetTextSize(0.3);
61   l->Draw();
62   
63   c->cd();
64   TPad* body = new TPad("body", "Body", 0, 0, 1, .95, 0, 0);
65   body->SetTopMargin(0.05);
66   body->SetRightMargin(0.05);
67   body->Divide(2, (nVtx+1)/2, 0.001, 0.001);
68   body->Draw();
69
70   return body;
71 }
72
73 //____________________________________________________________________
74 void
75 CompareSecMaps(const char* fn1,   const char* fn2, 
76                const char* n1=0,  const char* n2=0)
77 {
78
79   // --- Load libraries ----------------------------------------------
80   // gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
81   
82   // --- Open files --------------------------------------------------
83   const char* nam1 = n1;
84   const char* nam2 = n2;
85   if (!n1) nam1 = fn1;                          
86   if (!n2) nam2 = fn2;
87
88   TFile* file1 = TFile::Open(fn1, "READ");
89   TFile* file2 = TFile::Open(fn2, "READ");
90
91   if (!file1) { 
92     Error("CompareSecMaps", "File %s cannot be opened for %s", fn1, n1);
93     return;
94   }
95
96   if (!file2) { 
97     Error("CompareSecMaps", "File %s cannot be opened for %s", fn2, n2);
98     return;
99   }
100   
101   // --- Find Objects ------------------------------------------------
102   const char* objName = AliForwardCorrectionManager::Instance()
103     .GetObjectName(AliForwardCorrectionManager::kSecondaryMap);
104   
105   AliFMDCorrSecondaryMap* obj1 = 
106     static_cast<AliFMDCorrSecondaryMap*>(file1->Get(objName));
107   AliFMDCorrSecondaryMap* obj2 = 
108     static_cast<AliFMDCorrSecondaryMap*>(file2->Get(objName));
109   
110   if (!obj1) {
111     Error("CompareSecMaps", "File %s does not contain an object named %s", 
112           fn1, objName);
113     return;
114   }
115   if (!obj2) {
116     Error("CompareSecMaps", "File %s does not contain an object named %s", 
117           fn2, objName);
118     return;
119   }
120   UShort_t nVtx = obj1->GetVertexAxis().GetNbins();
121
122   // --- Make canvas -------------------------------------------------
123   const char* pdfName = "secMapComparison.pdf";
124   gStyle->SetPalette(1);
125   gStyle->SetTitleX(.10);
126   gStyle->SetTitleY(.99);
127   gStyle->SetTitleW(.85);
128   gStyle->SetTitleH(.085);
129   gStyle->SetTitleFillColor(0);
130   gStyle->SetTitleBorderSize(0);
131
132   TCanvas* c = new TCanvas("c", "c", 800, TMath::Sqrt(2)*800);
133   c->SetFillColor(0);
134   
135   c->Print(Form("%s[", pdfName), "pdf");
136
137   for (UShort_t d = 1; d <= 3; d++) { 
138     UShort_t nR = (d == 1 ? 1 : 2);
139     for (UShort_t q = 0; q < nR; q++) { 
140       Char_t   r  = (q == 0 ? 'I' : 'O');
141       UShort_t nS = (q == 0 ?  20 :  40);
142
143       TPad* body = ClearCanvas(c, nVtx, d, r, nam1, nam2);
144       TList hists;
145       for (UShort_t v=1; v <= nVtx; v++) { 
146         TVirtualPad* p = body->cd(v);
147         // p->SetTopMargin(0.1);
148         // p->SetBottomMargin(0.05);
149         // p->SetRightMargin(0.05);
150         
151         TH2* ratio = Compare2Maps(d, r, v, *obj1, *obj2);
152         if (ratio->GetMaximum()-ratio->GetMinimum() > 10) 
153           p->SetLogz();
154
155         ratio->Draw("colz");
156         hists.AddAt(ratio, v-1);
157       }
158       c->Print(pdfName, Form("Title:FMD%d%c", d, r));
159       
160       body = ClearCanvas(c, nVtx, d, r, nam1, nam2);
161
162       for (UShort_t v=1; v <= nVtx; v++) { 
163         TVirtualPad* p    = body->cd(v);
164         TH2*         hist = static_cast<TH2*>(hists.At(v-1));
165         TH1*         prof = hist->ProjectionX();
166         prof->Scale(1. / nS);
167         prof->SetStats(0);
168         prof->SetMinimum(0.8);
169         prof->SetMaximum(1.2);
170
171         // prof->Draw("hist");
172         // prof->DrawCopy("e same");
173         prof->Draw();
174         prof->Fit("pol0","Q");
175
176         TF1* f = prof->GetFunction("pol0");
177
178         TLatex* l = new TLatex(0.5, 0.4, Form("A = %f #pm %f", 
179                                               f->GetParameter(0), 
180                                               f->GetParError(0)));
181         l->SetTextAlign(22);
182         l->SetNDC();
183         l->Draw();
184         l->DrawLatex(0.5, 0.3, Form("#chi^2/NDF = %f / %d = %f", 
185                                     f->GetChisquare(), 
186                                     f->GetNDF(), 
187                                     f->GetChisquare() / f->GetNDF()));
188         Double_t dist = TMath::Abs(1 - f->GetParameter(0));
189         l->DrawLatex(0.5, 0.35, Form("|1 - A| = %f %s #deltaA", 
190                                      dist, dist <= f->GetParError(0) ? 
191                                      "#leq" : ">")); 
192
193         TLine* l1 = new TLine(-4, 1, 6, 1);
194         l1->SetLineColor(kRed);
195         l1->SetLineStyle(2);
196         l1->Draw();
197       }
198
199       c->Print(pdfName, Form("Title:FMD%d%c profiles", d, r));
200     }
201   }
202
203   c->Print(Form("%s]", pdfName), "pdf");
204   file1->Close();
205   file2->Close();
206 }
207
208