Adding macro for charge

Projects/ana_e850/macro/Charge/FillHistoCharge.C

+TChain* chain;
+TCutG* cut10Be[8];
+
+/////////////////////////////////////////////////////////////////
+void LoadCut(string part){
+  string path = "../Mass/cut/"+ part + "/cut10Be.root";
+  TFile* fcut = new TFile(path.c_str());
+  for(int i=0; i<8; i++){
+   cut10Be[i] = (TCutG*)fcut->Get(Form("cut10Be_det%i",i+1));
+  }
+}
+
+/////////////////////////////////////////////////////////////////
+void FillHistoCharge(string part="part1"){
+
+  LoadCut(part);
+
+  chain = new TChain("PhysicsTree");
+  if(part=="part1"){
+    //chain->Add("../../root/analysis/run_28.root");
+    //chain->Add("../../root/analysis/run_29.root");
+    //chain->Add("../../root/analysis/run_30.root");
+    //chain->Add("../../root/analysis/run_32.root");
+    //chain->Add("../../root/analysis/run_36.root");
+    //chain->Add("../../root/analysis/run_37.root");
+    //chain->Add("../../root/analysis/run_39.root");
+    //chain->Add("../../root/analysis/run_41.root");
+    //chain->Add("../../root/analysis/run_42.root");
+    //chain->Add("../../root/analysis/run_43.root");
+    //chain->Add("../../root/analysis/run_44.root");
+    //chain->Add("../../root/analysis/run_45.root");
+    //chain->Add("../../root/analysis/run_46.root");
+    chain->Add("../../root/analysis/run_48.root");
+    //chain->Add("../../root/analysis/run_49.root");
+    //chain->Add("../../root/analysis/run_50.root");
+    //chain->Add("../../root/analysis/run_51.root");
+    //chain->Add("../../root/analysis/run_52.root");
+    //chain->Add("../../root/analysis/run_53.root");
+  }
+  else if(part=="part2"){
+    chain->Add("../../root/analysis/run_55.root");
+    chain->Add("../../root/analysis/run_56.root");
+    chain->Add("../../root/analysis/run_57.root");
+    chain->Add("../../root/analysis/run_58.root");
+    chain->Add("../../root/analysis/run_59.root");
+    chain->Add("../../root/analysis/run_60.root");
+    chain->Add("../../root/analysis/run_61.root");
+    chain->Add("../../root/analysis/run_62.root");
+    chain->Add("../../root/analysis/run_63.root");
+    chain->Add("../../root/analysis/run_64.root");
+    chain->Add("../../root/analysis/run_65.root");
+    chain->Add("../../root/analysis/run_66.root");
+    chain->Add("../../root/analysis/run_67.root");
+    chain->Add("../../root/analysis/run_68.root");
+    chain->Add("../../root/analysis/run_69.root");
+  }
+  double DeltaEcorr;
+  double Elab;
+  double Mass;
+  double Ex240Pu;
+  double VAMOS_TS_hour;
+  int Telescope;
+  int m_2alpha;
+  int FPMW_Section;
+  double FF_Z;
+
+  chain->SetBranchStatus("DeltaEcorr","true");
+  chain->SetBranchAddress("DeltaEcorr",&DeltaEcorr);
+  chain->SetBranchStatus("Elab","true");
+  chain->SetBranchAddress("Elab",&Elab);
+  chain->SetBranchStatus("FF_Mass13","true");
+  chain->SetBranchAddress("FF_Mass13",&Mass);
+  chain->SetBranchStatus("Ex240Pu","true");
+  chain->SetBranchAddress("Ex240Pu",&Ex240Pu);
+  chain->SetBranchStatus("VAMOS_TS_hour","true");
+  chain->SetBranchAddress("VAMOS_TS_hour",&VAMOS_TS_hour);
+  chain->SetBranchStatus("Telescope","true");
+  chain->SetBranchAddress("Telescope",&Telescope);
+  chain->SetBranchStatus("m_2alpha","true");
+  chain->SetBranchAddress("m_2alpha",&m_2alpha);
+  chain->SetBranchStatus("FPMW_Section","true");
+  chain->SetBranchAddress("FPMW_Section",&FPMW_Section);
+  chain->SetBranchStatus("FF_Z","true");
+  chain->SetBranchAddress("FF_Z",&FF_Z);
+
+
+
+  TH1F *hcharge[11];
+  for(int i=0; i<11; i++){
+    TString histo_name = Form("hcharge_%iMeV",i+6);
+    hcharge[i] = new TH1F(histo_name, histo_name,1000,20,80); 
+  }
+  TH1F* hChargeAll = new TH1F("hcharge_all","hcharge_all",1000,20,80);
+  TH2D* hChargeVsEx = new TH2D("hcharge_vs_ex","hcharge_vs_ex",20,0,20,1000,20,80);
+
+  int nentries = chain->GetEntries();
+  cout << "**** Number of entries to be treated: " << nentries << " ****" << endl; 
+  for(int i=0; i<nentries; i++){
+    chain->GetEntry(i);
+    if(i%100000==0){
+      cout << "\033[34m\r Processing tree..." << (double)i/nentries*100 << "\% done" << flush;
+    }
+
+    if(VAMOS_TS_hour>0){
+
+      hChargeAll->Fill(FF_Z);
+      if(Telescope>0) {
+
+        if(FPMW_Section>4 && cut10Be[Telescope-1]->IsInside(Elab,DeltaEcorr)){
+          hChargeVsEx->Fill(Ex240Pu,FF_Z);
+          if(Ex240Pu>5.5 && Ex240Pu<6.5) hcharge[0]->Fill(FF_Z);
+          else if(Ex240Pu>6.5 && Ex240Pu<7.5) hcharge[1]->Fill(FF_Z);
+          else if(Ex240Pu>7.5 && Ex240Pu<8.5) hcharge[2]->Fill(FF_Z);
+          else if(Ex240Pu>8.5 && Ex240Pu<9.5) hcharge[3]->Fill(FF_Z);
+          else if(Ex240Pu>9.5 && Ex240Pu<10.5) hcharge[4]->Fill(FF_Z);
+          else if(Ex240Pu>10.5 && Ex240Pu<11.5) hcharge[5]->Fill(FF_Z);
+          else if(Ex240Pu>11.5 && Ex240Pu<12.5) hcharge[6]->Fill(FF_Z);
+          else if(Ex240Pu>12.5 && Ex240Pu<13.5) hcharge[7]->Fill(FF_Z);
+          else if(Ex240Pu>13.5 && Ex240Pu<14.5) hcharge[8]->Fill(FF_Z);
+          else if(Ex240Pu>14.5 && Ex240Pu<15.5) hcharge[9]->Fill(FF_Z);
+          else if(Ex240Pu>15.5 && Ex240Pu<16.5) hcharge[10]->Fill(FF_Z);
+        }
+      }
+    }
+  }
+
+  TString filename = Form("histo_charge_%s.root",part.c_str());
+  TFile *ofile = new TFile(filename,"recreate");
+  for(int i=0; i<11; i++){
+    hcharge[i]->Write();
+  }
+  hChargeAll->Write();
+  hChargeVsEx->Write();
+  ofile->Close();
+
+
+
+}

