BmnRescale.cxx

Go to the documentation of this file.

#include "BmnRescale.h"
 
BmnRescale::BmnRescale(
        UInt_t periodId, UInt_t runId,
        Double_t lowThr,
        Double_t ClusterSizeThr,
        Int_t nBins) {
    fPeriodId = periodId;
    fRunId = runId;
    fSetup = (fRunId >= 2041 && fRunId <= 3588) ? kSRCSETUP : kBMNSETUP;
    fLowThr = lowThr;
    fClusterSizeThr = ClusterSizeThr;
    fNBins = nBins;
    CreateStripVectors();
}
 
BmnRescale::~BmnRescale() {
    FreeVector3D<BmnSigInfo>(fStripGemInfoMC);
    FreeVector3D<BmnSigInfo>(fStripSilInfoMC);
    FreeVector3D<BmnSigInfo>(fStripCSCInfoMC);
    FreeVector3D<BmnSigInfo>(fStripGemInfoEx);
    FreeVector3D<BmnSigInfo>(fStripSilInfoEx);
    FreeVector3D<BmnSigInfo>(fStripCSCInfoEx);
    FreeVector3D<TF1>(fStripGemRescale);
    FreeVector3D<TF1>(fStripSilRescale);
    FreeVector3D<TF1>(fStripCSCRescale);
    delete fGemStationSet;
    delete fSilStationSet;
    delete fCscStationSet;
}
 
void BmnRescale::CreateStripVectors() {
    TString xmlConfFileNameGEM;
    TString xmlConfFileNameSil;
    TString xmlConfFileNameCSC;
    // Select xml configs
    switch (fPeriodId) {
        case 7:
            if (fSetup == kBMNSETUP) {
                xmlConfFileNameGEM = "GemRunSpring2018.xml";
                xmlConfFileNameSil = "SiliconRunSpring2018.xml";
                xmlConfFileNameCSC = "CSCRunSpring2018.xml";
            } else {
                xmlConfFileNameGEM = "GemRunSRCSpring2018.xml";
                xmlConfFileNameSil = "SiliconRunSRCSpring2018.xml";
                xmlConfFileNameCSC = "CSCRunSRCSpring2018.xml";
            }
            break;
        case 6:
            xmlConfFileNameGEM = "GemRunSpring2017.xml";
            xmlConfFileNameSil = "SiliconRunSpring2017.xml";
            break;
        default:
            printf("Error! Unknown config!\n");
            return;
            break;
 
    }
    // create 3d vectors for functions for all stations/modules/planes
    if (xmlConfFileNameGEM.Length() > 0) {
        xmlConfFileNameGEM = TString(getenv("VMCWORKDIR")) + "/parameters/gem/XMLConfigs/" + xmlConfFileNameGEM;
        printf("xmlConfFileName %s\n", xmlConfFileNameGEM.Data());
        fGemStationSet = new BmnGemStripStationSet(xmlConfFileNameGEM);
        FillInfoVector<BmnGemStripStationSet, BmnGemStripStation, BmnGemStripModule>(
                fGemStationSet, fStripGemInfoMC);
        FillInfoVector<BmnGemStripStationSet, BmnGemStripStation, BmnGemStripModule>(
                fGemStationSet, fStripGemInfoEx);
    }
    if (xmlConfFileNameSil.Length() > 0) {
        xmlConfFileNameSil = TString(getenv("VMCWORKDIR")) + "/parameters/silicon/XMLConfigs/" + xmlConfFileNameSil;
        printf("xmlConfFileName %s\n", xmlConfFileNameSil.Data());
        fSilStationSet = new BmnSiliconStationSet(xmlConfFileNameSil);
        FillInfoVector<BmnSiliconStationSet, BmnSiliconStation, BmnSiliconModule>(
                fSilStationSet, fStripSilInfoMC);
        FillInfoVector<BmnSiliconStationSet, BmnSiliconStation, BmnSiliconModule>(
                fSilStationSet, fStripSilInfoEx);
 
    }
    if (xmlConfFileNameCSC.Length() > 0) {
        xmlConfFileNameCSC = TString(getenv("VMCWORKDIR")) + "/parameters/csc/XMLConfigs/" + xmlConfFileNameCSC;
        printf("xmlConfFileName %s\n", xmlConfFileNameCSC.Data());
        fCscStationSet = new BmnCSCStationSet(xmlConfFileNameCSC);
        FillInfoVector<BmnCSCStationSet, BmnCSCStation, BmnCSCModule>(
                fCscStationSet, fStripCSCInfoMC);
        FillInfoVector<BmnCSCStationSet, BmnCSCStation, BmnCSCModule>(
                fCscStationSet, fStripCSCInfoEx);
    }
}
 
