PWGLF/FORWARD/analysis2/corrs/CompareCentralSecMaps.C

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