AliPhysics  d565ceb (d565ceb)
CompareCentralSecMaps.C
Go to the documentation of this file.
1 
9 //____________________________________________________________________
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_PHYSICS/PWGLF/FORWARD/analysis2/scripts/LoadLibs.C");
30  // gROOT->LoadMacro("$ALICE_PHYSICS/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;
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 //
virtual AliCorrectionManagerBase * GetManager() const
double Double_t
Definition: External.C:58
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
Definition: External.C:212
TH2D * GetCorrection(Double_t v) const
void CompareCentralSecMaps(const char *fn1, const char *fn2, const char *n1=0, const char *n2=0, bool load=true)
Definition: External.C:220
const TAxis & GetVertexAxis() const
unsigned short UShort_t
Definition: External.C:28