Projects/ana_e850/macro/Charge/FitCharge.C

+void Fit(TH1F* hcharge, string Energy);
+TFile* ofile;
+int rebin=1;
+
+//////////////////////////////////////////////////////
+void FitCharge(){
+
+  TFile* ifile = new TFile("histo_charge_part1.root","read");
+
+  ofile = new TFile("histo_charge_fitted.root","recreate");
+
+  TH1F* hall = (TH1F*)ifile->FindObjectAny("hcharge_all");
+  Fit(hall,"AllEnergy");
+  /*for(int i=0; i<11; i++){
+    int Ex = i+6;
+    if(i>4) rebin=2;
+    else rebin=1;
+
+    string Energy = to_string(Ex) + "MeV";
+    TString histo_name = Form("hcharge_%iMeV",i+6);
+    TH1F* h1 = (TH1F*)ifile->FindObjectAny(histo_name);
+    Fit(h1,Energy);
+  }*/
+
+  ofile->Close();
+}
+
+//////////////////////////////////////////////////////
+void Fit(TH1F* hcharge,string Energy){
+  TGraphErrors* gsigma = new TGraphErrors();
+  hcharge->Rebin(rebin);
+  hcharge->Draw();
+
+  int Zmin = 32;
+  int Zmax = 65;
+
+
+  int Z = Zmin;
+
+  int NumberOfZ = Zmax - Zmin;
+  Double_t para[3*NumberOfZ];
+
+  TF1* g[NumberOfZ];
+  double Integral[NumberOfZ];
+  double Integral_err[NumberOfZ];
+  double total_integral = 0;
+  double Yield[NumberOfZ];
+
+  for(int i=0; i<NumberOfZ; i++){
+    g[i] = new TF1(Form("g%i",i),"gaus",Z-0.4,Z+0.4);   
+    if(i==0){
+      g[i]->SetParameter(0,100);
+      g[i]->SetParameter(1,Zmin);
+      g[i]->SetParameter(2,0.3);
+    }
+    else if(i>0){
+      g[i]->SetParameter(0,g[i-1]->GetParameter(0));
+      g[i]->SetParameter(1,Z);
+      g[i]->SetParameter(2,g[i-1]->GetParameter(2));
+    }
+ 
+    g[i]->SetParLimits(1,Z-0.2,Z+0.2);
+    g[i]->SetParLimits(2,0.22,0.4);
+
+   
+    hcharge->Fit(g[i],"qBR");
+
+    g[i]->Draw("lsame");
+    g[i]->GetParameters(&para[3*i]);
+
+    Z++;
+  }
+
+
+  TString total_func = "gaus(0)";
+  for(int i=1; i<NumberOfZ; i++){
+    TString gaus_func = Form("gaus(%i)",3*i);
+    total_func += "+" + gaus_func;
+  }
+
+  TF1* total = new TF1("total",total_func,Zmin-0.5, Zmin+NumberOfZ-0.5);
+  total->SetParameters(para);
+  Z = Zmin;
+  for(int i=0; i<NumberOfZ; i++){
+    total->SetParameter(3*i,g[i]->GetParameter(0));
+    total->SetParameter(3*i+1,g[i]->GetParameter(1));
+    total->SetParameter(3*i+2,g[i]->GetParameter(2));
+    total->SetParLimits(3*i+1,Z-0.2,Z+0.2);
+    total->SetParLimits(3*i+2,g[i]->GetParameter(2)*0.85,g[i]->GetParameter(2)*1.01);
+    Z++;
+  }
+  total->SetLineColor(4);
+  total->SetNpx(1000);
+
+  hcharge->Fit(total,"RBq");
+
+  Z = Zmin;
+  //TF1* one_gauss = new TF1("one_gaus","gaus");
+  for(int i=0; i<NumberOfZ; i++){
+      
+    double Amplitude = total->GetParameter(3*i);
+    double mean = total->GetParameter(3*i+1);
+    double sigma = total->GetParameter(3*i+2);
+    double Amplitude_err = total->GetParError(3*i);
+    double mean_err = total->GetParError(3*i+1);
+    double sigma_err = total->GetParError(3*i+2);
+ 
+    /*double Amplitude = g[i]->GetParameter(0);
+    double mean = g[i]->GetParameter(1);
+    double sigma = g[i]->GetParameter(2);
+    double Amplitude_err = g[i]->GetParError(0);
+    double mean_err = g[i]->GetParError(1);
+    double sigma_err = g[i]->GetParError(2);*/
+    
+    //one_gauss->SetParameter(0,Amplitude);
+    //one_gauss->SetParameter(1,mean);
+    //one_gauss->SetParameter(2,sigma);
+    //one_gauss->SetParError(0,Amplitude_err);
+    //one_gauss->SetParError(1,mean_err);
+    //one_gauss->SetParError(2,sigma_err);
+
+    cout << "************************" << endl;
+    cout << "Z=  " << Z << endl;
+    //cout << "Amplitude= " << Amplitude << "+/-" << Amplitude_err << endl;
+    //cout << "Mean= " << mean  << endl;
+    cout << "Sigma= " << sigma << "+/-" << sigma_err << endl;
+
+    gsigma->SetPoint(i,Z,sigma);
+    gsigma->SetPointError(i,0,sigma_err);
+    
+
+    Integral[i] = (12.5/rebin)*Amplitude*sigma*sqrt(2*TMath::Pi());
+    Integral_err[i] = (12.5/rebin)*sqrt(2*TMath::Pi()*(pow(sigma*Amplitude_err,2) + pow(Amplitude*sigma_err,2)));
+    //Integral_err[i] = sqrt(Integral[i]);
+    
+
+    total_integral += Integral[i];
+    Z++;
+  }
+
+  double histo_integral = hcharge->Integral();
+  hcharge->Write();
+  total->Write();
+
+  string output_filename = "Charge_Yield_" + Energy + ".dat";
+  ofstream ofile;
+  ofile.open(output_filename.c_str());
+
+  // Yield calculation //
+  Z = Zmin;
+  for(int i=0; i<NumberOfZ; i++){
+    //Yield[i] = Integral[i]/total_integral*200;
+    Yield[i] = Integral[i]/histo_integral*200;
+    //double yield_err = Integral_err[i]/total_integral*200;
+    double yield_err = Integral_err[i]/histo_integral*200;
+    ofile << Z << " " << Yield[i] << " " << yield_err << endl;
+    Z++;
+  }
+  ofile.close();
+
+  gsigma->SetMarkerStyle(8);
+  gsigma->Write();
+
+  TCanvas* csig = new TCanvas("csig","csig",800,800);
+  csig->cd();
+  gsigma->Draw("ap");
+
+}