Scripts/MaterialMapping/Mat_map.C

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include <TROOT.h>
0010
0011 #include <string>
0012 #include <fstream>
0013 #include <iostream>
0014 #include <sstream>
0015
0016 using namespace std ;
0017
0018 /// Function used to make ratio plot.
0019
0020 void Draw_ratio(TCanvas* c, TProfile* h1, TProfile* h2, TLegend* leg, std::string name){
0021
0022   // Upper plot will be in pad1
0023   TPad *pad1 = new TPad("pad1", "pad1", 0, 0.3, 1, 1.0);
0024
0025   pad1->SetBottomMargin(0); // Upper and lower plot are joined
0026   pad1->SetGridx();         // Vertical grid
0027   pad1->Draw();             // Draw the upper pad: pad1
0028   pad1->cd();               // pad1 becomes the current pad
0029   h1->Draw("E");
0030   h2->Draw("E SAME");
0031   leg->Draw("SAME");
0032
0033   double max_hist [4];
0034   max_hist[0] = h1->GetMaximum()+h1->GetMaximum()*0.05;
0035   max_hist[1] = h2->GetMaximum()+h2->GetMaximum()*0.05;
0036
0037   h1->GetYaxis()->SetTitleOffset(1.5);
0038   h1->GetYaxis()->SetRangeUser(0, *max_element( begin(max_hist),end(max_hist) ) );
0039
0040   // lower plot will be in pad
0041   c->cd();          // Go back to the main canvas before defining pad2
0042   TPad *pad2 = new TPad("pad2", "pad2", 0, 0.05, 1, 0.3);
0043   pad2->SetTopMargin(0);
0044   pad2->SetBottomMargin(0.2);
0045   pad2->SetGridx(); // vertical grid
0046   pad2->Draw();
0047   pad2->cd();       // pad2 becomes the current pad
0048
0049   TLine line = TLine(h1->GetXaxis()->GetXmin(),1,h1->GetXaxis()->GetXmax(),1);
0050   line.SetLineColor(kRed);
0051   line.SetLineWidth(1);
0052
0053   // Define the ratio plot
0054   TProfile *h5 = (TProfile*)h2->Clone( (name+"_ratio").c_str() );
0055   h5->SetLineColor(kBlack);
0056
0057   h5->SetStats(0);      // No statistics on lower plot
0058   h5->Divide(h1);
0059
0060   double maxi = min( max( std::abs(h5->GetMinimum()-0.1*h5->GetMinimum()),h5->GetMaximum()+0.1*h5->GetMaximum() ), 10. );
0061
0062   h5->SetMinimum( 0.5 );  // Define Y ..
0063   h5->SetMaximum( 1.1 ); // .. range
0064
0065   h5->SetMarkerStyle(7);
0066   h5->Draw("hist P");       // Draw the ratio plot
0067   line.Draw("SAME");
0068
0069   // Y axis ratio plot settings
0070   h5->GetYaxis()->SetTitle("Ratio Validation/Geantino ");
0071   h5->GetYaxis()->SetNdivisions(505);
0072   h5->GetYaxis()->SetTitleSize(15);
0073   h5->GetYaxis()->SetTitleFont(43);
0074   h5->GetYaxis()->SetTitleOffset(1.55);
0075   h5->GetYaxis()->SetLabelFont(43);
0076   h5->GetYaxis()->SetLabelSize(15);
0077
0078   // X axis ratio plot settings
0079   h5->GetXaxis()->SetTitleSize(20);
0080   h5->GetXaxis()->SetTitleFont(43);
0081   h5->GetXaxis()->SetTitleOffset(3);
0082   h5->GetXaxis()->SetLabelFont(43);
0083   h5->GetXaxis()->SetLabelSize(15);
0084
0085   std::string name_axis = h1->GetXaxis()->GetTitle();
0086
0087   c->Print( (name+"/Ratio_Val_geant_mat_X0_"+name_axis+".pdf").c_str());
0088
0089   delete h5;
0090
0091 }
0092
0093 /// Compare two set of material tracks (for example one obtain with propagator and material map and one with geantino scan).
0094 /// Draw the ammont of material (in X0) encounter by tracks as function of eta and phi.
0095 /// Plot the ratio between the two set to help identify inconsistency.
0096
0097 void Mat_map(std::string Val = "", std::string geantino = "", std::string name = ""){
0098
0099   gStyle->SetOptStat(0);
0100   gStyle->SetOptTitle(0);
0101
0102   TProfile * Val_X0_Eta = new TProfile("Val_X0_Eta","Val_X0_Eta",160,-4,4);
0103   TProfile * Val_X0_Phi = new TProfile("Val_X0_Phi","Val_X0_Phi",160,-4,4);
0104   TH2F * Val_X0_Eta_spread = new TH2F("Val_X0_Eta_spread","Val_X0_Eta_spread",160,-4,4,160,0,4);
0105   TH2F * Val_X0_Phi_spread = new TH2F("Val_X0_Phi_spread","Val_X0_Phi_spread",160,-4,4,160,0,4);
0106   TH2F * Val_X0 = new TH2F("Val_X0","Val_X0",160,-4,4,160,-4,4);
0107
0108   TProfile * geantino_X0_Eta = new TProfile("geantino_X0_Eta","geantino_X0_Eta",160,-4,4);
0109   TProfile * geantino_X0_Phi = new TProfile("geantino_X0_Phi","geantino_X0_Phi",160,-4,4);
0110   TH2F * geantino_X0_Eta_spread = new TH2F("geantino_X0_Eta_spread","geantino_X0_Eta_spread",160,-4,4,160,0,4);
0111   TH2F * geantino_X0_Phi_spread = new TH2F("geantino_X0_Phi_spread","geantino_X0_Phi_spread",160,-4,4,160,0,4);
0112   TH2F * geantino_X0 = new TH2F("geantino_X0","geantinol_X0",160,-4,4,160,-4,4);
0113
0114   TH1F * comp_X0_Eta = new TH1F("comp_X0_Eta","comp_X0_Eta",160,-4,4);
0115   TH1F * comp_X0_Phi = new TH1F("comp_X0_Phi","comp_X0_Phi",160,-4,4);
0116   TH2F * comp_X0 = new TH2F("comp_X0","comp_X0",160,-4,4,160,-4,4);
0117
0118   TChain *Val_file = new TChain("material-tracks");
0119   TChain *geantino_file = new TChain("material-tracks");
0120
0121   // Define line corresponding to the different eta value
0122   TLine *eta_0 = new TLine(0,-1200,0,1200);
0123   eta_0->SetLineColor(kRed);
0124
0125   TLine *eta_1p = new TLine(-1250,-1064,1250,1064);
0126   eta_1p->SetLineColor(kRed);
0127   TLine *eta_2p = new TLine(-3000,-827,3000,827);
0128   eta_2p->SetLineColor(kRed);
0129   TLine *eta_3p = new TLine(-3000,-300,3000,300);
0130   eta_3p->SetLineColor(kRed);
0131   TLine *eta_4p = new TLine(-3000,-110,3000,110);
0132   eta_4p->SetLineColor(kRed);
0133
0134   TLine *eta_1n = new TLine(-1250,1064,1250,-1064);
0135   eta_1n->SetLineColor(kRed);
0136   TLine *eta_2n = new TLine(-3000,827,3000,-827);
0137   eta_2n->SetLineColor(kRed);
0138   TLine *eta_3n = new TLine(-3000,300,3000,-300);
0139   eta_3n->SetLineColor(kRed);
0140   TLine *eta_4n = new TLine(-3000,110,3000,-110);
0141   eta_4n->SetLineColor(kRed);
0142
0143   if(Val != ""){
0144     Val_file->Add(Val.c_str());
0145
0146     // 2D map for Validation input
0147     TCanvas *VM = new TCanvas("VM","Validation Map") ;
0148     Val_file->Draw("mat_y:mat_z","std::abs(mat_x)<1");
0149
0150     eta_0->Draw("Same");
0151     eta_1p->Draw("Same");
0152     eta_1n->Draw("Same");
0153     eta_2p->Draw("Same");
0154     eta_2n->Draw("Same");
0155     eta_3p->Draw("Same");
0156     eta_3n->Draw("Same");
0157     eta_4p->Draw("Same");
0158     eta_4n->Draw("Same");
0159
0160     VM->Print( (name+"/Val_mat_map.png").c_str());
0161     //VM->Print( (name+"/Val_mat_map.pdf").c_str());
0162
0163     // X0 as function of Eta for Validation input
0164     TCanvas *VM_X0_Eta = new TCanvas("VM_X0_Eta","Validation X0 Eta") ;
0165     Val_file->Draw("t_X0:v_eta>>Val_X0_Eta","","profile");
0166     Val_X0_Eta->SetMarkerStyle(7);
0167     Val_X0_Eta->Draw("HIST PC");
0168     Val_X0_Eta->GetXaxis()->SetTitle("Eta");
0169     Val_X0_Eta->GetYaxis()->SetTitle("X0");
0170     VM_X0_Eta->Print( (name+"/Val_mat_Eta_X0.pdf").c_str());
0171
0172     // X0 as function of Eta for Validation input
0173     TCanvas *VM_X0_Eta_spread = new TCanvas("VM_X0_Eta_spread","Validation X0 Eta") ;
0174     Val_file->Draw("t_X0:v_eta>>Val_X0_Eta_spread","","");
0175     Val_X0_Eta_spread->GetXaxis()->SetTitle("Eta");
0176     Val_X0_Eta_spread->GetYaxis()->SetTitle("X0");
0177     VM_X0_Eta_spread->Print( (name+"/Val_X0_Eta_spread.pdf").c_str());
0178
0179     // X0 as function of Phi for Validation input
0180     TCanvas *VM_X0_Phi = new TCanvas("VM_X0_Phi","Validation X0 Phi") ;
0181     Val_file->Draw("t_X0:v_phi>>Val_X0_Phi","","profile");
0182     Val_X0_Phi->SetMarkerStyle(7);
0183     Val_X0_Phi->Draw("HIST PC");
0184     Val_X0_Phi->GetXaxis()->SetTitle("Phi");
0185     Val_X0_Phi->GetYaxis()->SetTitle("X0");
0186     VM_X0_Phi->Print( (name+"/Val_mat_Phi_X0.pdf").c_str());
0187
0188     // X0 as function of Phi for Validation input
0189     TCanvas *VM_X0_Phi_spread = new TCanvas("VM_X0_Phi_spread","Validation X0 Phi") ;
0190     Val_file->Draw("t_X0:v_phi>>Val_X0_Phi_spread","","");
0191     Val_X0_Phi_spread->GetXaxis()->SetTitle("Phi");
0192     Val_X0_Phi_spread->GetYaxis()->SetTitle("X0");
0193     VM_X0_Phi_spread->Print( (name+"/Val_mat_Phi_X0_spread.pdf").c_str());
0194
0195     // 2D map of X0 for Validation input
0196     TCanvas *VM_2D = new TCanvas("VM_2D","Validation X0 2D") ;
0197     Val_file->Draw("v_phi:v_eta:t_X0>>Val_X0","","COLZ");
0198     Val_X0->GetXaxis()->SetTitle("Eta");
0199     Val_X0->GetYaxis()->SetTitle("Phi");
0200     Val_X0->GetZaxis()->SetTitle("X0");
0201     VM_2D->Print( (name+"/Val_mat_X0.png").c_str());
0202   }
0203
0204   if(geantino != ""){
0205     geantino_file->Add(geantino.c_str());
0206
0207     // 2D map for Geantino input
0208     TCanvas *GM = new TCanvas("GM","Geantino Map") ;
0209     geantino_file->Draw("mat_y:mat_z","std::abs(mat_x)<1");
0210
0211     eta_0->Draw("Same");
0212     eta_1p->Draw("Same");
0213     eta_1n->Draw("Same");
0214     eta_2p->Draw("Same");
0215     eta_2n->Draw("Same");
0216     eta_3p->Draw("Same");
0217     eta_3n->Draw("Same");
0218     eta_4p->Draw("Same");
0219     eta_4n->Draw("Same");
0220
0221     GM->Print( (name+"/geant_mat_map.png").c_str());
0222     //GM->Print( (name+"/geant_mat_map.pdf").c_str());
0223
0224     // X0 as function of Eta for Geantino input
0225     TCanvas *GM_X0_Eta = new TCanvas("GM_X0_Eta","Geantino X0 Eta") ;
0226     geantino_file->Draw("t_X0:v_eta>>geantino_X0_Eta","","profile");
0227     geantino_X0_Eta->SetMarkerStyle(7);
0228     geantino_X0_Eta->Draw("HIST PC");
0229     geantino_X0_Eta->GetXaxis()->SetTitle("Eta");
0230     geantino_X0_Eta->GetYaxis()->SetTitle("X0");
0231     GM_X0_Eta->Print( (name+"/geant_mat_Eta_X0.pdf").c_str());
0232
0233     // X0 as function of Eta for Geantino input
0234     TCanvas *GM_X0_Eta_spread = new TCanvas("GM_X0_Eta_spread","Geantino X0 Eta") ;
0235     geantino_file->Draw("t_X0:v_eta>>geantino_X0_Eta_spread","","");
0236     geantino_X0_Eta_spread->GetXaxis()->SetTitle("Eta");
0237     geantino_X0_Eta_spread->GetYaxis()->SetTitle("X0");
0238     GM_X0_Eta_spread->Print( (name+"/geant_mat_Eta_X0_spread.pdf").c_str());
0239
0240     // X0 as function of Phi for Geantino input
0241     TCanvas *GM_X0_Phi = new TCanvas("GM_X0_Phi","Geantino X0 Phi") ;
0242     geantino_file->Draw("t_X0:v_phi>>geantino_X0_Phi","","profile");
0243     geantino_X0_Phi->SetMarkerStyle(7);
0244     geantino_X0_Phi->Draw("HIST PC");
0245     geantino_X0_Phi->GetXaxis()->SetTitle("Phi");
0246     geantino_X0_Phi->GetYaxis()->SetTitle("X0");
0247     GM_X0_Phi->Print( (name+"/geant_mat_Phi_X0.pdf").c_str());
0248
0249     // X0 as function of Phi for Geantino input
0250     TCanvas *GM_X0_Phi_spread = new TCanvas("GM_X0_Phi_spread","Geantino X0 Phi") ;
0251     geantino_file->Draw("t_X0:v_phi>>geantino_X0_Phi_spread","","");
0252     geantino_X0_Phi_spread->GetXaxis()->SetTitle("Phi");
0253     geantino_X0_Phi_spread->GetYaxis()->SetTitle("X0");
0254     GM_X0_Phi_spread->Print( (name+"/geant_mat_Phi_X0_spread.pdf").c_str());
0255
0256     // 2D map of X0 for Geantino input
0257     TCanvas *GM_2D = new TCanvas("GM_2D","Geantino X0 2D") ;
0258     geantino_file->Draw("v_phi:v_eta:t_X0>>geantino_X0","","COLZ");
0259     geantino_X0->GetXaxis()->SetTitle("Eta");
0260     geantino_X0->GetYaxis()->SetTitle("Phi");
0261     geantino_X0->GetZaxis()->SetTitle("X0");
0262     GM_2D->Print( (name+"/geant_mat_X0.png").c_str());
0263   }
0264
0265   if(Val != "" && geantino != ""){
0266
0267     Val_X0_Eta->SetMarkerColor(kBlack);
0268     Val_X0_Eta->SetLineColor(kBlack);
0269     geantino_X0_Eta->SetMarkerColor(kRed);
0270     geantino_X0_Eta->SetLineColor(kRed);
0271
0272     Val_X0_Phi->SetMarkerColor(kBlack);
0273     Val_X0_Phi->SetLineColor(kBlack);
0274     geantino_X0_Phi->SetMarkerColor(kRed);
0275     geantino_X0_Phi->SetLineColor(kRed);
0276
0277     TLegend* leg = new TLegend(0.1,0.15,0.25,0.30);
0278     leg->AddEntry(Val_X0_Eta,"Validation X0");
0279     leg->AddEntry(geantino_X0_Eta,"Geantino X0");
0280
0281     TLegend* leg2 = new TLegend(0.1,0.15,0.25,0.30);
0282     leg2->AddEntry(Val_X0_Phi,"Validation X0");
0283     leg2->AddEntry(geantino_X0_Phi,"Geantino X0");
0284
0285     // X0 differences as function of eta of the Validation and Geantino input
0286     comp_X0_Eta->Add(Val_X0_Eta);
0287     comp_X0_Eta->Add(geantino_X0_Eta,-1);
0288     comp_X0_Eta->Divide(Val_X0_Eta);
0289     comp_X0_Eta->Scale(100);
0290     comp_X0_Eta->GetXaxis()->SetTitle("Eta");
0291     comp_X0_Eta->GetYaxis()->SetTitle("(Validation X0 - geantino X0) / Validation X0 [%]");
0292     TCanvas *Diff_Eta = new TCanvas("Diff_Eta","Comparison geantino_Validation Eta") ;
0293     comp_X0_Eta->Draw();
0294     Diff_Eta->Print( (name+"/Comp_Val_geant_mat_X0_Eta.pdf").c_str());
0295
0296     // X0 comparison as function of eta of the Validation and Geantino input
0297     TCanvas *Comp_Eta_spread = new TCanvas("Comp_Eta_spread","Comparison geantino_Validation Eta") ;
0298     Val_X0_Eta_spread->Draw();
0299     geantino_X0_Eta_spread->SetMarkerColor(kRed);
0300     geantino_X0_Eta_spread->Draw("SAME");
0301     leg->Draw("SAME");
0302     Comp_Eta_spread->Print( (name+"/Comp_Val_geant_mat_X0_Eta_spread.pdf").c_str());
0303
0304     // X0 differences as function of phi of the Validation and Geantino input
0305     comp_X0_Phi->Add(Val_X0_Phi);
0306     comp_X0_Phi->Add(geantino_X0_Phi,-1);
0307     comp_X0_Phi->Divide(Val_X0_Phi);
0308     comp_X0_Phi->Scale(100);
0309     comp_X0_Phi->GetXaxis()->SetTitle("Phi");
0310     comp_X0_Phi->GetYaxis()->SetTitle("(Validation X0 - geantino X0) / Validation X0 [%]");
0311     TCanvas *Diff_Phi = new TCanvas("Diff_Phi","Comparison geantino_Validation") ;
0312     comp_X0_Phi->Draw();
0313     Diff_Phi->Print( (name+"/Comp_Val_geant_mat_X0_Phi.pdf").c_str());
0314
0315     // X0 comparison as function of eta of the Validation and Geantino input
0316     TCanvas *Comp_Phi_spread = new TCanvas("Comp_Phi_spread","Comparison geantino_Validation Phi") ;
0317     Val_X0_Phi_spread->Draw();
0318     geantino_X0_Phi_spread->SetMarkerColor(kRed);
0319     geantino_X0_Phi_spread->Draw("SAME");
0320     leg2->Draw("SAME");
0321     Comp_Phi_spread->Print( (name+"/Comp_Val_geant_mat_X0_Phi_spread.pdf").c_str());
0322
0323     Float_t score = 0;
0324     for(int i=0; i<Val_X0->GetXaxis()->GetNbins(); i++){
0325       score += (Val_X0->GetBinContent(i+1)-geantino_X0->GetBinContent(i+1))*
0326       (Val_X0->GetBinContent(i+1)-geantino_X0->GetBinContent(i+1))/
0327       geantino_X0->GetBinError(i+1)*geantino_X0->GetBinError(i+1);
0328     }
0329     std::cout << "Score : " << score << std::endl;
0330
0331     TText *t = new TText(0,10, ("Chi^2 score : " + to_string(score)).c_str());
0332     t->Draw();
0333
0334     // X0 ratio plot as function of eta of the Validation and Geantino input
0335     TCanvas *Comp_Eta = new TCanvas("Comp_Eta","Ratio geantino_Validation Eta") ;
0336     Val_X0_Eta->SetMarkerStyle(7);
0337     geantino_X0_Eta->SetMarkerStyle(7);
0338     Val_X0_Eta->SetMarkerColor(kBlack);
0339     geantino_X0_Eta->SetMarkerColor(kRed);
0340     Val_X0_Eta->SetLineColor(kBlack);
0341     geantino_X0_Eta->SetLineColor(kRed);
0342
0343     Draw_ratio(Comp_Eta, geantino_X0_Eta, Val_X0_Eta, leg, name);
0344
0345     // X0 ratio plot as function of phi of the Validation and Geantino input
0346     TCanvas *Comp_Phi = new TCanvas("Comp_Phi","Ratio geantino_Validation Phi") ;
0347     Val_X0_Phi->SetMarkerStyle(7);
0348     geantino_X0_Phi->SetMarkerStyle(7);
0349     Val_X0_Phi->SetMarkerColor(kBlack);
0350     geantino_X0_Phi->SetMarkerColor(kRed);
0351     Val_X0_Phi->SetLineColor(kBlack);
0352     geantino_X0_Phi->SetLineColor(kRed);
0353
0354     Draw_ratio(Comp_Phi, geantino_X0_Phi, Val_X0_Phi, leg, name);
0355   }
0356
0357   return;
0358
0359 }