h3l_pair_selector.C

Go to the documentation of this file.

 
#include "TFile.h"
#include "TH1F.h"
#include "TTreeReader.h"
#include "TTreeReaderArray.h"
#include "TTreeReaderValue.h"
 
#include <TSystem.h>
// #include "BmnHypNuclPair.h"
 
using namespace std;
 
Float_t invMassMin = 2.95;
Float_t invMassMax = 3.05;
 
// best by me for EXP
Float_t pathMin = 0.3;
Float_t pathMax = 35;
Float_t dca2Min = 5.7;
Float_t dca2Max = 14;
Float_t dca12Min = 0.0;
Float_t dca12Max = 5.0;
Float_t mom1Min = 4.5;
Float_t mom1Max = 15.5;
Float_t mom2Min = 0.3;
Float_t mom2Max = 1.1;
Float_t dca1Min = 0.00;
Float_t dca1Max = 1000.0;
Float_t dca0Min = 0.11;
Float_t dca0Max = 2.5;
Int_t fsd1Min = 0;
Int_t fsd1Max = 4;
Int_t fsd2Min = 0;
Int_t fsd2Max = 4;
Int_t gem1Min = 5;
Int_t gem1Max = 7;
Int_t gem2Min = 1;
Int_t gem2Max = 7;
 
// best by me for MC
// Float_t pathMin = 1.1;
// Float_t pathMax = 35;
// Float_t dca2Min = 3.6;
// Float_t dca2Max = 100;
// Float_t dca12Min = 0.0;
// Float_t dca12Max = 6.0;
// Float_t mom1Min = 4.5;
// Float_t mom1Max = 10.5;
// Float_t mom2Min = 0.0;
// Float_t mom2Max = 1.1;
// Float_t dca1Min = 0.00;
// Float_t dca1Max = 1000.0;
// Float_t dca0Min = 0.08;
// Float_t dca0Max = 0.65;
// Int_t fsd1Min =  0 ;
// Int_t fsd1Max =  4 ;
// Int_t fsd2Min =  0 ;
// Int_t fsd2Max =  4 ;
// Int_t gem1Min =  3 ;
// Int_t gem1Max =  7 ;
// Int_t gem2Min =  0 ;
// Int_t gem2Max =  7 ;
 
