{
//%%%%%%%%%%%%%%   Original data file from the experiment  %%%%%%%%%%%%%%%%%//

  TString fname="56CoFe0808_IPAC_R10.root";
  TString tree_name="TREE";

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//
//%%%%%%%%%%%%%%%%%% Values considered for Different Cuts %%%%%%%%%%%%%%%%%%//
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//

 Double_t B_up1 = 62;              Double_t B_up2 = 68;      // Ext. Mag. field 'UP' limits  

 Double_t GE1_min = 1;            Double_t GE1_max = 8192;  // Sensible energy range for G1
 Double_t GE2_min = 1;            Double_t GE2_max = 8192;  // Sensible energy range for G2
 Double_t GE3_min = 1;            Double_t GE3_max = 8192;  // Sensible energy range for G3
 Double_t GE4_min = 1;            Double_t GE4_max = 8192;  // Sensible energy range for G4

//%%%%%%%%% TDC limits for Random subtraction %%%%%%%%%%%%%%%%%%%//

 Double_t GT12_Pk1 = 4050;           Double_t GT12_Pk2 = 4090;     // Peak limits
 Double_t GT12_L1 = 3800;            Double_t GT12_L2 = 4000;     // Left background limits
 Double_t GT12_R1 = 4150;            Double_t GT12_R2 = 4450;     // Right background limits

 Double_t GT13_Pk1 = 4050;           Double_t GT13_Pk2 = 4110;     // Peak limits
 Double_t GT13_L1 = 3800;            Double_t GT13_L2 = 4000;     // Left background limits
 Double_t GT13_R1 = 4150;            Double_t GT13_R2 = 4450;     // Right background limits

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//
//%%%%%%%%%%%%%%   Total cuts for data sorting  %%%%%%%%%%%%%%%%%%%%%%%%%%%//
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//
//
//%%%%%%%%%% Magnetic Field Direction Cut %%%%%%%%%%%%%%%%%//
 
  TCut b1 = Form("X.B>%g && X.B<%g",B_up1,B_up2);          

//%%%%%%%%%%%%%%%% Total Gamma energy Cut %%%%%%%%%%%%%%%%%%%//

  TCut ra_GE1 = Form("X.GE1>%g && X.GE1<%g",GE1_min,GE1_max);

  TCut ra_GE2 = Form("X.GE2>%g && X.GE2<%g",GE2_min,GE2_max);

  TCut ra_GE3 = Form("X.GE3>%g && X.GE3<%g",GE3_min,GE3_max);

  TCut ra_GE4 = Form("X.GE4>%g && X.GE4<%g",GE4_min,GE4_max);

//%%%% Gamma TDC Cuts for Random Subtraction %%%%%%%//

  TCut t12 = Form("X.GT12>%g && X.GT12<%g",GT12_Pk1,GT12_Pk2);     
  TCut t12L = Form("X.GT12>%g && X.GT12<%g",GT12_L1,GT12_L2);
  TCut t12R = Form("X.GT12>%g && X.GT12<%g",GT12_R1,GT12_R2);

  TCut t13 = Form("X.GT13>%g && X.GT13<%g",GT13_Pk1,GT13_Pk2);
  TCut t13L = Form("X.GT13>%g && X.GT13<%g",GT13_L1,GT13_L2);
  TCut t13R = Form("X.GT13>%g && X.GT13<%g",GT13_R1,GT13_R2);

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//
//%%%%%%%%%%%%%  Individaul cuts for Gamma Detectros %%%%%%%%%%%%%%%%%%%%%%//
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//

//%%%% Field UP %%%%%%%//

 TCut g1t12_up = ra_GE1 && t12 && b1;           TCut g2t12_up = ra_GE2 && t12 && b1;   //
 TCut g1t12L_up = ra_GE1 && t12L && b1;         TCut g2t12L_up = ra_GE2 && t12L && b1; //
 TCut g1t12R_up = ra_GE1 && t12R && b1;         TCut g2t12R_up = ra_GE2 && t12R && b1;  //

 TCut g1t13_up = ra_GE1 && t13 && b1;           TCut g3t13_up = ra_GE3 && t13 && b1;   //
 TCut g1t13L_up = ra_GE1 && t13L && b1;         TCut g3t13L_up = ra_GE3 && t13L && b1; //
 TCut g1t13R_up = ra_GE1 && t13R && b1;         TCut g3t13R_up = ra_GE3 && t13R && b1;  //

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//
//%%%%%%%%%%%%%%   Normalization Factors for Random Subtaction %%%%%%%%%%%%%%%%//
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//

    //{Normalization Factor = (peak width in TDC)/(background width in TDC)}

//%%%% Top & Bottom %%%%%%%//

 Double_t nr12 = ((GT12_Pk2 - GT12_Pk1 +1) / ((GT12_L2 - GT12_L1 +1) + (GT12_R2 - GT12_R1 +1)));

 Double_t nr13 = ((GT13_Pk2 - GT13_Pk1 +1) / ((GT13_L2 - GT13_L1 +1) + (GT13_R2 - GT13_R1 +1)));

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//
//%%%%%%%%%%%%%%%%  Sorted Matrices for Detector # 1,2,3,4  %%%%%%%%%%%%%%%%//
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//

  TString f12U    = "Umat_G12_56CoFe_IPAC_r10.root";     TString h12U    = "2UG12";
  TString f12U_TC = "Umat_G12_56CoFe_IPAC_TC_r10.root";  TString h12U_TC = "2UG12_TC";
  TString f12U_C  = "Umat_G12_56CoFe_IPAC_C_r10.root";   TString h12U_C  = "2UG12_C";

  TString f13U    = "Umat_G13_56CoFe_IPAC_r10.root";     TString h13U    = "2UG13";
  TString f13U_TC = "Umat_G13_56CoFe_IPAC_TC_r10.root";  TString h13U_TC = "2UG13_TC";
  TString f13U_C  = "Umat_G13_56CoFe_IPAC_C_r10.root";   TString h13U_C  = "2UG13_C";

//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&//
//&&&&&&&&&&&&&&&&&&&&  FINAL DATA SORTING  &&&&&&&&&&&&&&&&&&&&&&&&&&&&&//
//&&&&&&&&&&&&&&&&&&&&&&&  Field =  UP  &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&//

//%%%%%%%%%%%%%%%%%%%%%  Detector # G12 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//
//
//%%%%%% True + Chance(UP) %%%%%%%//
//
 TFile* f=TFile::Open( fname.Data(), "READ" );
 TTree *T = (TTree*)f.Get("tree_name");
//
 TCanvas *C1 = new TCanvas("C1","Fitting canvas",10,10,1400,1200);
 C1->Divide(2,2);
//
 C1->cd(1);
 TH2F *h12 = new TH2F("h12","",8192,-0.5,8191.5,8192,-0.5,8191.5);
// 
 TREE->Draw("X.GE1:X.GE2>>h12",g1t12_up||g2t12_up);
//
 h12->GetXaxis()->SetTitle("EG1");
 h12->GetYaxis()->SetTitle("EG3");
//
 TFile* f_out=TFile::Open( f12U_TC.Data(), "RECREATE" );
 h12->SetName(h12U_TC.Data());
 f_out->cd();
 h12->Write();
 f_out->Close();
//
//%%%%%% Chance(UP) %%%%%%%//
//
 TFile* f=TFile::Open( fname.Data(), "READ" );
 TTree *T = (TTree*)f.Get("tree_name");
//
 C1->cd(2);
 TH2F *hL12 = new TH2F("hL12","",8192,-0.5,8191.5,8192,-0.5,8191.5);
//
 TREE->Draw("X.GE1:X.GE2>>hL12",g1t12L_up||g2t12L_up);
//
 hL12->GetXaxis()->SetTitle("GE1");
 hL12->GetYaxis()->SetTitle("GE2");
//
//
 TFile* f=TFile::Open( fname.Data(), "READ" );
 TTree *T = (TTree*)f.Get("tree_name");
//
 C1->cd(3);
 TH2F *hR12 = new TH2F("hR12","",8192,-0.5,8191.5,8192,-0.5,8191.5);
//
 TREE->Draw("X.GE1:X.GE2>>hR12",g1t12R_up||g2t12R_up);
//
// hR->Draw("LEGO");
 hR12->GetXaxis()->SetTitle("GE1");
 hR12->GetYaxis()->SetTitle("GE2");
//
// total random background
//
 C1->cd(4);
 TH2F *h12c = (TH2F*)hL12->Clone("");
 h12c->Add(hR12);
 h12c->Draw();
 h12c->GetXaxis()->SetTitle("GE1");
 h12c->GetYaxis()->SetTitle("GE2");

//
 TFile* f_out=TFile::Open( f12U_C.Data(), "RECREATE" );
 h12c->SetName(h12U_C.Data());
 f_out->cd();
 h12c->Write();
 f_out->Close();
//
//%%%%%% Random Subtraction (UP) %%%%%%%//
//
//  C1->cd(5);
 TH2F *h_out = (TH2F*)h12->Clone("");
 h_out ->SetDirectory(0);
 h_out->Add(h12c,-1.0*nr12);
// h_out->Draw();
//
 TFile* f_out=TFile::Open(f12U.Data(), "RECREATE" );
 h_out->SetName(h12U.Data());
 f_out->cd();
 h_out->Write();
 f_out->Close();
 f->Close();
//
//
//%%%%%%%%%%%%%%%%%%%%%  Detector # G13  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//
//
//%%%%%% True + Chance(UP) %%%%%%%//
//
 TFile* f=TFile::Open( fname.Data(), "READ" );
 TTree *T = (TTree*)f.Get("tree_name");
//
 TCanvas *C2 = new TCanvas("C2","Fitting canvas",10,10,1400,1200);
 C2->Divide(2,2);
//
 C2->cd(1);
 TH2F *h13 = new TH2F("h13","",8192,-0.5,8191.5,8192,-0.5,8191.5);
// 
 TREE->Draw("X.GE1:X.GE3>>h13",g1t13_up||g3t13_up);
//
 h13->GetXaxis()->SetTitle("EG1");
 h13->GetYaxis()->SetTitle("EG3");
//
 TFile* f_out=TFile::Open( f13U_TC.Data(), "RECREATE" );
 h13->SetName(h13U_TC.Data());
 f_out->cd();
 h13->Write();
 f_out->Close();
//
//%%%%%% Chance(UP) %%%%%%%//
//
 TFile* f=TFile::Open( fname.Data(), "READ" );
 TTree *T = (TTree*)f.Get("tree_name");
//
 C2->cd(2);
 TH2F *hL13 = new TH2F("hL13","",8192,-0.5,8191.5,8192,-0.5,8191.5);
//
 TREE->Draw("X.GE1:X.GE3>>hL13",g1t13L_up||g3t13L_up);
//
 hL13->GetXaxis()->SetTitle("GE1");
 hL13->GetYaxis()->SetTitle("GE3");
//
//
 TFile* f=TFile::Open( fname.Data(), "READ" );
 TTree *T = (TTree*)f.Get("tree_name");
//
 C2->cd(3);
 TH2F *hR13 = new TH2F("hR13","",8192,-0.5,8191.5,8192,-0.5,8191.5);
//
 TREE->Draw("X.GE1:X.GE3>>hR13",g1t13R_up||g3t13R_up);
//
// hR->Draw("LEGO");
 hR13->GetXaxis()->SetTitle("GE1");
 hR13->GetYaxis()->SetTitle("GE3");
//
// total random background
//
 C2->cd(4);
 TH2F *h13c = (TH2F*)hL13->Clone("");
 h13c->Add(hR13);
 h13c->Draw();
 h13c->GetXaxis()->SetTitle("GE1");
 h13c->GetYaxis()->SetTitle("GE3");
//
 TFile* f_out=TFile::Open( f13U_C.Data(), "RECREATE" );
 h13c->SetName(h13U_C.Data());
 f_out->cd();
 h13c->Write();
 f_out->Close();
//
//%%%%%% Random Subtraction (UP) %%%%%%%//
//
//  C2->cd(5);
 TH2F *h_out = (TH2F*)h13->Clone("");
 h_out ->SetDirectory(0);
 h_out->Add(h13c,-1.0*nr13);
// h_out->Draw();
//
 TFile* f_out=TFile::Open(f13U.Data(), "RECREATE" );
 h_out->SetName(h13U.Data());
 f_out->cd();
 h_out->Write();
 f_out->Close();
 f->Close();
//
}