Renames and new scripts
[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
13 //____________________________________________________________________
14 void
15 CompareSecMaps(const char* fn1,   const char* fn2, 
16                const char* n1=0,  const char* n2=0,
17                bool load=true)
18 {
19
20   // --- Load Utilities ----------------------------------------------
21   if (load) {
22     gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
23     gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/CompareCorrs.C");
24   }
25
26   // --- Get Objects -------------------------------------------------
27   TObject* o1 = 0;
28   TObject* o2 = 0;
29   GetObjects(AliForwardCorrectionManager::kSecondaryMap, fn1, fn2, o1, o2);
30   if (!o1 || !o2) return; 
31   AliFMDCorrSecondaryMap* obj1 = static_cast<AliFMDCorrSecondaryMap*>(o1);
32   AliFMDCorrSecondaryMap* obj2 = static_cast<AliFMDCorrSecondaryMap*>(o2);
33   UShort_t nVtx = obj1->GetVertexAxis().GetNbins();
34
35   // --- Make canvas -------------------------------------------------
36   Canvas* c = new Canvas("secMapComparison", "Ratio of secondary maps", n1, n2);
37   c->Open();
38
39   // --- Loop over the data ------------------------------------------
40   for (UShort_t d = 1; d <= 3; d++) { 
41     UShort_t nR = (d == 1 ? 1 : 2);
42     for (UShort_t q = 0; q < nR; q++) { 
43       Char_t   r  = (q == 0 ? 'I' : 'O');
44       UShort_t nS = (q == 0 ?  20 :  40);
45
46       // --- Make 2D ratios ------------------------------------------
47       c->Clear(nVtx, d, r);
48       TList hists;
49       for (UShort_t v=1; v <= nVtx; v++) { 
50         TVirtualPad* p = c->cd(v);
51         
52         TH2* h1 = obj1->GetCorrection(d, r, v);
53         TH2* h2 = obj1->GetCorrection(d, r, v);
54   
55         if (!h1) { 
56           Error("CompareSecMaps", 
57                 "Map for FMD%d%c, vtxbin %3d not found in first", 
58                 d, r, v);
59           continue;
60         }
61         if (!h2) { 
62           Error("CompareSecMaps", 
63                 "Map for FMD%d%c, vtxbin %3d not found in second", 
64                 d, r, v);
65           continue;
66         }
67   
68         Double_t vl    = obj1->GetVertexAxis().GetBinLowEdge(v);
69         Double_t vh    = obj1->GetVertexAxis().GetBinUpEdge(v);
70         TH2*     ratio = 
71           static_cast<TH2*>(h1->Clone(Form("tmpFMD%d%c_%3d",d,r,v)));
72         ratio->SetName(Form("FMD%d%c_vtx%03d_ratio", d, r, v));
73         ratio->SetTitle(Form("%+5.1f<v_{z}<%-+5.1f", vl, vh));
74         ratio->Divide(h2);
75         ratio->SetStats(0);
76         ratio->SetDirectory(0);
77         ratio->SetZTitle("ratio");
78
79         if (ratio->GetMaximum()-ratio->GetMinimum() > 10) 
80           p->SetLogz();
81
82         ratio->Draw("colz");
83         hists.AddAt(ratio, v-1);
84       }
85       c->Print(d, r);
86       
87       // --- Make 1D profiles ----------------------------------------
88       c->Clear(nVtx, d, r);
89       for (UShort_t v=1; v <= nVtx; v++) { 
90         c->cd(v);
91         TH2* hist = static_cast<TH2*>(hists.At(v-1));
92         TH1* prof = hist->ProjectionX();
93         prof->Scale(1. / nS);
94         prof->SetStats(0);
95         prof->SetMinimum(0.8);
96         prof->SetMaximum(1.2);
97
98         prof->Draw();
99         prof->Fit("pol0","Q");
100
101         TF1* f = prof->GetFunction("pol0");
102
103         TLatex* l = new TLatex(0.5, 0.4, Form("A = %f #pm %f", 
104                                               f->GetParameter(0), 
105                                               f->GetParError(0)));
106         l->SetTextAlign(22);
107         l->SetNDC();
108         l->Draw();
109         l->DrawLatex(0.5, 0.3, Form("#chi^2/NDF = %f / %d = %f", 
110                                     f->GetChisquare(), 
111                                     f->GetNDF(), 
112                                     f->GetChisquare() / f->GetNDF()));
113         Double_t dist = TMath::Abs(1 - f->GetParameter(0));
114         l->DrawLatex(0.5, 0.35, Form("|1 - A| = %f %s #deltaA", 
115                                      dist, dist <= f->GetParError(0) ? 
116                                      "#leq" : ">")); 
117
118         TLine* l1 = new TLine(-4, 1, 6, 1);
119         l1->SetLineColor(kRed);
120         l1->SetLineStyle(2);
121         l1->Draw();
122       }
123
124       c->Print(d, r, "profiles");
125     }
126   }
127
128   // --- Close stuff -------------------------------------------------
129   c->Close();
130   // file1->Close();
131   // file2->Close();
132 }
133
134   
135 //____________________________________________________________________
136 //
137 // EOF
138 //