PWGLF/FORWARD/analysis2/corrs/CompareSecMaps.C

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