PhotonConversion/ana/conversionRate.C

0001
0002 TChain* handleFile(string name, string extension, string treename, unsigned int filecount){
0003   TChain *all = new TChain(treename.c_str());
0004   string temp;
0005   for (int i = 0; i < filecount; ++i)
0006   {
0007
0008     ostringstream s;
0009     s<<i;
0010     temp = name+string(s.str())+extension;
0011     all->Add(temp.c_str());
0012   }
0013   return all;
0014 }
0015
0016 /**makes a TEff holding the hists for a uniform conversion rate
0017  * @param out_file file data is saved to this is used for read and write needs converted and truth pT spectra to be inside these are created by photonEff
0018  * if the pT spectra are not inside the file @param ttree and @param allTree will be used to calculate these spectra
0019  * @param ttree tree for converted data
0020  * @param allTree tree for nonconverted data
0021  * @return NULL if the spectra are not in the existing data and the passes trees are not valid*/
0022 TEfficiency* makepTRes(TFile* out_file,TChain* ttree=NULL,TTree* allTree=NULL){
0023   //check if the spectra are in the file
0024   out_file->ReOpen("READ");
0025   if(out_file->Get("converted_photon_truth_pT")&&out_file->Get("all_photon_truth_pT"))
0026     return new TEfficiency(*(TH1F*)out_file->Get("converted_photon_truth_pT"),*(TH1F*)out_file->Get("all_photon_truth_pT"));
0027   //if they are not in the file check if the trees are NULL
0028   else if(!ttree||!allTree){
0029     return NULL;
0030   }
0031   out_file->ReOpen("UPDATE");
0032   float pT;
0033   float tpT;
0034   float track_pT;
0035   ttree->SetBranchAddress("photon_pT",&pT);
0036   ttree->SetBranchAddress("tphoton_pT",&tpT);
0037
0038   vector<float> *allpT=NULL;
0039   allTree->SetBranchAddress("photon_pT",&allpT);
0040
0041   TH1F *pTeffPlot = new TH1F("#frac{#it{p}^{T}}{#it{p}_{#it{truth}}^{T}}","",40,-2,2);
0042   TH2F *pTefffuncPlot = new TH2F("pT_resolution_to_truthpt","",40,1,35,40,-1.5,1.5);
0043   TH1F *converted_pTspec = new TH1F("converted_photon_truth_pT","",20,5,25);
0044   TH1F *all_pTspec = new TH1F("all_photon_truth_pT","",20,5,25);
0045   //TH1F *trackpTDist = new TH1F("truthpt","",40,0,35);
0046   pTeffPlot->Sumw2();
0047   converted_pTspec->Sumw2();
0048   all_pTspec->Sumw2();
0049   pTefffuncPlot->Sumw2();
0050   //TODO need to turn off other branches
0051   //make the pT spectra fo converted
0052   for (int event = 0; event < ttree->GetEntries(); ++event)
0053   {
0054     ttree->GetEvent(event);
0055     if(pT>0)pTeffPlot->Fill(pT/tpT);
0056     converted_pTspec->Fill(tpT);
0057     if(pT>0)pTefffuncPlot->Fill(tpT,pT/tpT);
0058     //trackpTDist->Fill(track_pT);
0059   }
0060   //make the pT spectra for unconverted
0061   for (int event = 0; event < allTree->GetEntries(); ++event)
0062   {
0063     allTree->GetEvent(event);
0064     for(auto i : *allpT){
0065       all_pTspec->Fill(i);
0066     }
0067   }
0068   //format and save the data
0069   TEfficiency* uni_rate = new TEfficiency(*converted_pTspec,*all_pTspec);
0070   uni_rate->SetName("uni_rate");
0071   uni_rate->Write();
0072   pTeffPlot->Scale(1./ttree->GetEntries(),"width");
0073   pTefffuncPlot->Scale(1./ttree->GetEntries());
0074   TProfile* resProfile = pTefffuncPlot->ProfileX("func_prof",5,30);
0075   resProfile->Write();
0076   //trackpTDist->Scale(1./ttree->GetEntries(),"width");
0077   out_file->Write();
0078   ttree->ResetBranchAddresses();
0079   return uni_rate;
0080 }
0081
0082 unsigned totalMB(string path, string extension, string hardname, string softname){
0083   string name = path+hardname+extension;
0084   TFile *hardFile=new TFile(name.c_str(),"READ");
0085   name=path+softname+extension;
0086   TFile *softFile=new TFile(name.c_str(),"READ");
0087   TTree* noPhoton = (TTree*) softFile->Get("nophotonTree");
0088   unsigned r=noPhoton->GetEntries()*2000000;
0089   noPhoton = (TTree*) hardFile->Get("nophotonTree");
0090   r+=noPhoton->GetEntries()*2000000;
0091   return r;
0092 }
0093
0094 /**@param rate the uniform conversion rate
0095  * @param file file to get the pythia pT spectra from
0096  * @return the pT dependent conversion rate*/
0097 TH1F* calculateRate(TEfficiency* rate,TFile* file,unsigned nMB){
0098   //get the combined pythiaspec from the file then clone it to rate
0099   TH1F* conversion_rate = (TH1F*)((TH1F*) file->Get("combinedpythia"))->Clone("rate");
0100   //make a histogram for the uniform rate
0101   TH1* uni_rate = (TH1F*)rate->GetPassedHistogram()->Clone("uni_rate");
0102   uni_rate->Divide(rate->GetTotalHistogram());
0103   conversion_rate->Multiply(uni_rate);
0104   conversion_rate->Scale(1./nMB);
0105   file->ReOpen("UPDATE");
0106   conversion_rate->Write();
0107   return conversion_rate;
0108 }
0109
0110 void derivitvePlot(TH1F* finalrate){
0111   TH1F *dplot = new TH1F("derivative","",finalrate->GetNbinsX(),5,25);
0112   for (int i = 1; i < finalrate->GetNbinsX(); ++i)
0113   {
0114     double error;
0115     dplot->SetBinContent(i,finalrate->IntegralAndError(i,finalrate->GetNbinsX(),error));
0116     dplot->SetBinError(i,error);
0117   }
0118   dplot->Write();
0119 }
0120
0121 void conversionRate(){
0122   TFile *out_file = new TFile("effplots.root","UPDATE");
0123   TEfficiency* uni_rate=makepTRes(out_file);
0124   string treeExtension = ".root";
0125   if(!uni_rate){
0126     string treePath = "/sphenix/user/vassalli/gammasample/truthconversiononlineanalysis";
0127     unsigned int nFiles=200;
0128     TChain *ttree = handleFile(treePath,treeExtension,"cutTreeSignal",nFiles);
0129     TChain *observations = handleFile(treePath,treeExtension,"observTree",nFiles);
0130     cout<<"Got tree: "<<ttree->GetName()<<" and "<<ttree->GetEntries()<<" entries"<<endl;
0131     cout<<"Got tree: "<<observations->GetName()<<" and "<<observations->GetEntries()<<" entries"<<endl;
0132     uni_rate=makepTRes(out_file,ttree,observations);
0133   }
0134   out_file->ReOpen("READ");
0135   string treePath="/sphenix/user/vassalli/minBiasPythia/";
0136   string softname = "softana";
0137   string hardname = "hard4ana";
0138   auto conversion_rate=calculateRate(uni_rate,out_file,totalMB(treePath,treeExtension,softname,hardname));
0139   derivitvePlot(conversion_rate);
0140 }