PWG2/FORWARD/analysis2/scripts/CompareCentralSecMaps.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 CompareCentralSecMaps(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   TFile file1(fn1);
  30   TFile file2(fn2);
  31
  32   AliCentralMultiplicityTask task("tmp");
  33
  34   o1 = file1.Get(task.GetManager().GetSecMapName());
  35   o2 = file2.Get(task.GetManager().GetSecMapName());
  36   if (!o1 || !o2) return;
  37   AliCentralCorrSecondaryMap* obj1 = static_cast<AliCentralCorrSecondaryMap*>(o1);
  38   AliCentralCorrSecondaryMap* obj2 = static_cast<AliCentralCorrSecondaryMap*>(o2);
  39   UShort_t nVtx = obj1->GetVertexAxis().GetNbins();
  40
  41   // --- Make canvas -------------------------------------------------
  42   // Canvas* c = new Canvas("CentralSecMapComparison", "Ratio of central secondary maps", n1, n2);
  43   // c->Open();
  44
  45   // --- Loop over the data ------------------------------------------
  46   //for (UShort_t d = 1; d <= 3; d++) {
  47   //  UShort_t nR = (d == 1 ? 1 : 2);
  48   //  for (UShort_t q = 0; q < nR; q++) {
  49   //   Char_t   r  = (q == 0 ? 'I' : 'O');
  50   //    UShort_t nS = (q == 0 ?  20 :  40);
  51
  52       // --- Make 2D ratios ------------------------------------------
  53       //c->Clear(nVtx, d, r);
  54   TCanvas* c2D = new TCanvas("2Dcomparison","2Dcomparison",800,1000);
  55   TCanvas* c1D = new TCanvas("1Dcomparison","1Dcomparison",800,1000);
  56   c2D->Divide(2,5);
  57   c1D->Divide(2,5);
  58   c2D->cd();
  59
  60   TList hists;
  61   for (UShort_t v=1; v <= nVtx; v++) {
  62     //  TVirtualPad* p = c->cd(v);
  63     c2D->cd(v);
  64     TH2* h1 = obj1->GetCorrection( v);
  65     TH2* h2 = obj2->GetCorrection( v);
  66     //std::cout<<h1<<"  "<<h2<<std::endl;
  67     /*  if (!h1) {
  68         Error("CompareSecMaps",
  69         "Map for SPD, vtxbin %3d not found in first",
  70         v);
  71         continue;
  72         }
  73         if (!h2) {
  74         Error("CompareSecMaps",
  75         "Map for SPD, vtxbin %3d not found in second",
  76         v);
  77         continue;
  78         }
  79     */
  80     Double_t vl    = obj1->GetVertexAxis().GetBinLowEdge(v);
  81     Double_t vh    = obj1->GetVertexAxis().GetBinUpEdge(v);
  82     TH2*     ratio =
  83       static_cast<TH2*>(h1->Clone(Form("tmpSPD_%3d",v)));
  84     ratio->SetName(Form("SPD_vtx%03d_ratio",  v));
  85     ratio->SetTitle(Form("%+5.1f<v_{z}<%-+5.1f", vl, vh));
  86     ratio->Divide(h2);
  87     ratio->SetStats(0);
  88     ratio->SetDirectory(0);
  89     ratio->SetZTitle("ratio");
  90
  91         //if (ratio->GetMaximum()-ratio->GetMinimum() > 10)
  92         //  p->SetLogz();
  93
  94     ratio->DrawCopy("colz");
  95     hists.AddAt(ratio, v-1);
  96
  97   }
  98       // c->Print(d, r);
  99
 100       // --- Make 1D profiles ----------------------------------------
 101       //c->Clear(nVtx, d, r);
 102
 103   c1D->cd();
 104   for (UShort_t v=1; v <= nVtx; v++) {
 105     //c->cd(v);
 106     c1D->cd(v);
 107     TH2* hist = static_cast<TH2*>(hists.At(v-1));
 108     TH1D* prof = hist->ProjectionX();
 109     prof->Clear();
 110     for(Int_t i=1; i<=hist->GetNbinsX();i++) {
 111
 112       Float_t sum = 0;
 113       Float_t error = 0;
 114       Int_t   n   = 0;
 115       for(Int_t j=1; j<=hist->GetNbinsY();j++) {
 116         if(hist->GetBinContent(i,j) > 0) {
 117           sum = sum+hist->GetBinContent(i,j);
 118           error = error + TMath::Power(hist->GetBinError(i,j),2);
 119           n++;
 120         }
 121
 122       }
 123       if(n>0) {
 124         sum = sum/(Float_t)n;
 125         error = TMath::Sqrt(error) / (Float_t)n;
 126         prof->SetBinContent(i,sum);
 127         prof->SetBinError(i,error);
 128       }
 129     }
 130
 131     //    prof->Scale(0.05);
 132     prof->SetStats(0);
 133     prof->SetMinimum(0.8);
 134     prof->SetMaximum(1.2);
 135
 136
 137     prof->Fit("pol0","Q");
 138     prof->DrawCopy();
 139
 140     TF1* f = prof->GetFunction("pol0");
 141
 142     TLatex* l = new TLatex(0.5, 0.4, Form("A = %f #pm %f",
 143                                           f->GetParameter(0),
 144                                           f->GetParError(0)));
 145     l->SetTextAlign(22);
 146     l->SetNDC();
 147     l->Draw();
 148     l->DrawLatex(0.5, 0.3, Form("#chi^2/NDF = %f / %d = %f",
 149                         f->GetChisquare(),
 150                         f->GetNDF(),
 151                         f->GetChisquare() / f->GetNDF()));
 152     Double_t dist = TMath::Abs(1 - f->GetParameter(0));
 153     l->DrawLatex(0.5, 0.35, Form("|1 - A| = %f %s #deltaA",
 154          dist, dist <= f->GetParError(0) ?
 155                  "#leq" : ">"));
 156
 157     TLine* l1 = new TLine(-4, 1, 6, 1);
 158     l1->SetLineColor(kRed);
 159     l1->SetLineStyle(2);
 160     l1->Draw();
 161
 162   }
 163
 164   //c->Print(d, r, "profiles");
 165   c2D->Print("comparisonSPD.pdf(");
 166   c1D->Print("comparisonSPD.pdf)");
 167
 168
 169   // --- Close stuff -------------------------------------------------
 170   //c->Close();
 171   // file1->Close();
 172   // file2->Close();
 173 }
 174
 175
 176 //____________________________________________________________________
 177 //
 178 // EOF
 179 //