TF1* BmnRescale::GetRescaleFunc(TString name, TF1 *mc, TF1 *ex) {
    if (!mc || !ex)
        return nullptr;
    TF1 *funcRescale = new TF1(name,
            [mc, ex](Double_t *x, Double_t * p) {
                Double_t r = mc->Eval(x[0]);
                Double_t f = ex->GetX(r); //, ex->GetXmin(), ex->GetXmax(), 1.E-15, 100, kTRUE);
                //                printf("x = %f  mc = %f  y = %f\n", x[0], r, f);
                return f;
            }, mc->GetXmin(), mc->GetXmax(), 0);
    return funcRescale;
}
 
BmnStatus BmnRescale::LoadDigitDistributionsFromFile(
        TString fileName,
        vector<TString> digiNames,
        vector<vector<vector<vector<BmnSigInfo* > > > > digiVecs) {
    vector<TClonesArray*> ars;
    TChain* chain = new TChain("bmndata");
    chain->Add(fileName.Data());
    Long64_t NEventsMC = chain->GetEntries();
    printf("%s recorded entries = %lld\n", fileName.Data(), NEventsMC);
    for (TString &digiName : digiNames) {
        TClonesArray* arDigi = nullptr;
        chain->SetBranchAddress(digiName.Data(), &arDigi);
        ars.push_back(arDigi);
    }
    for (size_t i = 0; i < digiVecs.size(); i++)
        CreateDigitDistributions(chain, ars[i], digiVecs[i]);
    delete chain;
    return kBMNSUCCESS;
}
 
BmnStatus BmnRescale::LoadClusterDistributionsFromFile(
        TString fileName,
        vector<TString> hitNames,
        vector<vector<vector<vector<BmnSigInfo* > > > > clusterVecs) {
    vector<TClonesArray*> hitArr;
    TChain* chain = new TChain("bmndata");
    chain->Add(fileName.Data());
    Long64_t NEvents = chain->GetEntries();
    printf("%s recorded entries = %lld\n", fileName.Data(), NEvents);
    for (TString &hitName : hitNames) {
        TClonesArray* arDigi = nullptr;
        chain->SetBranchAddress(hitName.Data(), &arDigi);
        hitArr.push_back(arDigi);
    }
    //    for (Int_t i = 0; i < fNArs; i++)
    //        CreateClusterDistributions(chain, ars[i], clusterVecs[i]);
    CreateClusterDistributions<
            BmnGemStripHit,
            BmnGemStripStationSet,
            BmnGemStripStation,
            BmnGemStripModule,
            BmnGemStripLayer>(chain, hitArr[0], fGemStationSet, clusterVecs[0]);
    CreateClusterDistributions<
            BmnSiliconHit,
            BmnSiliconStationSet,
            BmnSiliconStation,
            BmnSiliconModule,
            BmnSiliconLayer>(chain, hitArr[1], fSilStationSet, clusterVecs[1]);
    CreateClusterDistributions<
            BmnCSCHit,
            BmnCSCStationSet,
            BmnCSCStation,
            BmnCSCModule,
            BmnCSCLayer>(chain, hitArr[2], fCscStationSet, clusterVecs[2]);
    delete chain;
    return kBMNSUCCESS;
}
 
