24 #include <TDatabasePDG.h>
56 int main(
int argc,
char* argv[]){
74 cerr <<
"Detector sepcified as " << qapars.
detstring
75 <<
" not recognized or empty." << endl;
83 if ( buildresult !=0 ){
84 cerr <<
"Failed to build a tree from " << qapars.
txtfilename << endl;
90 TString smearedname = rootname;
91 smearedname.ReplaceAll (
".root",
".smeared.root" );
96 SmearTree( detector, rootname.Data(), smearedname.Data());
101 TChain* inTree =
new TChain(
"EICTree");
102 inTree->Add(rootname);
103 inTree->AddFriend(
"Smeared",smearedname);
108 inTree->SetBranchAddress(
"event",&inEvent);
109 inTree->SetBranchAddress(
"eventS",&inEventS);
122 map<int,pidqacollection> qabook;
124 if ( qapars.
pids.size() == 0 ) qapars.
pids = { 11, 211, 321, 2212, 2112 };
130 for(
long iEvent=0; iEvent<inTree->GetEntries(); iEvent++){
133 if(inTree->GetEntry(iEvent) <=0)
break;
134 if(iEvent%10000 == 0) cout <<
"Event " << iEvent << endl;
146 for(
int j=0; j<inEventS->GetNTracks(); j++){
151 const Particle* inParticle = inEvent->GetTrack(j);
154 if ( inParticle->
GetStatus() != 1 )
continue;
157 if(inParticleS == NULL)
continue;
172 qapars.usedevents = inTree->GetEntries();
173 PlotQA( qapars, eventqa, qabook );
185 if ( eventqa.
Q2_NM && inEventS->
GetQ2()>0) eventqa.
Q2_NM->Fill ( std::log10(inEvent->
GetQ2()), std::log10(inEventS->
GetQ2()));
197 if ( eventqa.
y_NM && inEventS->
GetY()>0) eventqa.
y_NM->Fill ( inEvent->
GetY(), inEventS->
GetY());
209 if ( eventqa.
x_NM && inEventS->
GetX()>0) eventqa.
x_NM->Fill ( std::log10(inEvent->
GetX()), std::log10(inEventS->
GetX()));
229 for (
auto& pidcoll : qabook ){
230 auto&
pid = pidcoll.first;
231 auto& coll = pidcoll.second;
235 auto delP = (inParticle->
GetP() - inParticleS->
GetP()) / inParticle->
GetP();
236 coll.DelP_th->Fill(inParticle->
GetTheta(), delP);
237 coll.DelP_eta->Fill(inParticle->
GetEta(), delP);
241 auto delE = (inParticle->
GetE() - inParticleS->
GetE()) / inParticle->
GetE();
242 coll.DelE_E->Fill(inParticle->
GetE(), delE);
243 coll.DelE_th->Fill(inParticle->
GetTheta(), delE);
244 coll.DelE_eta->Fill(inParticle->
GetEta(), delE);
249 coll.dEta_p->Fill(inParticle->
GetP(), inParticle->
GetEta() - inParticleS->
GetEta());
253 coll.dPhi_p->Fill(inParticle->
GetP(), inParticle->
GetPhi() - inParticleS->
GetPhi());
262 vector<string> arguments(argv + 1, argv + argc);
266 for (
auto parg = arguments.begin() ; parg!=arguments.end() ; ++parg){
271 }
else if ( arg ==
"-o" ){
272 if (++parg == arguments.end() ){ argsokay=
false;
break; }
274 }
else if ( arg ==
"-i" ){
275 if (++parg ==arguments.end() ){ argsokay=
false;
break; }
277 }
else if ( arg ==
"-N" ){
278 if (++parg==arguments.end() ){ argsokay=
false;
break; }
279 qapars.
nevents=std::stoi(parg->data());
280 }
else if ( arg ==
"-addpid" ){
281 if ( ++parg == arguments.end() ){ argsokay=
false;
break; }
282 qapars.
pids.push_back(std::stoi(parg->data()));
283 }
else if ( arg ==
"-det" ){
284 if (++parg == arguments.end() ){ argsokay=
false;
break; }
287 if ( TString(qapars.
detstring).Contains(
"MATRIX") && TString(qapars.
detstring).Contains(
"FF")){
288 if (++parg == arguments.end() ){ argsokay=
false;
break; }
291 if ( TString(qapars.
detstring).Contains(
"CORE") && !TString(qapars.
detstring).Contains(
"B")){
292 if (++parg == arguments.end() ){ argsokay=
false;
break; }
300 }
catch (
const std::exception&
e){
301 cerr <<
"Caught exception during argument parsing: "
307 cerr <<
"usage: " << argv[0] << endl
308 <<
" [-i txtfilename] (Lund-style file)" << endl
309 <<
" [-o OutFileBase] (extension will be added)" << endl
310 <<
" [-N Nevents] (<0 for all)" << endl
311 <<
" [-addpid pid] (can be called multiple times)" << endl
312 <<
" [-det detstring] matrix02B3, matrix02B15, matrix01, matrixff [beam_mom_nn], handbook, perfect, beast, ephenix, zeus, jleic (capitalization does not matter.)" << endl
314 throw std::runtime_error(
"Not a valid list of options");
328 if ( rootname.EndsWith(
".gz", TString::kIgnoreCase) ||
329 rootname.EndsWith(
".zip", TString::kIgnoreCase) )
330 rootname.Replace(rootname.Last(
'.'), rootname.Length(),
"");
333 if (rootname.Last(
'.') > -1) {
334 rootname.Replace(rootname.Last(
'.'), rootname.Length(),
"");
344 rootname = gSystem->BaseName( rootname );
345 rootname.Prepend( qapars.
outpath );
346 cout <<
" ======================= " << endl;
347 cout <<
" Transforming input file " << endl
349 <<
" into root file " << endl
355 gStyle->SetHistLineColor(kRed);
375 float thmax = TMath::Pi();
386 float detamin = -0.1;
391 float phimax = TMath::TwoPi();
394 float dphimin = -0.1;
398 for (
auto pid : qapars.
pids ){
402 qabook[
pid].DelE_E =
new TH2D( s,s+
";E;#Delta E/E", ebins, emin, emax, debins, demin, demax);
405 qabook[
pid].dPhi_p =
new TH2D( s,s+
";p;#Delta#phi", pbins, pmin, pmax, dphibins, dphimin, dphimax );
408 qabook[
pid].DelP_th =
new TH2D( s,s+
";#theta;#Delta p/p", thbins, thmin, thmax, dpbins, dpmin, dpmax);
411 qabook[
pid].DelE_th =
new TH2D( s,s+
";#theta;#Delta E/E", thbins, thmin, thmax, debins, demin, demax);
414 qabook[
pid].dTh_p =
new TH2D( s,s+
";p;#Delta#theta", pbins, pmin, pmax, dthbins, dthmin, dthmax );
417 qabook[
pid].DelP_eta =
new TH2D( s,s+
";#eta;#Delta p/p", etabins, etamin, etamax, dpbins, dpmin, dpmax);
420 qabook[
pid].DelE_eta =
new TH2D( s,s+
";#eta;#Delta E/E", etabins, etamin, etamax, debins, demin, demax);
423 qabook[
pid].dEta_p =
new TH2D( s,s+
";p;#Delta#eta", pbins, pmin, pmax, detabins, detamin, detamax );
440 eventqa.
Q2_NM =
new TH2D( s,s+
";log Q^{2};log Q^{2}_{NM}", logQ2bins, std::log10(Q2min), std::log10(Q2max), logQ2bins, std::log10(Q2min), std::log10(Q2max));
443 eventqa.
Q2_JB =
new TH2D( s,s+
";log Q^{2};log Q^{2}_{JB}", logQ2bins, std::log10(Q2min), std::log10(Q2max), logQ2bins, std::log10(Q2min), std::log10(Q2max));
446 eventqa.
Q2_DA =
new TH2D( s,s+
";log Q^{2};log Q^{2}_{DA}", logQ2bins, std::log10(Q2min), std::log10(Q2max), logQ2bins, std::log10(Q2min), std::log10(Q2max));
450 eventqa.
delQ2_NM =
new TH2D( s,s+
";log Q^{2};(Q^{2}_{NM}-Q^{2})/Q2", logQ2bins, std::log10(Q2min), std::log10(Q2max), delQ2bins, 0, 1.2);
453 eventqa.
delQ2_JB =
new TH2D( s,s+
";log Q^{2};(Q^{2}_{JB}-Q^{2})/Q2", logQ2bins, std::log10(Q2min), std::log10(Q2max), delQ2bins, 0, 1.2);
456 eventqa.
delQ2_DA =
new TH2D( s,s+
";log Q^{2};(Q^{2}_{DA}-Q^{2})/Q2", logQ2bins, std::log10(Q2min), std::log10(Q2max), delQ2bins, 0, 1.2);
464 eventqa.
y_NM =
new TH2D( s,s+
";y;y_{NM}", ybins, ymin, ymax, ybins, ymin, ymax);
467 eventqa.
y_JB =
new TH2D( s,s+
";y;y_{JB}", ybins, ymin, ymax, ybins, ymin, ymax);
470 eventqa.
y_DA =
new TH2D( s,s+
";y;y_{DA}", ybins, ymin, ymax, ybins, ymin, ymax);
474 eventqa.
dely_NM =
new TH2D( s,s+
";y;(y_{NM}-y)/y", ybins, ymin, ymax, delybins, 0, 1.2);
477 eventqa.
dely_JB =
new TH2D( s,s+
";y;(y_{JB}-y)/y", ybins, ymin, ymax, delybins, 0, 1.2);
480 eventqa.
dely_DA =
new TH2D( s,s+
";y;(y_{DA}-y)/y", ybins, ymin, ymax, delybins, 0, 1.2);
509 eventqa.
x_NM =
new TH2D( s,s+
";log x;log x_{NM}", logxbins, std::log10(xmin), std::log10(xmax), logxbins, std::log10(xmin), std::log10(xmax));
512 eventqa.
x_JB =
new TH2D( s,s+
";log x;log x_{JB}", logxbins, std::log10(xmin), std::log10(xmax), logxbins, std::log10(xmin), std::log10(xmax));
515 eventqa.
x_DA =
new TH2D( s,s+
";log x;log x_{DA}", logxbins, std::log10(xmin), std::log10(xmax), logxbins, std::log10(xmin), std::log10(xmax));
519 eventqa.
delx_NM =
new TH2D( s,s+
";log x;(x_{NM}-x)/x", logxbins, std::log10(xmin), std::log10(xmax), delxbins, 0, 1.2);
522 eventqa.
delx_JB =
new TH2D( s,s+
";log x;(x_{JB}-x)/x", logxbins, std::log10(xmin), std::log10(xmax), delxbins, 0, 1.2);
525 eventqa.
delx_DA =
new TH2D( s,s+
";log x;(x_{DA}-x)/x", logxbins, std::log10(xmin), std::log10(xmax), delxbins, 0, 1.2);
533 gStyle->SetStatX(0.25);
534 gStyle->SetStatW(0.15);
535 gStyle->SetStatY(0.9);
536 gStyle->SetStatH(0.15);
547 gErrorIgnoreLevel = kWarning;
557 if ( eventqa.
y_NM ) {
558 eventqa.
y_NM->Draw(
"colz");
562 if ( eventqa.
x_NM ) {
563 eventqa.
x_NM->Draw(
"colz");
567 if ( eventqa.
Q2_NM ) {
568 eventqa.
Q2_NM->Draw(
"colz");
573 if ( eventqa.
y_DA ) {
574 eventqa.
y_DA->Draw(
"colz");
578 if ( eventqa.
x_DA ) {
579 eventqa.
x_DA->Draw(
"colz");
583 if ( eventqa.
Q2_DA ) {
584 eventqa.
Q2_DA->Draw(
"colz");
589 if ( eventqa.
y_JB ) {
590 eventqa.
y_JB->Draw(
"colz");
594 if ( eventqa.
x_JB ) {
595 eventqa.
x_JB->Draw(
"colz");
599 if ( eventqa.
Q2_JB ) {
600 eventqa.
Q2_JB->Draw(
"colz");
657 gStyle->SetStatX(0.55);
658 for (
auto& pidcoll : qabook ){
659 auto&
pid = pidcoll.first;
660 auto& coll = pidcoll.second;
664 coll.DelP_th->Draw(
"colz");
667 coll.DelP_eta->Draw(
"colz");
668 coll.DelP_eta->ProfileX(
"_px",1,-1,
"s")->Draw(
"same");
671 gStyle->SetStatX(0.25);
672 coll.DelE_th->Draw(
"colz");
675 coll.DelE_eta->Draw(
"colz");
676 coll.DelE_eta->ProfileX(
"_px",1,-1,
"s")->Draw(
"same");
678 gStyle->SetStatX(0.55);
680 coll.DelE_E->Draw(
"colz");
681 coll.DelE_E->ProfileX(
"_px",1,-1,
"s")->Draw(
"same");
684 coll.dTh_p->Draw(
"colz");
687 coll.dEta_p->Draw(
"colz");
688 coll.dEta_p->ProfileX(
"_px",1,-1,
"s")->Draw(
"same");
691 coll.dPhi_p->Draw(
"colz");
692 coll.dPhi_p->ProfileX(
"_px",1,-1,
"s")->Draw(
"same");
697 gErrorIgnoreLevel = kInfo;