LightFlavorRatios/swimming_correction/process_swimming.C

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