LightFlavorRatios/swimming_correction/pythia_lifetime.C

0001 #include <phool/PHRandomSeed.h>
0002 #include <gsl/gsl_rng.h>
0003 #include <gsl/gsl_randist.h>
0004
0005 #include <iostream>
0006
0007 #include "TFile.h"
0008 #include "TTree.h"
0009 #include "TChain.h"
0010 #include "TH1F.h"
0011
0012 const float mass = 1.115;
0013 const float mass_resolution = 0.01; // absolute uncertainty [GeV]
0014 const float p_resolution_coef = 0.002; // sigma_p/p = coef*p
0015 const float vertex_resolution = 0.01; // absolute uncertainty [cm]
0016
0017 void pythia_lifetime(int process=1)
0018 {
0019   const uint seed = PHRandomSeed();
0020   std::unique_ptr<gsl_rng> m_rng;
0021   m_rng.reset(gsl_rng_alloc(gsl_rng_mt19937));
0022   gsl_rng_set(m_rng.get(), seed);
0023
0024   TChain* t = new TChain("ntp_gtrack");
0025   if(process==0)
0026   {
0027     t->Add("/sphenix/tg/tg01/hf/frawley/sims/pythia_fromhits_detroit/pythia_reco_3122_[0-9]_g4svtx_eval.root");
0028     t->Add("/sphenix/tg/tg01/hf/frawley/sims/pythia_fromhits_detroit/pythia_reco_3122_[0-9][0-9]_g4svtx_eval.root");
0029   }
0030   else
0031   {
0032     t->Add(("/sphenix/tg/tg01/hf/frawley/sims/pythia_fromhits_detroit/pythia_reco_3122_"+std::to_string(process)+"[0-9][0-9]_g4svtx_eval.root").c_str());
0033   }
0034
0035   TFile* fout = new TFile(("/sphenix/user/mjpeters/analysis/LightFlavorRatios/swimming_correction/output/fout_"+std::to_string(process)+".root").c_str(),"RECREATE");
0036
0037   int segment;
0038   int event_out;
0039
0040   float lifetime_out;
0041   float lifetime_labframe_out;
0042
0043   float mother_vx;
0044   float mother_vy;
0045   float mother_vz;
0046   float mother_vt;
0047
0048   float track1_vx;
0049   float track1_vy;
0050   float track1_vz;
0051   float track1_vt;
0052
0053   float track2_vx;
0054   float track2_vy;
0055   float track2_vz;
0056   float track2_vt;
0057
0058   float mother_pt;
0059   float mother_pz;
0060   float track1_pt;
0061   float track1_pz;
0062   float track2_pt;
0063   float track2_pz;
0064
0065   float decaylength_out;
0066   float decaylength_err_out;
0067   float p_err_out;
0068   float decaytime_err_out;
0069   float smear_amount_out;
0070
0071   bool mother_isReco;
0072   bool track1_isReco;
0073   bool track2_isReco;
0074   bool mother_isPrimary;
0075
0076   TTree* t_out = new TTree("lifetime_study","truth candidate lifetime");
0077   t_out->Branch("segment",&segment_out);
0078   t_out->Branch("event",&event_out);
0079   t_out->Branch("lifetime",&lifetime_out);
0080   t_out->Branch("lifetime_labframe",&lifetime_labframe_out);
0081
0082   t_out->Branch("mother_vx",&mother_vx);
0083   t_out->Branch("mother_vy",&mother_vy);
0084   t_out->Branch("mother_vz",&mother_vz);
0085   t_out->Branch("mother_vt",&mother_vt);
0086
0087   t_out->Branch("track1_vx",&track1_vx);
0088   t_out->Branch("track1_vy",&track1_vy);
0089   t_out->Branch("track1_vz",&track1_vz);
0090   t_out->Branch("track1_vt",&track1_vt);
0091
0092   t_out->Branch("track2_vx",&track2_vx);
0093   t_out->Branch("track2_vy",&track2_vy);
0094   t_out->Branch("track2_vz",&track2_vz);
0095   t_out->Branch("track2_vt",&track2_vt);
0096
0097   t_out->Branch("mother_pt",&mother_pt);
0098   t_out->Branch("mother_pz",&mother_pz);
0099   t_out->Branch("track1_pt",&track1_pt);
0100   t_out->Branch("track1_pz",&track1_pz);
0101   t_out->Branch("track2_pt",&track2_pt);
0102   t_out->Branch("track2_pz",&track2_pz);
0103
0104   t_out->Branch("decaylength",&decaylength_out);
0105   t_out->Branch("decaylength_err",&decaylength_err_out);
0106   t_out->Branch("p_err",&p_err_out);
0107   t_out->Branch("lifetime_err",&decaytime_err_out);
0108   t_out->Branch("smear_amount",&smear_amount_out);
0109
0110   t_out->Branch("mother_isReco",&mother_isReco);
0111   t_out->Branch("track1_isReco",&track1_isReco);
0112   t_out->Branch("track2_isReco",&track2_isReco);
0113   t_out->Branch("mother_isPrimary",&mother_isPrimary);
0114
0115   //std::string indir = "/sphenix/tg/tg01/hf/frawley/sims/pythia_fromhits_detroit/";
0116   TH1F* lifetime = new TH1F("lifetime","Truth lifetime;#tau [ns]",100,0.,0.5);
0117   TH1F* lifetime_ptcut = new TH1F("lifetime_ptcut","Truth lifetime after 160MeV daughter pT cut;#tau [ns]",100,0.,0.5);
0118   TH1F* lifetime_motherptcut = new TH1F("lifetime_motherptcut","Truth lifetime after 800MeV mother pT cut;#tau [ns]",100,0.,0.5);
0119   TH1F* lifetime_bothptcut = new TH1F("lifetime_bothptcut","Truth lifetime after both 160MeV daughter pT cut and 800MeV mother pT cut;#tau [ns]",100,0.,0.5);
0120   TH1F* lifetime_reco = new TH1F("lifetime_reco","Truth lifetime for reconstructible mothers;#tau [ns]",100,0.,0.5);
0121   TH1F* lifetime_fromlength = new TH1F("lifetime_fromlength","Truth lifetime from decay length;#tau [ns]",100,0.,0.5);
0122   TH1F* lifetime_diffmethods = new TH1F("lifetime_diffmethods","Truth lifetime difference from decay time - decay length;#Delta#tau [ns]",100,-0.02,0.02);
0123   TH1F* lifetime_smeared = new TH1F("lifetime_smeared","Truth lifetime after resolution smearing;#tau [ns]",100,0.,0.5);
0124   TH1F* lifetime_smear_amount = new TH1F("lifetime_smear_amount","Amount of smearing applied;#Delta#tau [ns]",100,-0.05,0.05);
0125   TH1F* lifetime_resolution = new TH1F("lifetime_resolution","Lifetime resolution;#sigma_{#tau} [ns]",100,0.,0.1);
0126
0127   TH1F* lifetime_labframe = new TH1F("lifetime_labframe","Lab-frame truth lifetime;t [ns]",100,0.,0.5);
0128   TH1F* lifetime_labframe_ptcut = new TH1F("lifetime_labframe_ptcut","Lab-frame truth lifetime after 160MeV daughter pT cut",100,0.,0.5);
0129   TH1F* lifetime_labframe_motherptcut = new TH1F("lifetime_labframe_motherptcut","Lab-frame truth lifetime after 800MeV mother pT cut",100,0.,0.5);
0130   TH1F* lifetime_labframe_bothptcut = new TH1F("lifetime_labframe_bothptcut","Lab-frame truth lifetime after both 160MeV daughter pT cut and 800MeV mother pT cut",100,0.,0.5);
0131
0132   TH2F* lifetime_vsmotherpt = new TH2F("lifetime_vsmotherpt","Truth lifetime vs. mother pT;t [ns];mother pT [GeV]",100,0.,0.5,100,0.,3.);
0133   TH2F* lifetime_vstrack1pt = new TH2F("lifetime_vstrack1pt","Truth lifetime vs. track 1 pT;t [ns];track 1 pT [GeV]",100,0.,0.5,100,0.,3.);
0134   TH2F* lifetime_vstrack2pt = new TH2F("lifetime_vstrack2pt","Truth lifetime vs. track 2 pT;t [ns];track 2 pT [GeV]",100,0.,0.5,100,0.,3.);
0135
0136   TH2F* lifetime_labframe_vsmotherpt = new TH2F("lifetime_labframe_vsmotherpt","Lab-frame truth lifetime vs. mother pT;t [ns];mother pT [GeV]",100,0.,0.5,100,0.,3.);
0137   TH2F* lifetime_labframe_vstrack1pt = new TH2F("lifetime_labframe_vstrack1pt","Lab-frame truth lifetime vs. track 1 pT;t [ns];track 1 pT [GeV]",100,0.,0.5,100,0.,3.);
0138   TH2F* lifetime_labframe_vstrack2pt = new TH2F("lifetime_labframe_vstrack2pt","Lab-frame truth lifetime vs. track 2 pT;t [ns];track 2 pT [Gev]",100,0.,0.5,100,0.,3.);
0139
0140   //for(int i=0; i<100; i++)
0141   //{
0142   //  std::cout << "file " << i << std::endl;
0143   //  std::string filename = indir+"pythia_reco_3122_"+std::to_string(i)+"_g4svtx_eval.root";
0144   //  TFile* f = TFile::Open(filename.c_str());
0145   //  TTree* t = (TTree*)f->Get("ntp_gtrack");
0146     float event;
0147     float gtrackid;
0148     float gflavor;
0149     float gparentflavor;
0150     float gparentid;
0151     float gvx;
0152     float gvy;
0153     float gvz;
0154     float gvt;
0155     float trackid;
0156     float gpt;
0157     float gpz;
0158     t->SetBranchAddress("event",&event);
0159     t->SetBranchAddress("gtrackID",&gtrackid);
0160     t->SetBranchAddress("gflavor",&gflavor);
0161     t->SetBranchAddress("gparentflavor",&gparentflavor);
0162     t->SetBranchAddress("gparentid",&gparentid);
0163     t->SetBranchAddress("gvx",&gvx);
0164     t->SetBranchAddress("gvy",&gvy);
0165     t->SetBranchAddress("gvz",&gvz);
0166     t->SetBranchAddress("gvt",&gvt);
0167     t->SetBranchAddress("trackID",&trackid);
0168     t->SetBranchAddress("gpt",&gpt);
0169     t->SetBranchAddress("gpz",&gpz);
0170
0171     std::vector<int> parent_ids;
0172     std::vector<int> parent_flavors;
0173     std::vector<float> primary_gvt;
0174     std::vector<float> parent_pt;
0175     std::vector<float> parent_pz;
0176     std::vector<float> parent_p;
0177     std::vector<std::array<float,3>> parent_vxyz;
0178
0179     std::vector<std::vector<int>> daughter_gtrackids;
0180     std::vector<std::vector<float>> daughter_trackids;
0181     std::vector<std::vector<int>> daughter_flavors;
0182     std::vector<std::vector<float>> daughter_pt;
0183     std::vector<std::vector<float>> daughter_pz;
0184     std::vector<std::vector<std::array<float,3>>> secondary_vxyz;
0185     std::vector<std::vector<float>> secondary_vts;
0186
0187     int current_file = 0;
0188     int total_nMothers = 0;
0189     int total_nDaughters = 0;
0190     int total_nValid = 0;
0191     int total_aboveptcut = 0;
0192     int total_abovemotherptcut = 0;
0193     int total_abovebothptcut = 0;
0194     int total_nReconstructible = 0;
0195
0196     int current_event = 0;
0197     int current_event_start = 0;
0198     int nMothers = 0;
0199     int nDaughters = 0;
0200
0201     // main event loop
0202     for(int j=0;j<t->GetEntries();j++)
0203     {
0204       t->GetEntry(j);
0205
0206       // check if event changed
0207       if(event!=current_event)
0208       {
0209
0210         if(event==0)
0211         {
0212           std::cout << "file " << current_file << std::endl;
0213           std::cout << "total mothers " << total_nMothers << std::endl;
0214           std::cout << "total daughters " << total_nDaughters << std::endl;
0215           std::cout << "total valid " << total_nValid << " above pt cut " << total_aboveptcut << " above mother pt cut " << total_abovemotherptcut << " above both pt cuts " << total_abovebothptcut << " reconstructible " << total_nReconstructible << std::endl;
0216           current_file++;
0217         }
0218         //std::cout << "event " << event << std::endl;
0219         //std::cout << "found " << nMothers << " mothers" << std::endl;
0220
0221         // process full event
0222         current_event = event;
0223         int current_event_end = j-1;
0224         // find daughters
0225         for(int k=current_event_start; k<=current_event_end; k++)
0226         {
0227           t->GetEntry(k);
0228           if((int)fabs(round(gparentflavor))==3122)
0229           {
0230             auto mother_it = std::find(parent_ids.begin(),parent_ids.end(),gparentid);
0231             if(mother_it != parent_ids.end())
0232             {
0233               const int mother_ind = std::distance(parent_ids.begin(),mother_it);
0234               daughter_gtrackids[mother_ind].push_back(round(gtrackid));
0235               daughter_trackids[mother_ind].push_back(trackid);
0236               daughter_flavors[mother_ind].push_back(round(gflavor));
0237               daughter_pt[mother_ind].push_back(gpt);
0238               daughter_pz[mother_ind].push_back(gpz);
0239               secondary_vxyz[mother_ind].push_back({gvx,gvy,gvz});
0240               secondary_vts[mother_ind].push_back(gvt);
0241               nDaughters++;
0242               total_nDaughters++;
0243             }
0244           }
0245         }
0246
0247         //std::cout << "found " << nDaughters << " daughters" << std::endl;
0248
0249         int nValid = 0;
0250         int nAbovePtCut = 0;
0251         int nAboveMotherPtCut = 0;
0252         int nAboveBothPtCut = 0;
0253         int nReconstructible = 0;
0254
0255         // filter for valid (has two daughters of correct flavor) and reconstructed (daughters both have non-NaN trackIDs) candidates
0256         // decay time = secondary vt - primary vt
0257         for(int lambda_ind = 0; lambda_ind<parent_ids.size(); lambda_ind++)
0258         {
0259           const int lambda_parity = (parent_flavors[lambda_ind]>0)? 1:-1;
0260
0261           // find daughter proton, pi, if they exist
0262           auto daughter_p_it = std::find(daughter_flavors[lambda_ind].begin(),daughter_flavors[lambda_ind].end(),2212*lambda_parity);
0263           auto daughter_pi_it = std::find(daughter_flavors[lambda_ind].begin(),daughter_flavors[lambda_ind].end(),-211*lambda_parity);
0264
0265           if(daughter_p_it != daughter_flavors[lambda_ind].end() && daughter_pi_it != daughter_flavors[lambda_ind].end())
0266           {
0267             nValid++;
0268             total_nValid++;
0269
0270             segment_out = process*100+current_file;
0271             event_out = current_event;
0272
0273             const int daughter_p_ind = std::distance(daughter_flavors[lambda_ind].begin(),daughter_p_it);
0274             const int daughter_pi_ind = std::distance(daughter_flavors[lambda_ind].begin(),daughter_pi_it);
0275
0276             mother_vx = parent_vxyz[lambda_ind][0];
0277             mother_vy = parent_vxyz[lambda_ind][1];
0278             mother_vz = parent_vxyz[lambda_ind][2];
0279             mother_vt = primary_gvt[lambda_ind];
0280
0281             track1_vx = secondary_vxyz[lambda_ind][daughter_p_ind][0];
0282             track1_vy = secondary_vxyz[lambda_ind][daughter_p_ind][1];
0283             track1_vz = secondary_vxyz[lambda_ind][daughter_p_ind][2];
0284             track1_vt = secondary_vts[lambda_ind][daughter_p_ind];
0285
0286             track2_vx = secondary_vxyz[lambda_ind][daughter_pi_ind][0];
0287             track2_vy = secondary_vxyz[lambda_ind][daughter_pi_ind][1];
0288             track2_vz = secondary_vxyz[lambda_ind][daughter_pi_ind][2];
0289             track2_vt = secondary_vts[lambda_ind][daughter_pi_ind];
0290
0291             const float lab_decay_time = secondary_vts[lambda_ind][daughter_p_ind] - primary_gvt[lambda_ind];
0292             lifetime_labframe->Fill(lab_decay_time);
0293             lifetime_labframe_vsmotherpt->Fill(lab_decay_time,parent_pt[lambda_ind]);
0294             lifetime_labframe_vstrack1pt->Fill(lab_decay_time,daughter_pt[lambda_ind][daughter_p_ind]);
0295             lifetime_labframe_vstrack2pt->Fill(lab_decay_time,daughter_pt[lambda_ind][daughter_pi_ind]);
0296
0297             const float betagamma = parent_p[lambda_ind]/1.115; // = sqrt(gamma^2 - 1)
0298             const float gamma = sqrt(pow(betagamma,2)+1.);
0299             const float proper_decay_time = lab_decay_time/gamma;
0300             //std::cout << "proper time: " << proper_decay_time << std::endl;
0301
0302             //std::cout << "PV: (" << parent_vxyz[lambda_ind][0] << ", " << parent_vxyz[lambda_ind][1] << ", " << parent_vxyz[lambda_ind][2] << "), SV (" << gvx << ", " << gvy << ", " << gvz << ")" << std::endl;
0303
0304             const float decay_length = sqrt(pow(secondary_vxyz[lambda_ind][daughter_p_ind][0] - parent_vxyz[lambda_ind][0],2) +
0305                                             pow(secondary_vxyz[lambda_ind][daughter_p_ind][1] - parent_vxyz[lambda_ind][1],2) +
0306                                             pow(secondary_vxyz[lambda_ind][daughter_p_ind][2] - parent_vxyz[lambda_ind][2],2));
0307             const float decay_length_err = 1./decay_length*sqrt(6.*pow(vertex_resolution,2));
0308             //std::cout << "decay length: " << decay_length << " +- " << decay_length_err << std::endl;
0309
0310             const float p_err = p_resolution_coef*pow(parent_p[lambda_ind],2);
0311
0312             const float speed = 2.99792458e1; // ctau [cm] -> ctau [ns]
0313
0314             const float proper_decay_time_fromdecaylength = 1./speed * mass*decay_length/parent_p[lambda_ind];
0315             const float proper_decay_time_fromdecaylength_err = 1./speed * sqrt(pow(decay_length/parent_p[lambda_ind],2)*pow(mass_resolution,2) +
0316                                                                                 pow(mass/parent_p[lambda_ind],2)*pow(decay_length_err,2) +
0317                                                                                 pow(mass*decay_length/pow(parent_p[lambda_ind],2),2)*pow(p_err,2));
0318
0319             lifetime_resolution->Fill(proper_decay_time_fromdecaylength_err);
0320
0321             //std::cout << "proper time from decay length: " << proper_decay_time_fromdecaylength << " +- " << proper_decay_time_fromdecaylength_err << std::endl;
0322
0323             const float proper_decay_time_resolution_smear = gsl_ran_gaussian(m_rng.get(),proper_decay_time_fromdecaylength_err);
0324             lifetime_smear_amount->Fill(proper_decay_time_resolution_smear);
0325             const float proper_decay_time_smeared = fabs(proper_decay_time_fromdecaylength + proper_decay_time_resolution_smear);
0326             //std::cout << "proper time, smeared: " << proper_decay_time_smeared << std::endl;
0327
0328             lifetime_out = proper_decay_time;
0329             lifetime_labframe_out = lab_decay_time;
0330             decaylength_out = decay_length;
0331             decaylength_err_out = decay_length_err;
0332             p_err_out = p_err;
0333             lifetime_err_out = proper_decay_time_fromdecaylength_err;
0334             smear_amount_out = proper_decay_time_resolution_smear;
0335
0336             track1_isReco = !std::isnan(daughter_trackids[lambda_ind][daughter_p_ind]);
0337             track2_isReco = !std::isnan(daughter_trackids[lambda_ind][daughter_pi_ind]);
0338             mother_isReco = track1_isReco && track2_isReco;
0339
0340             lifetime->Fill(proper_decay_time);
0341             lifetime_vsmotherpt->Fill(proper_decay_time,parent_pt[lambda_ind]);
0342             lifetime_vstrack1pt->Fill(proper_decay_time,daughter_pt[lambda_ind][daughter_p_ind]);
0343             lifetime_vstrack2pt->Fill(proper_decay_time,daughter_pt[lambda_ind][daughter_pi_ind]);
0344             lifetime_fromlength->Fill(proper_decay_time_fromdecaylength);
0345             lifetime_diffmethods->Fill(proper_decay_time - proper_decay_time_fromdecaylength);
0346             lifetime_smeared->Fill(proper_decay_time_smeared);
0347             if(daughter_pt[lambda_ind][daughter_p_ind]>0.16 && daughter_pt[lambda_ind][daughter_pi_ind]>0.16)
0348             {
0349               nAbovePtCut++;
0350               total_aboveptcut++;
0351               lifetime_ptcut->Fill(proper_decay_time);
0352               lifetime_labframe_ptcut->Fill(lab_decay_time);
0353               if(parent_pt[lambda_ind]>0.8)
0354               {
0355                 nAboveBothPtCut++;
0356                 total_abovebothptcut++;
0357                 lifetime_bothptcut->Fill(proper_decay_time);
0358                 lifetime_labframe_bothptcut->Fill(lab_decay_time);
0359               }
0360             }
0361             if(parent_pt[lambda_ind]>0.8)
0362             {
0363               nAboveMotherPtCut++;
0364               total_abovemotherptcut++;
0365               lifetime_motherptcut->Fill(proper_decay_time);
0366               lifetime_labframe_motherptcut->Fill(lab_decay_time);
0367             }
0368             if(!std::isnan(daughter_trackids[lambda_ind][daughter_p_ind]) && !std::isnan(daughter_trackids[lambda_ind][daughter_pi_ind]))
0369             {
0370               nReconstructible++;
0371               total_nReconstructible++;
0372               lifetime_reco->Fill(proper_decay_time);
0373             }
0374           }
0375         }
0376
0377         //std::cout << "of these, " << nValid << " were valid, " << nAbovePtCut << " were above pt cut, and " << nReconstructible << " were reconstructible" << std::endl;
0378
0379         // clean up for next event
0380
0381         for(int i=0; i<parent_ids.size(); i++)
0382         {
0383           daughter_gtrackids[i].clear();
0384           daughter_trackids[i].clear();
0385           daughter_flavors[i].clear();
0386           daughter_pt[i].clear();
0387           daughter_pz[i].clear();
0388           secondary_vxyz[i].clear();
0389           secondary_vts[i].clear();
0390         }
0391         daughter_gtrackids.clear();
0392         daughter_trackids.clear();
0393         daughter_flavors.clear();
0394         daughter_pt.clear();
0395         daughter_pz.clear();
0396         secondary_vxyz.clear();
0397         secondary_vts.clear();
0398         parent_ids.clear();
0399         parent_flavors.clear();
0400         primary_gvt.clear();
0401         parent_p.clear();
0402         parent_pt.clear();
0403         parent_pz.clear();
0404         parent_vxyz.clear();
0405
0406         nMothers = 0;
0407         nDaughters = 0;
0408         nValid = 0;
0409         nAbovePtCut = 0;
0410         nAboveMotherPtCut = 0;
0411         nAboveBothPtCut = 0;
0412         nReconstructible = 0;
0413         current_event_start = j;
0414         t->GetEntry(j);
0415       }
0416
0417       // find mothers
0418       if(gflavor==3122)
0419       {
0420         parent_ids.push_back(round(gtrackid));
0421         parent_flavors.push_back(round(gflavor));
0422         primary_gvt.push_back(gvt);
0423         parent_p.push_back(sqrt(pow(gpt,2)+pow(gpz,2)));
0424         parent_pt.push_back(gpt);
0425         parent_pz.push_back(gpz);
0426         parent_vxyz.push_back({gvx,gvy,gvz});
0427         daughter_gtrackids.emplace_back();
0428         daughter_trackids.emplace_back();
0429         daughter_flavors.emplace_back();
0430         daughter_pt.emplace_back();
0431         secondary_vxyz.emplace_back();
0432         secondary_vts.emplace_back();
0433         nMothers++;
0434         total_nMothers++;
0435       }
0436     }
0437   TCanvas* ctruth = new TCanvas("ctruth","ctruth",600,600);
0438   lifetime->Scale(1./lifetime->GetEntries());
0439   lifetime_reco->Scale(1./lifetime_reco->GetEntries());
0440   lifetime_ptcut->Scale(1./lifetime_ptcut->GetEntries());
0441   lifetime->Draw();
0442   lifetime_reco->SetLineColor(kRed);
0443   lifetime_reco->Draw("SAME");
0444   lifetime_ptcut->SetLineColor(kBlack);
0445   lifetime_ptcut->Draw("SAME");
0446   ctruth->SetLogy();
0447
0448   TCanvas* cdiff = new TCanvas("cdiff","cdiff",600,600);
0449   lifetime_diffmethods->Draw();
0450
0451   TCanvas* csmear = new TCanvas("csmear","csmear",600,600);
0452   lifetime_fromlength->Scale(1./lifetime_fromlength->GetEntries());
0453   lifetime_smeared->Scale(1./lifetime_fromlength->GetEntries());
0454   lifetime_fromlength->SetMaximum(1.1*std::max(lifetime_fromlength->GetMaximum(),lifetime_smeared->GetMaximum()));
0455   lifetime_fromlength->SetLineColor(kBlack);
0456   lifetime_fromlength->Draw();
0457   lifetime_smeared->SetLineColor(kRed);
0458   lifetime_smeared->Draw("SAME");
0459   csmear->SetLogy();
0460
0461   lifetime->Write();
0462   lifetime_reco->Write();
0463   lifetime_ptcut->Write();
0464   lifetime_motherptcut->Write();
0465   lifetime_bothptcut->Write();
0466   lifetime_fromlength->Write();
0467   lifetime_diffmethods->Write();
0468   lifetime_smeared->Write();
0469   lifetime_smear_amount->Write();
0470   lifetime_resolution->Write();
0471
0472   lifetime_labframe->Write();
0473   lifetime_labframe_ptcut->Write();
0474   lifetime_labframe_motherptcut->Write();
0475   lifetime_labframe_bothptcut->Write();
0476
0477   lifetime_vsmotherpt->Write();
0478   lifetime_vstrack1pt->Write();
0479   lifetime_vstrack2pt->Write();
0480   lifetime_labframe_vsmotherpt->Write();
0481   lifetime_labframe_vstrack1pt->Write();
0482   lifetime_labframe_vstrack2pt->Write();
0483
0484 }