LightFlavorRatios/corrections/process_swimming.C

0001 void fill_histogram(TH1F* h, std::vector<float>* tps)
0002 {
0003   // even-index turning points are no->yes
0004   // odd-index turning points are yes->no
0005   for(int i=1; i<h->GetNbinsX(); i++)
0006   {
0007     float bincenter = h->GetBinCenter(i);
0008     int tp_index = -1;
0009     // find which turning point index applies to this bin
0010     for(int j=0; j<tps->size(); j++)
0011     {
0012       if(bincenter>tps->at(j))
0013       {
0014         tp_index = j;
0015       }
0016     }
0017     if(tp_index != -1 && (tp_index % 2) == 0)
0018     {
0019       h->Fill(bincenter);
0020     }
0021   }
0022 }
0023
0024 void process_swimming(std::string infile)
0025 {
0026   gStyle->SetOptStat(0);
0027
0028   TFile* f = TFile::Open(infile.c_str());
0029   TTree* t = (TTree*)f->Get("DecayTree");
0030
0031   std::vector<float>* tp_lifetime = nullptr;
0032   std::vector<float>* tp_decaylength = nullptr;
0033
0034   t->SetBranchAddress("turningPoints_TAU",&tp_lifetime);
0035   t->SetBranchAddress("turningPoints_FD",&tp_decaylength);
0036
0037   TH1F* decaylength_plot = new TH1F("h_decaylength","Normalized sum of decay-length step functions; 2D decay length [cm]",1000,0.,5.5);
0038   TH1F* lifetime_plot = new TH1F("h_lifetime","Normalized sum of lifetime step functions; decay time [ps]",1000,0.,500.);
0039
0040   for(int i=0; i<t->GetEntries(); i++)
0041   {
0042     if(i % 1000 == 0) std::cout << i << std::endl;
0043     t->GetEntry(i);
0044
0045     std::sort(tp_decaylength->begin(),tp_decaylength->end());
0046     std::sort(tp_lifetime->begin(),tp_lifetime->end());
0047 /*
0048     std::cout << "FD size " << tp_decaylength->size() << std::endl;
0049     for(float fd : *tp_decaylength) std::cout << fd << " ";
0050     std::cout << std::endl;
0051     std::cout << "TAU size " << tp_lifetime->size() << std::endl;
0052     for(float lt : *tp_lifetime) std::cout << lt << " ";
0053     std::cout << std::endl;
0054 */
0055     fill_histogram(lifetime_plot,tp_lifetime);
0056     fill_histogram(decaylength_plot,tp_decaylength);
0057   }
0058
0059   decaylength_plot->Scale(1./t->GetEntries());
0060   lifetime_plot->Scale(1./t->GetEntries());
0061   decaylength_plot->GetYaxis()->SetRangeUser(0.,1.);
0062   lifetime_plot->GetYaxis()->SetRangeUser(0.,1.);
0063
0064   TCanvas* cd = new TCanvas("fd","fd",600,600);
0065   decaylength_plot->Draw("L");
0066   TCanvas* ct = new TCanvas("ft","ft",600,600);
0067   lifetime_plot->Draw("L");
0068 }