BmnStatus BmnRescale::LoadClustersOfTrackDistributionsFromFile(
        TString fileName,
        TString globalTrackName,
        vector<TString> subTrackNames,
        vector<TString> hitNames,
        vector<vector<vector<vector<BmnSigInfo* > > > > clusterVecs) {
    TClonesArray* globalTrackArr;
    vector<TClonesArray*> subTrackArr;
    vector<TClonesArray*> hitArr;
    TChain* chain = new TChain("bmndata");
    chain->Add(fileName.Data());
    Long64_t NEvents = chain->GetEntries();
    printf("%s recorded entries = %lld\n", fileName.Data(), NEvents);
 
    chain->SetBranchAddress(globalTrackName.Data(), &globalTrackArr);
 
    for (TString &subTrackName : subTrackNames) {
        TClonesArray* arDigi = nullptr;
        if (subTrackName.Length() > 0)
            chain->SetBranchAddress(subTrackName.Data(), &arDigi);
        subTrackArr.push_back(arDigi);
    }
    for (TString &hitName : hitNames) {
        TClonesArray* arDigi = nullptr;
        chain->SetBranchAddress(hitName.Data(), &arDigi);
        hitArr.push_back(arDigi);
    }
    //    for (Int_t i = 0; i < fNArs; i++)
    //        CreateClusterDistributions(chain, ars[i], clusterVecs[i]);
    CreateClusterDistributions<
            BmnGemStripHit,
            BmnGemStripStationSet,
            BmnGemStripStation,
            BmnGemStripModule,
            BmnGemStripLayer>(chain, globalTrackArr, subTrackArr[0], hitArr[0], fGemStationSet, clusterVecs[0]);
    CreateClusterDistributions<
            BmnSiliconHit,
            BmnSiliconStationSet,
            BmnSiliconStation,
            BmnSiliconModule,
            BmnSiliconLayer>(chain, globalTrackArr, subTrackArr[1], hitArr[1], fSilStationSet, clusterVecs[1]);
    CreateClusterDistributions<
            BmnCSCHit,
            BmnCSCStationSet,
            BmnCSCStation,
            BmnCSCModule,
            BmnCSCLayer>(chain, globalTrackArr, nullptr, hitArr[2], fCscStationSet, clusterVecs[2]);
    delete chain;
    return kBMNSUCCESS;
}
 
BmnStatus BmnRescale::CreateRescales(TString fileNameMC, TString fileNameEx) {
    if (fileNameMC == "" || fileNameEx == "") {
        cout << "Files need to be specified!" << endl;
        return kBMNERROR;
    }
    TString GlobalTrackName = "BmnGlobalTrack";
    vector<TString> mcDigiNames = {"BmnGemStripDigit", "BmnSiliconDigit", "BmnCSCDigit"};
    auto mcVecs = {fStripGemInfoMC, fStripSilInfoMC, fStripCSCInfoMC};
    LoadDigitDistributionsFromFile(fileNameMC, mcDigiNames, mcVecs);
 
    vector<TString> subTrackNames = {"BmnGemTrack", "BmnSiliconTrack", ""};
    vector<TString> hitNames = {"BmnGemStripHit", "BmnSiliconHit", "BmnCSCHit"};
    //    vector<TString> expDigiNames = {"SILICON", "GEM", "CSC"};
    //    LoadDigitDistributionsFromFile(fileNameEx, expDigiNames, expVecs);
    auto expVecs = {fStripGemInfoEx, fStripSilInfoEx, fStripCSCInfoEx};
    LoadClustersOfTrackDistributionsFromFile(fileNameEx, GlobalTrackName, subTrackNames, hitNames, expVecs);
 
    FillRescaleVector<BmnGemStripStationSet, BmnGemStripStation, BmnGemStripModule>(fGemStationSet, fStripGemInfoMC, fStripGemInfoEx, fStripGemRescale);
    FillRescaleVector<BmnSiliconStationSet, BmnSiliconStation, BmnSiliconModule>(fSilStationSet, fStripSilInfoMC, fStripSilInfoEx, fStripSilRescale);
    FillRescaleVector<BmnCSCStationSet, BmnCSCStation, BmnCSCModule>(fCscStationSet, fStripCSCInfoMC, fStripCSCInfoEx, fStripCSCRescale);
 
    return kBMNSUCCESS;
}
 