void h3l_pair_selector(TString inputFileName = "/home/merz/data/h3l_sim/h3lPairs_XeCsI_3.8AGeV_mb_4mln.root",
                       TString outputFileName = "/home/merz/data/h3l_sim/h3lSelect_XeCsI_3.8AGeV_mb_4mln.root")
{
 
    if (gSystem->AccessPathName(inputFileName))
        return;
 
    cout << "File: " << inputFileName << endl;
    TFile* inFile = TFile::Open(inputFileName);
 
    TTreeReader myReader("bmndata", inFile);
    TTreeReaderArray<BmnHypNuclPair> pairs(myReader, "pairs");
 
    TFile outFile(outputFileName, "recreate");
    vector<Float_t> Path;
    vector<Float_t> Dca0;
    vector<Float_t> Dca1;
    vector<Float_t> Dca2;
    vector<Float_t> Dca12;
    vector<Float_t> He3Mass;
    vector<Float_t> Dedx;
    vector<Float_t> InvMass;
    vector<Int_t> NHitsSil1;
    vector<Int_t> NHitsGem1;
    vector<Int_t> NHitsSil2;
    vector<Int_t> NHitsGem2;
    vector<Float_t> Px1;
    vector<Float_t> Py1;
    vector<Float_t> Pz1;
    vector<Float_t> Px2;
    vector<Float_t> Py2;
    vector<Float_t> Pz2;
    vector<Int_t> NTrInEvent;
 
    TTree t1("H3L", "H3L");
    t1.Branch("path", &Path);
    t1.Branch("dca0", &Dca0);
    t1.Branch("dca1", &Dca1);
    t1.Branch("dca2", &Dca2);
    t1.Branch("dca12", &Dca12);
    t1.Branch("He3Mass", &He3Mass);
    t1.Branch("dedx", &Dedx);
    t1.Branch("invMass", &InvMass);
    t1.Branch("nHitsSil1", &NHitsSil1);
    t1.Branch("nHitsGem1", &NHitsGem1);
    t1.Branch("nHitsSil2", &NHitsSil2);
    t1.Branch("nHitsGem2", &NHitsGem2);
    t1.Branch("Px1", &Px1);
    t1.Branch("Py1", &Py1);
    t1.Branch("Pz1", &Pz1);
    t1.Branch("Px2", &Px2);
    t1.Branch("Py2", &Py2);
    t1.Branch("Pz2", &Pz2);
    t1.Branch("nTrInEvent", &NTrInEvent);
 
    int evId = 0;
 
    while (myReader.Next()) {
        Path.clear();
        Dca0.clear();
        Dca1.clear();
        Dca2.clear();
        Dca12.clear();
        He3Mass.clear();
        Dedx.clear();
        InvMass.clear();
        NHitsSil1.clear();
        NHitsGem1.clear();
        NHitsSil2.clear();
        NHitsGem2.clear();
        Px1.clear();
        Py1.clear();
        Pz1.clear();
        Px2.clear();
        Py2.clear();
        Pz2.clear();
        NTrInEvent.clear();
 
        if (evId % 10000 == 0)
            printf("evId = %d\n", evId);
        evId++;
        for (BmnHypNuclPair pair : pairs) {
            Float_t mass = pair.GetInvMass();
            if (mass > invMassMax)
                continue;
            if (mass < invMassMin)
                continue;
            Float_t path = pair.GetPath();
            Float_t dedx = pair.GetDedx();
            Float_t mom1 = pair.GetMom1();
            Float_t dca0 = pair.GetDCA0();
            Float_t massTof = pair.GetM2();
            Float_t dca12 = pair.GetDCA12() / pair.GetdDCA12();
            Float_t dca1 = pair.GetDCA1() / pair.GetdDCA1();
            Float_t dca2 = pair.GetDCA2() / pair.GetdDCA2();
            Float_t mom2 = pair.GetMom2();
            Int_t nH1fsd = pair.GetNHitsPart1("SILICON");
            Int_t nH2fsd = pair.GetNHitsPart2("SILICON");
            Int_t nH1gem = pair.GetNHitsPart1("GEM");
            Int_t nH2gem = pair.GetNHitsPart2("GEM");
            Float_t px1 = pair.GetPx1();
            Float_t py1 = pair.GetPy1();
            Float_t pz1 = pair.GetPz1();
            Float_t px2 = pair.GetPx2();
            Float_t py2 = pair.GetPy2();
            Float_t pz2 = pair.GetPz2();
            Int_t nTr = pair.GetNTrInEvent();
 
            if (nH1fsd < fsd1Min || nH1fsd > fsd1Max)
                continue;
            if (nH2fsd < fsd2Min || nH2fsd > fsd2Max)
                continue;
            if (nH1gem < gem1Min || nH1gem > gem1Max)
                continue;
            if (nH2gem < gem2Min || nH2gem > gem2Max)
                continue;
            if (path < pathMin || path > pathMax)
                continue;
            // if (dca1 < dca1Min || dca1 > dca1Max) continue;
            if (dca2 < dca2Min || dca2 > dca2Max)
                continue;
            if (dca12 < dca12Min || dca12 > dca12Max)
                continue;
            if (mom1 < mom1Min || mom1 > mom1Max)
                continue;
            if (mom2 < mom2Min || mom2 > mom2Max)
                continue;
            if (dca0 < dca0Min || dca0 > dca0Max)
                continue;
 
            TLorentzVector lPos, lNeg;
            lPos.SetXYZM(px1, py1, pz1, 2.8094135);
            lNeg.SetXYZM(px2, py2, pz2, 0.1395704);
            TLorentzVector lHyp = lPos + lNeg;
            Float_t pt = lHyp.Pt();
            if (pt < 0.2)
                continue;
 
            Path.push_back(path);
            Dca0.push_back(dca0);
            Dca1.push_back(dca1);
            Dca2.push_back(dca2);
            Dca12.push_back(dca12);
            He3Mass.push_back(massTof);
            Dedx.push_back(dedx);
            InvMass.push_back(mass);
            NHitsSil1.push_back(nH1fsd);
            NHitsGem1.push_back(nH1gem);
            NHitsSil2.push_back(nH2fsd);
            NHitsGem2.push_back(nH2gem);
            Px1.push_back(px1);
            Py1.push_back(py1);
            Pz1.push_back(pz1);
            Px2.push_back(px2);
            Py2.push_back(py2);
            Pz2.push_back(pz2);
            NTrInEvent.push_back(nTr);
        }
        t1.Fill();
    }
 
    inFile->Close();
 
    t1.Write();
    outFile.Close();
}