h3lAna.C

Go to the documentation of this file.

 
#include "TFile.h"
#include "TH1F.h"
#include "TTreeReader.h"
#include "TTreeReaderValue.h"
 
#include <TLorentzVector.h>
#include <TSystem.h>
// #include "BmnHypNuclPair.h"
 
using namespace std;
using namespace TMath;
 
Float_t invMassMin = 2.95;
Float_t invMassMax = 3.05;
 
const Int_t nPath = 1;
const Int_t nDca2 = 1;
const Int_t nDca12 = 1;
const Int_t nDca0 = 1;
 
// EXP max by S
// vector<Float_t> pathMin = {1.0};
// vector<Float_t> pathMax = {35};
// vector<Float_t> dca2Min = {0.0};
// vector<Float_t> dca2Max = {100};
// vector<Float_t> dca12Min = {0.0};
// vector<Float_t> dca12Max = {9};
// vector<Float_t> dca0Min = {0.08};
// vector<Float_t> dca0Max = {0.65};
 
// about 58% of max S and 19% of max B (SIGN = 8.49, dSIGN = 1.37, S = 1262, dS = 207, B = 20854, dB = 144)
//  vector<Float_t> pathMin = {1.4};
//  vector<Float_t> pathMax = {35};
//  vector<Float_t> dca2Min = {3.6};
//  vector<Float_t> dca2Max = {100};
//  vector<Float_t> dca12Min = {0.0};
//  vector<Float_t> dca12Max = {6.0};
//  vector<Float_t> dca0Min = {0.08};
//  vector<Float_t> dca0Max = {0.65};
 
// MC max S
// vector<Float_t> pathMin = {1.2};
// vector<Float_t> pathMax = {35};
// vector<Float_t> dca2Min = {0.4};
// vector<Float_t> dca2Max = {100};
// vector<Float_t> dca12Min = {0.0};
// vector<Float_t> dca12Max = {18};
// vector<Float_t> dca0Min = {0.08};
// vector<Float_t> dca0Max = {1.05};
 
// MC by EXP (projection through param scan) (SIGN = 7.8, dSIGN = 1.1, S = 136, dS = 22, B = 164, dB = 13)
//  vector<Float_t> pathMin = {1.1};
//  vector<Float_t> pathMax = {35};
//  vector<Float_t> dca2Min = {3.6};
//  vector<Float_t> dca2Max = {100};
//  vector<Float_t> dca12Min = {0.0};
//  vector<Float_t> dca12Max = {6.0};
//  vector<Float_t> dca0Min = {0.08};
//  vector<Float_t> dca0Max = {0.65};
 
vector<Float_t> pathMin = {0.3};
vector<Float_t> pathMax = {35};
vector<Float_t> dca2Min = {5.7};
vector<Float_t> dca2Max = {14};
vector<Float_t> dca12Min = {0.0};
vector<Float_t> dca12Max = {5.0};
vector<Float_t> dca0Min = {0.11};
vector<Float_t> dca0Max = {2.5};
 
const Int_t nMom1 = 1;
vector<Float_t> mom1Min = {4.5};
vector<Float_t> mom1Max = {15.5};
const Int_t nMom2 = 1;
vector<Float_t> mom2Min = {0.3};
vector<Float_t> mom2Max = {1.1};
 
const Int_t nFsd1 = 1;
vector<Int_t> fsd1Min = {0};
vector<Int_t> fsd1Max = {4};
const Int_t nFsd2 = 1;
vector<Int_t> fsd2Min = {0};
vector<Int_t> fsd2Max = {4};
const Int_t nGem1 = 1;
vector<Int_t> gem1Min = {5};
vector<Int_t> gem1Max = {7};
const Int_t nGem2 = 1;
vector<Int_t> gem2Min = {1};
vector<Int_t> gem2Max = {7};
 
