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