//macro to set momentum cuts
//this cut should be applied after back-propagation
//date: 22/02/2014
//author: A.Karavdina

//void MomCut(TString path = "/panda/myResults/G4_ph_lists/BOX_2_12mrad_G4_FTF_BERT_"){
  void MomCut(TString path = "/panda/myResults/DPMoutputDec2013/FTF/XYCut/"){
  // const int npmom=5;
 //  //  const int npmom=1;
 //  double Pbeam[npmom]={1.5,4.06,8.9,11.91,15};
 //  TString Pbeam_str[npmom] = {"mom_1_5","mom_4_06","mom_8_9","mom_11_91","mom_15"};
 // TString Pbeam_str_name[npmom] = {"P_{beam} = 1.5 GeV/c","P_{beam} = 4.06 GeV/c","P_{beam} = 8.9 GeV/c","P_{beam} = 11.91 GeV/c","P_{beam} = 15 GeV/c"};
 //  // double Thres[npmom]={685,228,103,86,77};//resolution in murad
 //  double Thres[npmom]={601,212,100,83,74};//resolution in murad [NEW official]

 const int npmom=2;
  double Pbeam[npmom]={1.5,15};
  TString Pbeam_str[npmom] = {"mom_1_5","mom_15"};
 TString Pbeam_str_name[npmom] = {"P_{beam} = 1.5 GeV/c","P_{beam} = 15 GeV/c"};
  double Thres[npmom]={601,74};//resolution in murad [NEW official]

  TH2D *hthmom2[npmom];
  TH2D *hthmom3[npmom];
  TH2D *hthmom4[npmom];
  TH2D *hthmom5[npmom];
  for(int qb=0;qb<npmom;qb++){//loop over Pbeam
  // for(int qb=0;qb<1;qb++){//loop over Pbeam
  //  TString path_XYCuts = path+"XYCut_3sigma/";
    //    TString fname_XYCuts = path_XYCuts;
    TString fname_XYCuts = path;
    fname_XYCuts+=Pbeam_str[qb];
    //    fname_XYCuts+="/Lumi_TrksQA_SUM.root";
    fname_XYCuts+="/Lumi_TrksQA_SUM_el.root";
    TFile *fin_XYCuts = new TFile(fname_XYCuts,"READ");
    if(fin_XYCuts->IsZombie()) cout<<"File: "<<fname_XYCuts<<" =("<<endl;
    TTree *tTrks = (TTree *)fin_XYCuts->Get("cbmsim"); //Bkg data
    TString fhistname = "hthdmom_";
    fhistname+=Pbeam_str[qb];
    TString fhistname1 = fhistname;
    fhistname1 +="_1sig";
    TString fhistname2 = fhistname;
    fhistname2 +="_2sig";
    TString fhistname3 = fhistname;
    fhistname3 +="_3sig";
    TString fhistname4 = fhistname;
    fhistname4 +="_4sig";
    TString fhistname5 = fhistname;
    fhistname5 +="_5sig";
    TString vars = "1e3*(LMDTrackQ.fThetamc):abs(";
    vars +=Pbeam[qb];
    vars +="-LMDTrackQ.fMomrec)/";
    vars +=Pbeam[qb];

    hthmom2[qb] = new TH2D(fhistname2,";|P_{beam}-P_{REC}|/P_{beam};#theta_{MC}, mrad",1e5,0,1e0,1e2,0,15);
    // tTrks->Project(fhistname2,"1e3*(LMDTrackQ.fThetamc):abs(LMDTrackQ.fMommc-LMDTrackQ.fMomrec)/LMDTrackQ.fMommc","LMDTrackQ.fSecondary<0 && LMDTrackQ.fTrkRecStatus==0 && 1e3*(LMDTrackQ.fThetamc-LMDTrackQ.fThetarec)/0.6<2.");
    tTrks->Project(fhistname2,vars,"LMDTrackQ.fSecondary<0 && LMDTrackQ.fTrkRecStatus==0 && 1e3*(LMDTrackQ.fThetamc-LMDTrackQ.fThetarec)/0.6<2.");
    hthmom3[qb] = new TH2D(fhistname3,";|P_{MC}-P_{REC}|/P_{MC};#theta_{MC}, mrad",1e5,0,1e0,1e2,0,15);
    tTrks->Project(fhistname3,vars,"LMDTrackQ.fSecondary<0 && LMDTrackQ.fTrkRecStatus==0 && 1e3*(LMDTrackQ.fThetamc-LMDTrackQ.fThetarec)/0.6>3.");
    hthmom4[qb] = new TH2D(fhistname4,";|P_{MC}-P_{REC}|/P_{MC};#theta_{MC}, mrad",1e5,0,1e0,1e2,0,15);
    tTrks->Project(fhistname4,vars,"LMDTrackQ.fSecondary<0 && LMDTrackQ.fTrkRecStatus==0 && 1e3*(LMDTrackQ.fThetamc-LMDTrackQ.fThetarec)/0.6>4.");
    hthmom5[qb] = new TH2D(fhistname5,";|P_{MC}-P_{REC}|/P_{MC};#theta_{MC}, mrad",1e5,0,1e0,1e2,0,15);
    tTrks->Project(fhistname5,vars,"LMDTrackQ.fSecondary<0 && LMDTrackQ.fTrkRecStatus==0 && 1e3*(LMDTrackQ.fThetamc-LMDTrackQ.fThetarec)/0.6>5.");
    TCanvas *c1  = new TCanvas();
    // hthmom1[qb]->SetMarkerStyle(20);
    // hthmom1[qb]->SetMarkerColor(2);
    // hthmom1[qb]->Draw();
    hthmom2[qb]->SetMarkerStyle(25);
    hthmom2[qb]->SetMarkerColor(4);
    hthmom2[qb]->Draw("BOX");
    hthmom3[qb]->SetMarkerStyle(31);
    hthmom3[qb]->SetMarkerColor(3);
    hthmom3[qb]->Draw("same");
    hthmom4[qb]->SetMarkerStyle(2);
    hthmom4[qb]->SetMarkerColor(1);
    hthmom4[qb]->Draw("same");
    hthmom5[qb]->SetMarkerStyle(33);
    hthmom5[qb]->SetMarkerColor(5);
    hthmom5[qb]->Draw("same");
    TLegend *leg2 = new TLegend(0.75,0.75,0.99,0.99);
    //    leg2->SetHeader("BOX: 2 mrad < #theta_{MC} < 12 mrad");
    leg2->SetHeader("DPM: elastic only");
    leg2->SetTextAlign(12);
    leg2->SetTextSize(0.03);
    leg2->SetFillColor(0);
    //    leg2->AddEntry(hthmom1[qb],"|#theta_{MC}-#theta_{REC}|>1#sigma_{#theta}","p");
    leg2->AddEntry((TObject*)0,Pbeam_str_name[qb]," ");
    leg2->AddEntry(hthmom2[qb],"|#theta_{MC}-#theta_{REC}|<2#sigma_{#theta}","p");
    leg2->AddEntry(hthmom3[qb],"|#theta_{MC}-#theta_{REC}|>3#sigma_{#theta}","p");
    leg2->AddEntry(hthmom4[qb],"|#theta_{MC}-#theta_{REC}|>4#sigma_{#theta}","p");
    leg2->AddEntry(hthmom5[qb],"|#theta_{MC}-#theta_{REC}|>5#sigma_{#theta}","p");
    leg2->Draw();
    TString filename = fhistname;
    filename +=".root";
    c1.SaveAs(filename);
    //  fin_XYCuts->Close();
  }
}