void h3lAna(TString inputFileName = "/home/merz/data/h3l_sim/h3lPairs_XeCsI_3.8AGeV_mb_4mln.root",
            TString outputFileName = "/home/merz/data/h3l_sim/h3lPlots_XeCsI_3.8AGeV_mb_4mln.root")
{
 
    if (gSystem->AccessPathName(inputFileName))
        return;
 
    TH1F* hInvMass[nPath][nDca2][nDca12][nDca0][nMom1][nMom2][nFsd1][nFsd2][nGem1][nGem2];
    for (Int_t iPath = 0; iPath < nPath; iPath++) {
        for (Int_t iDca2 = 0; iDca2 < nDca2; iDca2++) {
            for (Int_t iDca12 = 0; iDca12 < nDca12; iDca12++) {
                for (Int_t iDca0 = 0; iDca0 < nDca0; iDca0++) {
                    for (Int_t iMom1 = 0; iMom1 < nMom1; iMom1++) {
                        for (Int_t iMom2 = 0; iMom2 < nMom2; iMom2++) {
                            for (Int_t iFsd1 = 0; iFsd1 < nFsd1; iFsd1++) {
                                for (Int_t iFsd2 = 0; iFsd2 < nFsd2; iFsd2++) {
                                    for (Int_t iGem1 = 0; iGem1 < nGem1; iGem1++) {
                                        for (Int_t iGem2 = 0; iGem2 < nGem2; iGem2++) {
                                            hInvMass[iPath][iDca2][iDca12][iDca0][iMom1][iMom2][iFsd1][iFsd2][iGem1]
                                                    [iGem2] = new TH1F(Form("hInvMass_%d_%d_%d_%d_%d_%d_%d_%d_%d_%d",
                                                                            iPath, iDca2, iDca12, iDca0, iMom1, iMom2,
                                                                            iFsd1, iFsd2, iGem1, iGem2),
                                                                       "; M; Count", 100, invMassMin, invMassMax);
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
 
    cout << "File: " << inputFileName << endl;
    TFile* myFile = TFile::Open(inputFileName);
 
    TTreeReader myReader("bmndata", myFile);
    TTreeReaderArray<BmnHypNuclPair> pairs(myReader, "pairs");
 
    int evId = 0;
 
    while (myReader.Next()) {
        if (evId % 50000 == 0)
            printf("evId = %d\n", evId);
        evId++;
        for (BmnHypNuclPair pair : pairs) {
 
            // TLorentzVector lPos, lNeg;
            // Float_t px1 = pair.GetPx1();
            // Float_t py1 = pair.GetPy1();
            // Float_t pz1 = pair.GetPz1();
            // Float_t tx1 = px1 / pz1;
            // Float_t ty1 = py1 / pz1;
            // Float_t mom1_new = mom1;   // - 50.0 * corrFit->Eval(mom1);
            // //cout << mom1 << "\t" << corrFit->Eval(mom1) << "\t" << mom1_new << endl;
            // Float_t pz1_new = mom1_new / TMath::Sqrt(tx1 * tx1 + ty1 * ty1 + 1);
            // Float_t px1_new = pz1_new * tx1;
            // Float_t py1_new = pz1_new * ty1;
            // Float_t px2 = pair.GetPx2();
            // Float_t py2 = pair.GetPy2();
            // Float_t pz2 = pair.GetPz2();
 
            // // if ((px1 * px2) > 0 || (py1 * py2) > 0)
            // //     continue;
 
            // //cout << px1 << "\t" << py1 << "\t" << pz1 << endl;
            // //cout << px2 << "\t" << py2 << "\t" << pz2 << endl;
            // //cout << px1_new << "\t" << py1_new << "\t" << pz1_new << endl;
            // 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.6)
            //     continue;
 
            // Float_t mass_new = TLorentzVector((lPos + lNeg)).Mag();
            Float_t mass = pair.GetInvMass();
            // Float_t mass_new = mass;
            // cout << mass << "\t" << mass_new - mass << "\t" << mass_new << endl;
            if (mass > invMassMax)
                continue;
            if (mass < invMassMin)
                continue;
            Float_t mom1 = pair.GetMom1();
            Float_t mom2 = pair.GetMom2();
            // Float_t massTof = pair.GetM2_700_1();
            //  if (massTof < m2CutLow->Eval(mom1_new) || massTof > m2CutUp->Eval(mom1_new))
            //      continue;
            Float_t path = pair.GetPath();
            Float_t dca0 = pair.GetDCA0();
            Float_t dca12 = pair.GetDCA12() / pair.GetdDCA12();
            Float_t dca2 = pair.GetDCA2() / pair.GetdDCA2();
            Float_t dca1 = pair.GetDCA1() / pair.GetdDCA1();
 
            Int_t fsd1 = pair.GetNHitsPart1("SILICON");
            Int_t fsd2 = pair.GetNHitsPart2("SILICON");
            Int_t gem1 = pair.GetNHitsPart1("GEM");
            Int_t gem2 = pair.GetNHitsPart2("GEM");
 
            for (Int_t iPath = 0; iPath < nPath; iPath++) {
                if (path < pathMin[iPath] || path > pathMax[iPath])
                    continue;
                for (Int_t iDca2 = 0; iDca2 < nDca2; iDca2++) {
                    if (dca2 < dca2Min[iDca2] || dca2 > dca2Max[iDca2])
                        continue;
                    for (Int_t iDca12 = 0; iDca12 < nDca12; iDca12++) {
                        if (dca12 < dca12Min[iDca12] || dca12 > dca12Max[iDca12])
                            continue;
                        for (Int_t iDca0 = 0; iDca0 < nDca0; iDca0++) {
                            if (dca0 < dca0Min[iDca0] || dca0 > dca0Max[iDca0])
                                continue;
                            for (Int_t iMom1 = 0; iMom1 < nMom1; iMom1++) {
                                if (mom1 < mom1Min[iMom1] || mom1 > mom1Max[iMom1])
                                    continue;
                                for (Int_t iMom2 = 0; iMom2 < nMom2; iMom2++) {
                                    if (mom2 < mom2Min[iMom2] || mom2 > mom2Max[iMom2])
                                        continue;
                                    for (Int_t iFsd1 = 0; iFsd1 < nFsd1; iFsd1++) {
                                        if (fsd1 < fsd1Min[iFsd1] || fsd1 > fsd1Max[iFsd1])
                                            continue;
                                        for (Int_t iFsd2 = 0; iFsd2 < nFsd2; iFsd2++) {
                                            if (fsd2 < fsd2Min[iFsd2] || fsd2 > fsd2Max[iFsd2])
                                                continue;
                                            for (Int_t iGem1 = 0; iGem1 < nGem1; iGem1++) {
                                                if (gem1 < gem1Min[iGem1] || gem1 > gem1Max[iGem1])
                                                    continue;
                                                for (Int_t iGem2 = 0; iGem2 < nGem2; iGem2++) {
                                                    if (gem2 < gem2Min[iGem2] || gem2 > gem2Max[iGem2])
                                                        continue;
                                                    hInvMass[iPath][iDca2][iDca12][iDca0][iMom1][iMom2][iFsd1][iFsd2]
                                                            [iGem1][iGem2]
                                                                ->Fill(mass);
                                                }
                                            }
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
 
    myFile->Close();
 
    TFile outFile(outputFileName, "RECREATE");
    for (Int_t iPath = 0; iPath < nPath; iPath++) {
        for (Int_t iDca2 = 0; iDca2 < nDca2; iDca2++) {
            for (Int_t iDca12 = 0; iDca12 < nDca12; iDca12++) {
                for (Int_t iDca0 = 0; iDca0 < nDca0; iDca0++) {
                    for (Int_t iMom1 = 0; iMom1 < nMom1; iMom1++) {
                        for (Int_t iMom2 = 0; iMom2 < nMom2; iMom2++) {
                            for (Int_t iFsd1 = 0; iFsd1 < nFsd1; iFsd1++) {
                                for (Int_t iFsd2 = 0; iFsd2 < nFsd2; iFsd2++) {
                                    for (Int_t iGem1 = 0; iGem1 < nGem1; iGem1++) {
                                        for (Int_t iGem2 = 0; iGem2 < nGem2; iGem2++) {
                                            hInvMass[iPath][iDca2][iDca12][iDca0][iMom1][iMom2][iFsd1][iFsd2][iGem1]
                                                    [iGem2]
                                                        ->Write();
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    outFile.Close();
}