void BmnRescale::CreateDigitDistributions(
        TTree *treeDig, TClonesArray *digits,
        vector<vector<vector<BmnSigInfo* > > > &infoVec) {
    printf("CreateDigitDistributions started for %s\n", digits->GetClass()->GetName());
    // Find bounds
    for (Long64_t i = 0; i < treeDig->GetEntries(); i++) {
        treeDig->GetEntry(i);
        for (Int_t iHit = 0; iHit < digits->GetEntriesFast(); iHit++) {
            BmnStripDigit *dig = static_cast<BmnStripDigit*> (digits->UncheckedAt(iHit));
            Int_t s = dig->GetStation();
            Int_t m = dig->GetModule();
            Int_t l = dig->GetStripLayer();
            BmnSigInfo* info = infoVec[s][m][l];
            Double_t val = dig->GetStripSignal();
            if (val > fLowThr) {
                if (info->maxVal < val)
                    info->maxVal = val;
                if (info->minVal > val)
                    info->minVal = val;
            }
        }
    }
    // Create histograms
    for (size_t iStation = 0; iStation < infoVec.size(); iStation++) {
        printf("iStation %zu\n", iStation);
        for (size_t iModule = 0; iModule < infoVec[iStation].size(); iModule++) {
            printf("\tiModule %zu\n", iModule);
            for (size_t iLayer = 0; iLayer < infoVec[iStation][iModule].size(); iLayer++) {
                printf("\t\tiLayer %zu\n", iLayer);
                BmnSigInfo* info = infoVec[iStation][iModule][iLayer];
                Double_t minVal = info->minVal;
                Double_t maxVal = info->maxVal;
                if (minVal > maxVal)
                    continue;
                printf("\t\t\tMin  = %f max = %f\n", minVal, maxVal);
                TString name = Form("Sig_mc_%s_%zu_%zu_%zu",
                        digits->GetName(), iStation, iModule, iLayer);
                info->hSig = new TH1D(name, name, fNBins,
                        minVal - 0.5 * (maxVal - minVal) / (Double_t) fNBins,
                        maxVal + 0.5 * (maxVal - minVal) / (Double_t) fNBins);
                name = Form("IntSignal_%s_%zu_%zu_%zu",
                        digits->GetName(), iStation, iModule, iLayer);
                info->hIntSig = new TH1D(name, name, fNBins,
                        minVal - 0.5 * (maxVal - minVal) / (Double_t) fNBins,
                        maxVal + 0.5 * (maxVal - minVal) / (Double_t) fNBins);
            }
        }
    }
    // Fill signal hists
    for (Long64_t i = 0; i < treeDig->GetEntries(); i++) {
        treeDig->GetEntry(i);
        for (Int_t iHit = 0; iHit < digits->GetEntriesFast(); iHit++) {
            BmnStripDigit *dig = static_cast<BmnStripDigit*> (digits->UncheckedAt(iHit));
            Int_t s = dig->GetStation();
            Int_t m = dig->GetModule();
            Int_t l = dig->GetStripLayer();
            BmnSigInfo* info = infoVec[s][m][l];
            if (!(info->hSig))
                continue;
            Double_t val = dig->GetStripSignal();
            if (val > fLowThr) {
                info->hSig->Fill(val);
            }
        }
    }
    // Fill integral hists
    for (size_t iStation = 0; iStation < infoVec.size(); iStation++) {
        for (size_t iModule = 0; iModule < infoVec[iStation].size(); iModule++) {
            for (size_t iLayer = 0; iLayer < infoVec[iStation][iModule].size(); iLayer++) {
                BmnSigInfo* info = infoVec[iStation][iModule][iLayer];
                if (!(info->hSig))
                    continue;
                Double_t bc = 0;
                for (Int_t i = 0; i < fNBins; i++) {
                    bc += info->hSig->GetBinContent(i);
                    info->hIntSig->SetBinContent(i, bc);
                }
                if (bc == 0)
                    continue;
                info->hIntSig->Scale(1 / bc);
                info->hIntSig->SetLineColor(kRed);
                info->func = new TF1(
                        Form("%s_%zu_%zu_%zu", digits->GetName(), iStation, iModule, iLayer),
                        [info, bc](Double_t *x, Double_t * p) {
                            //            Double_t val = (bc - x[0])/bc;//log(x[0]);
                            if (x[0] >= info->maxVal) return 1.0;
                            Int_t iBin = info->hIntSig->GetXaxis()->FindBin(x[0]);
                            Int_t jBin = iBin + ((x[0] > info->hIntSig->GetBinCenter(iBin)) ? 1 : -1);
                            Double_t v = info->hIntSig->GetBinContent(iBin) +
                                    (info->hIntSig->GetBinContent(jBin) - info->hIntSig->GetBinContent(iBin))*
                                    2 * abs(x[0] - info->hIntSig->GetBinCenter(iBin)) /
                                    (info->hIntSig->GetBinWidth(iBin) + info->hIntSig->GetBinWidth(jBin));
                            return v;
                        }, info->minVal, info->maxVal, 0);
                info->func->SetNpx(fNBins);
                info->func->SetLineColor(kBlue);
            }
        }
    }
}
 
TCanvas * BmnRescale::GetCanvasGem() {
    return GetCanvas<
            BmnGemStripStationSet,
            BmnGemStripStation,
            BmnGemStripModule>(fGemStationSet);
}
 
TCanvas * BmnRescale::GetCanvasSil() {
    return GetCanvas<
            BmnSiliconStationSet,
            BmnSiliconStation,
            BmnSiliconModule>(fSilStationSet);
}
 
TCanvas * BmnRescale::GetCanvasCSC() {
    return GetCanvas<
            BmnCSCStationSet,
            BmnCSCStation,
            BmnCSCModule>(fCscStationSet);
}