checkmuons.cpp
Go to the documentation of this file.
1 void checkmuons()
2 {
3  // CaliceStyle();
4 
5  //Compare between ptrue <-> preco and angle <-> anglereco
6 
7  TH1F* hPullMomentum = new TH1F("hPullMomentum", "hPullMomentum", 300, -1, 1);
8  hPullMomentum->SetLineWidth(2);
9  hPullMomentum->SetLineColor(kBlack);
10 
11  TH2F* hMomemtumTruevsReco = new TH2F("hMomemtumTruevsReco", "hMomemtumTruevsReco", 100, 0, 5, 100, 0, 5);
12 
13  TH1F* hPullAngle = new TH1F("hPullAngle", "hPullAngle", 100, -1, 1);
14  hPullAngle->SetLineWidth(2);
15  hPullAngle->SetLineColor(kBlack);
16 
17  TH2F* hAngleTruevsReco = new TH2F("hAngleTruevsReco", "hAngleTruevsReco", 100, -1, 1, 100, -1, 1);
18 
19  TH1F* hPullNeutrons = new TH1F("hPullNeutrons", "hPullNeutrons", 100, -1, 1);
20  hPullNeutrons->SetLineWidth(2);
21  hPullNeutrons->SetLineColor(kBlack);
22 
23  //if you are reading only one file, for debugging purposes: 1) comment out lines 24 through 57 and then lines 94 and 95,
24  //then 2) un-comment lines 59 and 60
25  const char *dirname="/pnfs/dune/persistent/users/ND_GAr/2020_06_21/neutrino/";
26 
27  const char *ext=".root";
28  std::vector<const char*> lista;
29  lista.clear();
30 
31  TChain *chain = new TChain("caf");
32 
33  TSystemDirectory dir(dirname, dirname);
34  TList *files = dir.GetListOfFiles();
35 
36 
37  TSystemFile *file;
38  TString fname;
39  TIter next(files);
40  while ((file = (TSystemFile*)next())) {
41  fname = file->GetName();
42  if (file->IsDirectory() && fname.EndsWith(ext)) continue;
43 
44  char result[255]; // array to hold the result.
45  strcpy(result,dirname); // copy string one into the result.
46  strcat(result,fname.Data()); // append string two to the result.
47 
48  lista.push_back(result);
49 
50  for(std::vector<const char*>::iterator it = lista.begin(); it != lista.end(); it++)
51  {
52 
53 
54  std::cout << *it << std::endl;
55  string a(*it);
56 
57  TFile * tf = new TFile ( Form("%s",a.c_str() ) );
58  chain->Add( Form("%s/caf", a.c_str()) );
59 
60  //TChain *chain = new TChain("caf");
61  //chain->Add("/pnfs/dune/persistent/users/ND_GAr/2020_06_21/neutrino/neutrino.nd_hall_mpd_only.volTPCGas.Ev336000.Ev336999.11141.caf.root");
62 
63  std::vector<float> *truep = 0;
64  std::vector<float> *preco = 0;
65  std::vector<int> *truepdg = 0;
66  std::vector<int> *recopid = 0;
67  std::vector<float> *angle = 0;
68  std::vector<float> *anglereco = 0;
69  std::vector<float> *erecon = 0;
70 
71  chain->SetBranchAddress("truep", &truep);
72  chain->SetBranchAddress("preco", &preco);
73  chain->SetBranchAddress("truepdg", &truepdg);
74  chain->SetBranchAddress("recopid", &recopid);
75  chain->SetBranchAddress("angle", &angle);
76  chain->SetBranchAddress("anglereco", &anglereco);
77  chain->SetBranchAddress("erecon", &erecon);
78 
79 
80  for(int itree = 0; itree < chain->GetEntries(); itree++)
81  {
82  chain->GetEntry(itree);
83 
84  for(int i = 0; i < truep->size(); i++)
85  {
86  //Muons
87  if( std::abs(truepdg->at(i)) == 13 && std::abs(recopid->at(i)) == 13){
88  hPullMomentum->Fill( (truep->at(i) - preco->at(i)) / truep->at(i) );
89  hMomemtumTruevsReco->Fill(truep->at(i), preco->at(i));
90  hPullAngle->Fill( (angle->at(i) - anglereco->at(i)) / angle->at(i) );
91  hAngleTruevsReco->Fill(angle->at(i), anglereco->at(i));
92  }
93  }
94  }
95  }
96  }
97  gStyle->SetOptStat(1110);
98  gStyle->SetOptFit(1);
99  hPullMomentum->Scale(1./hPullMomentum->Integral());
100 
101  hPullMomentum->SetStats(kTRUE);
102  hPullMomentum->Fit("gaus", "EMR+", "", -0.8, 0.8);
103  TCanvas *c1 = new TCanvas("c1", "Pull p", 800, 600);
104  hPullMomentum->GetXaxis()->SetTitle( "(p_{true} - p_{reco}) / p_{true}" );
105  hPullMomentum->GetYaxis()->SetTitle( "Normalized Events" );
106  hPullMomentum->Draw("hist");
107  hPullMomentum->GetFunction("gaus")->SetLineColor(kRed);
108  hPullMomentum->GetFunction("gaus")->SetLineStyle(2);
109  hPullMomentum->GetFunction("gaus")->Draw("same");
110  c1->SaveAs("PullMomentumMuons.pdf");
111 
112  TCanvas *c2 = new TCanvas("c2", "Mom truth vs reco", 800, 600);
113  gPad->SetRightMargin(0.12);
114  hMomemtumTruevsReco->GetXaxis()->SetTitle( "p_{true} [GeV]" );
115  hMomemtumTruevsReco->GetYaxis()->SetTitle( "p_{reco} [GeV]" );
116  hMomemtumTruevsReco->Draw("COLZ");
117  c2->SaveAs("MomentumMuons_TruevsReco.pdf");
118 
119  hPullAngle->Scale(1./hPullAngle->Integral());
120  hPullAngle->SetStats(kTRUE);
121 
122  TCanvas *c3 = new TCanvas("c3", "Pull angle", 800, 600);
123  hPullAngle->GetXaxis()->SetTitle( "(angle_{true} - angle_{reco}) / angle_{true}" );
124  hPullAngle->GetYaxis()->SetTitle( "Normalized Events" );
125  hPullAngle->Draw("hist");
126  c3->SaveAs("PullAngleMuons.pdf");
127 
128  TCanvas *c4 = new TCanvas("c4", "Angle truth vs reco", 800, 600);
129  gPad->SetRightMargin(0.12);
130  hAngleTruevsReco->GetXaxis()->SetTitle( "angle_{true} [GeV]" );
131  hAngleTruevsReco->GetYaxis()->SetTitle( "angle_{reco} [GeV]" );
132  hAngleTruevsReco->Draw("COLZ");
133  c4->SaveAs("AngleMuons_TruevsReco.pdf");
134 
135  // TCanvas *c5 = new TCanvas("c5", "Pull Neutrons", 800, 600);
136  // hPullNeutrons->GetXaxis()->SetTitle( "(E_{true} - E_{reco}) / E_{true}" );
137  // hPullNeutrons->GetYaxis()->SetTitle( "Normalized Events" );
138  // hPullNeutrons->Draw("hist");
139  // c5->SaveAs("PullEnergyNeutrons.pdf");
140 
141 }
intermediate_table::iterator iterator
static QCString result
void checkmuons()
Definition: checkmuons.cpp:1
string dir
Definition: tf_graph.h:23
list files
Definition: languages.py:9
T abs(T value)
const double a
QTextStream & endl(QTextStream &s)