fit_thresh_corrs.C

Go to the documentation of this file.

#include "TCanvas.h"
#include "TDecompChol.h"
#include "TDecompLU.h"
#include "TDecompSVD.h"
#include "TF1.h"
#include "TFile.h"
#include "TGraph.h"
#include "TH2.h"
#include "TMatrixD.h"
#include "TVectorD.h"
#include "TProfile.h"
#include "TTree.h"

// Order of fit in each dimension
const int Na = 12; // W
const int Nb = 6; // Cells
const int Nc = 4; // Position in module
const int Nflag = 3;

// NearDet settings
//const int cellMax = 32*3;
//const double cutW = 190;
//const double cutWMuonX = 90;
//const double cellCorrMin = 0.5;
//const double cellCorrMax = 2.0;
//const double wCorrMin = 0.5;
//const double wCorrMax = 2.0;

// TestBeam settings
const int cellMax = 32*2;
const double cutW = 114;
const double cutWMuonX = 90;//Needs to be set but is unused. Value is irrelevant.
const double cellCorrMin = 0.5;
const double cellCorrMax = 2.0;
const double wCorrMin = 0.5;
const double wCorrMax = 2.0;

// FarDet settings, uncomment as appropriate

//int cellMax = 32*12;
//double cutW = 750;
//const double cellCorrMin = 0.5;//2.0;
//const double cellCorrMax = 2.0;//2.6;
//const double wCorrMin = 0.5;//1.8;
//const double wCorrMax = 2.0;//3.0;

TH1F* singlecellcorr(TH2* h2, int bin, TString name)
{

  int nbinsx = h2->GetXaxis()->GetNbins();
  int minx = h2->GetXaxis()->GetXmin();
  int maxx = h2->GetXaxis()->GetXmax();
  TH1F* h1 = new TH1F (name, name+";W (cm);", nbinsx, minx, maxx);//h2->GetXaxis()->GetXmin(), h2->GetXaxis()->GetXmax);

  for(int i=0;i<nbinsx;i++)
    {
      h1->SetBinContent(i, h2->GetBinContent(i, bin));
    }
  return h1;

}

TVectorD solve(TH2* h2)
{
    // Linearize three dimensions into one
    TMatrixD A(Na*Nb*Nc+Nflag, Na*Nb*Nc+Nflag);
    // The A matrix has a symmetry property. Build up this vector instead
    TVectorD Avec(Na*Nb*Nc*8);
    TVectorD Nvec(Na*Nb*Nc);
    TMatrixD Aflag(Na*Nb*Nc+Nflag, Nflag);
    TVectorD vec(Na*Nb*Nc+Nflag);

    // Precompute values of pow(w/1000, a+d). 1000 is to scale to approximately
    // 0-1 so that matrix entries will all be comparably sized.
    double wads[h2->GetNbinsX()+2][Na*2];
    for(int iw = 0; iw < h2->GetNbinsX()+2; ++iw)
        for(int ad = 0; ad < Na*2; ++ad)
            wads[iw][ad] = pow(h2->GetXaxis()->GetBinCenter(iw)/1000, ad);

    // Precompute values of pow(cell/1000., b+e)
    double cbes[32*12][Nb*2];
    for(int cell = 0; cell < 32*12; ++cell)
        for(int be = 0; be < Nb*2; ++be)
            cbes[cell][be] = pow(cell/1000., be);

    // Precompute values of pow(mode, c+f)
    double mcfs[32][Nc*2];
    for(int mod = 0; mod < 32; ++mod)
        for(int cf = 0; cf < Nc*2; ++cf)
            mcfs[mod][cf] = pow(mod, cf);

    for(int iw = 0; iw < h2->GetNbinsX()+2; ++iw){
        //      std::cout << iw << std::endl;
        //      const double w = h2->GetXaxis()->GetBinCenter(iw);
        for(int ic = 0; ic < h2->GetNbinsY()+2; ++ic){
            // One for underflow, and one because we added one earlier to put a white
            // space around
            const int cell = ic-2;
            const int mod = cell%32;

            // I don't understand why start and end behave differently
            // physically. But empirically they do.
            int isStart = int(mod == 0);
            int isEnd = int(mod == 31);
            int isMid = int(mod == 15 || mod == 16);

            const double vi = h2->GetBinContent(iw, ic);
            // Don't fit to zeros
            if(!vi) continue;

            for(int a = 0; a < Na*2; ++a){
                const double wa = wads[iw][a]; //pow(w/1000, a);
                for(int b = 0; b < Nb*2; ++b){
                    const double cb = cbes[cell][b]; //pow(cell/1000., b);
                    for(int c = 0; c < Nc*2; ++c){
                        const double mc = mcfs[mod][c]; //pow(mod/10., c);

                        Avec[Nb*Nc*4*a+Nc*2*b+c] += wa*cb*mc;

                        if(a < Na && b < Nb && c < Nc){
                            vec[Nb*Nc*a+Nc*b+c] += vi*wa*cb*mc;

                            if(isStart) Aflag[Nb*Nc*a+Nc*b+c][0] += wa*cb*mc;
                            if(isEnd)   Aflag[Nb*Nc*a+Nc*b+c][1] += wa*cb*mc;
                            if(isMid)   Aflag[Nb*Nc*a+Nc*b+c][2] += wa*cb*mc;
                        }
                    } // end for c
                } // end for b
            } // end for a

            // Only need to do diagonals. Can't have both flags at once
            if(isStart) Aflag[Na*Nb*Nc+0][0] += 1;
            if(isEnd)   Aflag[Na*Nb*Nc+1][1] += 1;
            if(isMid)   Aflag[Na*Nb*Nc+2][2] += 1;

            vec[Na*Nb*Nc+0] += vi*isStart;
            vec[Na*Nb*Nc+1] += vi*isEnd;
            vec[Na*Nb*Nc+2] += vi*isMid;
        } // end for ic
    } // end for iw

    // Expand Avec out into the cells of A
    for(int a = 0; a < Na; ++a){
        for(int b = 0; b < Nb; ++b){
            for(int c = 0; c < Nc; ++c){
                for(int d = 0; d < Na; ++d){
                    for(int e = 0; e < Nb; ++e){
                        for(int f = 0; f < Nc; ++f){
                            A[Nb*Nc*d+Nc*e+f][Nb*Nc*a+Nc*b+c] = Avec[Nb*Nc*4*(a+d)+Nc*2*(b+e)+c+f];
                        }
                    }
                }
            }
        }
    }

    for(int flag = 0; flag < Nflag; ++flag){
        for(int abc = 0; abc < Na*Nb*Nc+Nflag; ++abc){
            const double val = Aflag[abc][flag];
            A[Na*Nb*Nc+flag][abc] = val;
            A[abc][Na*Nb*Nc+flag] = val;
        }
    }

    //    A.Print();
    //    vec.Print();

    TDecompSVD dc(A);
    bool junk;
    TVectorD p = dc.Solve(vec, junk);

    return p;
}

void fit_thresh_corrs()
{
  //std::cout << "SHOULD HAVE CHANGES" << std::endl;
    // Output from ThresholdAna_module
    
    // TFile* fin = new TFile("fdThresh_1.root");
    // TFile* fin = new TFile("/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/ndmc/Ideal_Conditions/threshold_2ndAna_fulldetector_ndmc_idealconditions.root");  // threshold and shadowing file. Output file of ThresholdAna_module.cc
    // TFile* fin = new TFile("/nova/ana/users/mcampbel/calibration_miniproduction/ndmc/consolidated/threshana_hadded.root"); 
    // TFile* fin = new TFile("/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/ndmc/Ideal_Conditions/threshold_2ndAna_fulldetector_ndmc_idealconditions.root");
    TFile* fin = new TFile("/pnfs/nova/persistent/users/kmulder/testbeam_rel_calib_28062019/thresholdCorr_mc/thresholdAna_mc_hadded.root");
    TTree* tr = new TTree("tree", "tree");
    int nViews = 2;//CHANGE TO 2 FOR FD, 4 FOR ND
    int view;
    tr->Branch("view", &view);
    int wpow, cellpow, modpow, startpow, midpow, endpow;
    double coeff;
    tr->Branch("wpow", &wpow);
    tr->Branch("cellpow", &cellpow);
    tr->Branch("modpow", &modpow);
    tr->Branch("startpow", &startpow);
    tr->Branch("midpow", &midpow);
    tr->Branch("endpow", &endpow);
    tr->Branch("coeff", &coeff);

    TFile* f_correction             = new TFile("threshold_comparison.root", "RECREATE"); // this file will be used to compare threshold and shielding correction factors.
    //TFile* fout = new TFile("thresh_fits.root", "RECREATE");

    double maxW = -1;
    const TString viewStrArr[4]={"xview","yview","muxview","muyview"};
    const TString viewStrUserArr[4]={"X View","Y View","Muon Catcher X View","Muon Catcher Y View"};
    for(view = 0; view < nViews; ++view){
        const TString viewStr = viewStrArr[view];
        const TString viewStrUser = viewStrUserArr[view];

        TH2* h2 = 0; // Wait until maxW is set

        for(int cell = 0; cell < cellMax; ++cell){
            TH2* h = (TH2*)fin->Get(TString::Format("thresh/ratioVsW_%s_cell%d", viewStr.Data(), cell));
            if(!h) continue;

            maxW = h->GetXaxis()->GetXmax();

            if(!h2) h2 = new TH2F("", viewStrUser+";W;Cell", 100, -maxW, +maxW, cellMax+2, -1, cellMax+1);

            TProfile* pr = h->ProfileX();

            for(int n = 0; n < pr->GetNbinsX(); ++n){
                // Weird stuff happens right at the edges. Don't include it in the fit,
                // because it compromises the main body.
                if(fabs(h2->GetXaxis()->GetBinCenter(n)) > cutW && view!=2) continue;
                if(fabs(h2->GetXaxis()->GetBinCenter(n)) > cutWMuonX && view==2) continue;

                // +2, one for overflow, one to put a white cell padding around
                h2->SetBinContent(n, cell+2, pr->GetBinContent(n));
            }
        }


        TH2* l2 = 0; 
        for(int cell = 0; cell < cellMax; ++cell){
          TH2* l = (TH2*)fin->Get(TString::Format("thresh/fPEpercmVsW_%s_cell%d", viewStr.Data(), cell));
          if(!l) continue;

          maxW = l->GetXaxis()->GetXmax();

          if(!l2) l2 = new TH2F("", viewStrUser+";W;Cell", 100, -maxW, +maxW, cellMax+2, -1, cellMax+1);

          TProfile* prl = l->ProfileX();
          for(int n = 0; n < prl->GetNbinsX(); ++n){
            // Weird stuff happens right at the edges. Don't include it in the fit,
            // because it compromises the main body.
            if(fabs(l2->GetXaxis()->GetBinCenter(n)) > cutW && view!=2) continue;
            if(fabs(l2->GetXaxis()->GetBinCenter(n)) > cutWMuonX && view==2) continue;

            // +2, one for overflow, one to put a white cell padding around
            l2->SetBinContent(n, cell+2, prl->GetBinContent(n));
          }
        }
        new TCanvas;
        l2->GetXaxis()->CenterTitle();
        l2->GetYaxis()->CenterTitle();
        l2->Draw("colz");
        // Don't need to see zeros
        //l2->SetMinimum(.5);
        //l2->SetMaximum(3.0);
        gPad->Print("plots/thresh_PEpercm_"+viewStr+".eps");


        new TCanvas;
        h2->GetXaxis()->CenterTitle();
        h2->GetYaxis()->CenterTitle();
        h2->Draw("colz");
        // Don't need to see zeros
        h2->SetMinimum(.5);
        h2->SetMaximum(3.0);
        gPad->Print("plots/thresh_corr_"+viewStr+".eps");

        new TCanvas;
        TString newplotname = "Cell31";
        TH1F* cell31plot = singlecellcorr(h2, 31, newplotname);
        cell31plot->Draw();
        gPad->Print("plots/thresh_corr_cell31"+viewStr+".eps");

        const TVectorD p = solve(h2);


        // Save coefficients to tree
        startpow = 0;
        midpow = 0;
        endpow = 0;
        for(int a = 0; a < Na; ++a){
            for(int b = 0; b < Nb; ++b){
                for(int c = 0; c < Nc; ++c){
                    wpow = a;
                    cellpow = b;
                    modpow = c;
                    coeff = p[Nb*Nc*a+Nc*b+c];
                    tr->Fill();
                }
            }
        }
        wpow = 0;
        cellpow = 0;
        modpow = 0;
        startpow = 1;
        coeff = p[Na*Nb*Nc+0];
        tr->Fill();
        startpow = 0;
        endpow = 1;
        coeff = p[Na*Nb*Nc+1];
        tr->Fill();
        endpow = 0;
        midpow = 1;
        coeff = p[Na*Nb*Nc+2];
        tr->Fill();

        // Make a map of what the fit looks like
        TH2* h2fit = new TH2F("", viewStrUser+";W;Cell", 100, -maxW, +maxW, cellMax+2, -1, cellMax+1);
        for(int iw = 0; iw < h2->GetNbinsX()+2; ++iw){
            const double w = h2->GetXaxis()->GetBinCenter(iw);
            for(int ic = 0; ic < h2->GetNbinsY()+2; ++ic){
                const int cell = ic-2;
                const int mod = cell%32;
                int isStart = int(mod == 0);
                int isEnd = int(mod == 31);
                int isMid = int(mod == 15 || mod == 16);

                double v = 0;
                for(int a = 0; a < Na; ++a){
                    for(int b = 0; b < Nb; ++b){
                        for(int c = 0; c < Nc; ++c){
                            v += p[Nb*Nc*a+Nc*b+c]*pow(w/1000, a)*pow(cell/1000., b)*pow(mod, c);
                        }
                    }
                }
                if(isStart) v += p[Na*Nb*Nc+0];
                if(isEnd) v += p[Na*Nb*Nc+1];
                if(isMid) v += p[Na*Nb*Nc+2];
                if(h2->GetBinContent(iw, ic)) h2fit->SetBinContent(iw, ic, v);
            }
        }

        new TCanvas;
        h2fit->GetXaxis()->CenterTitle();
        h2fit->GetYaxis()->CenterTitle();
        h2fit->Draw("colz");
        h2fit->SetMinimum(h2->GetMinimum());
        h2fit->SetMaximum(h2->GetMaximum());
        gPad->Print("plots/thresh_corr_fit_"+viewStr+".eps");

        new TCanvas;
        TH1* fx = h2fit->ProjectionX("a"+viewStr);
        int nbinsy = 0;
        for(int i = 0; i < h2->GetNbinsY()+2; ++i){
            if(h2fit->GetBinContent(h2->GetNbinsX()/2+1, i)) ++nbinsy;
        }
        fx->Scale(1./nbinsy);
        fx->SetLineColor(kRed);
        fx->GetYaxis()->SetTitle("Correction factor");
        fx->GetYaxis()->CenterTitle();
        fx->SetTitle(viewStrUser);
        fx->GetYaxis()->SetRangeUser(wCorrMin, wCorrMax);//1.1, 1.3
        fx->Draw("l");
        TH1* dx = h2->ProjectionX("b"+viewStr);
        dx->Scale(1./nbinsy);
        dx->SetLineColor(kBlack);
        dx->Draw("same");
        
        gPad->Print("plots/thresh_corr_proj_W_"+viewStr+".eps");
        fx->SetName("a"+viewStr);
        fx->Write();

        new TCanvas;
        TH1* fy = h2fit->ProjectionY("c"+viewStr);
        int nbinsx = 0;
        for(int i = 0; i < h2->GetNbinsX()+2; ++i){
          std::cout << "xbin : " << i << ", ybins/2 +1 : " << h2->GetNbinsY()/2+1 << ", content: " << h2fit->GetBinContent(i, h2->GetNbinsY()/2+1) << std::endl;
            if(h2fit->GetBinContent(i, h2->GetNbinsY()/2+1)) ++nbinsx;
        }
        std::cout << nbinsx << " nbinsx" << std::endl;
        fy->Scale(1./nbinsx);
        fy->SetLineColor(kRed);
        fy->GetYaxis()->SetTitle("Correction factor");
        fy->GetYaxis()->CenterTitle();
        fy->SetTitle(viewStrUser);
        fy->GetYaxis()->SetRangeUser(cellCorrMin, cellCorrMax);//1.1, 1.25
        fy->Draw("l");
        TH1* dy = h2->ProjectionY("d"+viewStr);
        dy->Scale(1./nbinsx);
        dy->SetLineColor(kBlack);
        dy->Draw("same");

        gPad->Print("plots/thresh_corr_proj_cell_"+viewStr+".eps");
        fy->SetName("c"+viewStr);
        fy->Write();

        new TCanvas;
        TH2* ratio = (TH2*)h2->Clone();
        ratio->Divide(h2fit);
        ratio->Draw("colz");
        ratio->SetMinimum(.97);
        ratio->SetMaximum(1.03);
        gPad->Print("plots/thresh_corr_ratio_"+viewStr+".eps");

        new TCanvas;
        TH1* px = ratio->ProjectionX(viewStr+"_px");
        px->Scale(1./nbinsy);
        px->GetYaxis()->SetRangeUser(.98, 1.02);
        px->Draw();
        gPad->Print("plots/thresh_corr_ratio_W_"+viewStr+".eps");

        new TCanvas;
        TH1* py = ratio->ProjectionY(viewStr+"_py");
        py->Scale(1./nbinsx);
        py->GetYaxis()->SetRangeUser(.98, 1.02);
        py->Draw();
        gPad->Print("plots/thresh_corr_ratio_cell_"+viewStr+".eps");
    } // end for view
//    f_correction->Close();

    TFile* fout = new TFile("thresh_fits.root", "RECREATE");
    tr->SetDirectory(fout);
    tr->Write();
    fout->Close();

    std::cout << "Wrote results to thresh_fits.root" << std::endl;
}