epic-lfhcal-tbana/OldStructure/Gauss.C

0001 #include <TCanvas.h>
0002 #include <TH1F.h>
0003 #include <TF1.h>
0004 #include <TRandom3.h>
0005
0006 const int gMaxChannels = 64;
0007 const int gMaxBoard    = 1;
0008 // track hits
0009 Long64_t gTrID;
0010 Double_t gTRtimeStamp;
0011 Long64_t* gBoard          = new Long64_t[gMaxChannels];
0012 Long64_t* gChannel        = new Long64_t[gMaxChannels];
0013 Long64_t* gLG             = new Long64_t[gMaxChannels];
0014 Long64_t* gHG             = new Long64_t[gMaxChannels];
0015 void SetBranchAddressesTree(TTree* inputTree){
0016     //**********************************
0017     //format
0018     //**********************************
0019
0020     //     *Tree    :tree      :                                                        *
0021     //     ******************************************************************************
0022     //     *Br    0 :t_stamp   : t_stamp/D                                              *
0023     //     *Br    1 :trgid     : trgid/L                                                *
0024     //     *Br    2 :board     : board[64]/L                                            *
0025     //     *Br    3 :channel   : channel[64]/L                                          *
0026     //     *Br    4 :LG        : LG[64]/L                                               *
0027     //     *Br    5 :HG        : HG[64]/L                                               *
0028     //******************************************************************************
0029
0030
0031     if (inputTree->GetBranchStatus("t_stamp") ){
0032       inputTree->SetBranchAddress("trgid",          &gTrID);
0033       inputTree->SetBranchAddress("t_stamp",          &gTRtimeStamp);
0034       inputTree->SetBranchAddress("board",            gBoard);
0035       inputTree->SetBranchAddress("channel",          gChannel);
0036       inputTree->SetBranchAddress("LG",               gLG);
0037       inputTree->SetBranchAddress("HG",               gHG);
0038     }
0039 }
0040
0041 void Gauss(const char* noisefile, const char* datafile, int runnumber, TString outputDir    = "NoiseSubData/MIP/")
0042 {
0043
0044
0045     outputDir = Form("%s/Run%05d",outputDir.Data(), runnumber );
0046     gSystem->Exec("mkdir -p "+outputDir);
0047     // Create a ROOT canvas to display the histogram and fit results
0048     TCanvas *c1 = new TCanvas("c1", "Gaussian Fit Example", 800, 600);
0049
0050     //When you dont have a noise sweep, you can fit the noise up to 40ADC. So the noise and data file are the same in this case except you still have to subtract at the tree level
0051     //This is the processed noise data.
0052     TFile* fnoise = new TFile(noisefile,"READ");
0053     //This is th unprocessed data file
0054     TFile* fdata = new TFile(datafile,"READ");
0055
0056     //These are noise histograms that exist in processed data file
0057     TH1F* nhistLG[gMaxBoard][gMaxChannels];//Histogram array
0058     TH1F* nhistHG[gMaxBoard][gMaxChannels];
0059
0060     TH1F* histLG[gMaxBoard][gMaxChannels];
0061     TH1F* histHG[gMaxBoard][gMaxChannels];
0062
0063     //noise fit histograms
0064     TH1F* nfithistLG[gMaxBoard][gMaxChannels];
0065     TH1F* nfithistHG[gMaxBoard][gMaxChannels];
0066
0067     //signal histogram
0068     TH1F* shistLG[gMaxBoard][gMaxChannels];
0069     TH1F* shistHG[gMaxBoard][gMaxChannels];
0070
0071     const int rows = gMaxBoard; // Number of rows
0072     const int cols = gMaxChannels; // Number of columns
0073     std::vector<std::vector<double>> gauFitSigma(rows, std::vector<double>(cols));
0074     std::vector<std::vector<double>> gauFitMean(rows, std::vector<double>(cols));
0075
0076     std::vector<std::vector<double>> gauFitSigmaLG(rows, std::vector<double>(cols));
0077     std::vector<std::vector<double>> gauFitMeanLG(rows, std::vector<double>(cols));
0078
0079     TF1* gauFit[gMaxBoard][gMaxChannels];
0080
0081
0082     //Create histograms noise, data and signal
0083     for (Int_t j = 0; j < gMaxBoard; j++)
0084     {
0085         for (Int_t i = 0; i < gMaxChannels; i++)
0086         {
0087             // Define the Gaussian function. Will define custom parameters after the fit.
0088             nfithistHG[j][i] = new TH1F(Form("hfit_HG_B%d_C%02d",j,i),"",4001, 0, 4001);
0089             nfithistLG[j][i] = new TH1F(Form("hfit_LG_B%d_C%02d",j,i),"",4001, 0, 4001);
0090             //Data histogram
0091             histHG[j][i] = new TH1F(Form("HG_TriggB%d_C%02d",j,i),"",4001,0,4001);//changing the name for each hist and adding to histogram list
0092             histHG[j][i]->RebinX(2);
0093             histLG[j][i] = new TH1F(Form("LG_B%d_C%02d",j,i),"",4001,0,4001);
0094             //Subtracted histogram
0095             shistHG[j][i] = new TH1F(Form("sh_HG_B%d_C%02d",j,i),"",4001,0,4001);//changing the name for each hist and adding to histogram list
0096             shistHG[j][i]->RebinX(2);
0097             shistLG[j][i] = new TH1F(Form("sh_LG_B%d_C%02d",j,i),"",4001,0,4001);
0098
0099         }
0100     }
0101
0102      TRandom3 randomGenerator;
0103     //Now fit the noise data and create fit histograms
0104     for (Int_t j = 0; j < gMaxBoard; j++){
0105         for (Int_t i = 0; i < gMaxChannels; i++){
0106             //noise histogram
0107             nhistHG[j][i] = (TH1F*)fnoise->Get(Form("h_HGNoise_B%d_C%02d",j,i));
0108             nhistHG[j][i]->RebinX(2);
0109 //            nhistHG[j][i]->Scale(1./nhistHG[j][i]->Integral());
0110
0111             cout<<nhistHG[j][i]->GetName()<<endl;
0112             cout<<nhistHG[j][i]->GetEntries()<<endl;
0113             cout<<nhistHG[j][i]->GetMaximum()<<endl;
0114
0115             TF1 *gaussianFit = new TF1("gaussianFit", "gaus");
0116
0117 //
0118             gaussianFit->SetParameter(0, nhistHG[j][i]->GetMaximum()); // Amplitude
0119             gaussianFit->SetParameter(1, nhistHG[j][i]->GetMean());   // Mean
0120             gaussianFit->SetParameter(2, nhistHG[j][i]->GetStdDev());   // Sigma
0121             nhistHG[j][i]->Fit("gaussianFit", "R","",0,100); //fit a gaussian to the noise histogram
0122
0123
0124             double_t amplitude = gaussianFit->GetParameter(0);
0125             double_t mean = gaussianFit->GetParameter(1);
0126             double_t sigma = gaussianFit->GetParameter(2);
0127
0128
0129            gauFitMean[j][i] = mean;
0130            gauFitSigma[j][i] = sigma;
0131            gauFit[j][i]= gaussianFit;
0132
0133
0134         }
0135     }
0136
0137     //Now open the data file, read it and subtract noise from signal on the fly
0138     //For the actual signal we have raw data
0139
0140     TTree *tt_event = (TTree*)fdata->Get("tree");
0141     SetBranchAddressesTree(tt_event);
0142     double nEntriesTree = tt_event->GetEntries();
0143     Long64_t startEvent = 0;
0144     Long64_t nEventsProcessed=0;
0145     for (Long64_t i=startEvent; i<nEntriesTree;i++)
0146     {
0147         // load current event
0148         tt_event->GetEntry(i);
0149         nEventsProcessed++;
0150         for(Int_t ch=0; ch<gMaxChannels; ch++)
0151         {
0152
0153             histHG[gBoard[ch]][gChannel[ch]]->Fill(gHG[ch]);
0154
0155             shistHG[gBoard[ch]][gChannel[ch]]->Fill((gHG[ch]-randomGenerator.Gaus(gauFitMean[gBoard[ch]][gChannel[ch]],gauFitSigma[gBoard[ch]][gChannel[ch]])));
0156
0157 //            cout<<gauFit[gBoard[ch]][gChannel[ch]]->GetRandom()<<endl;
0158
0159         }
0160 //
0161     }
0162
0163
0164     //    histHG[0][24]->SetLineColor(kBlue);
0165     //    histHG[0][24]->GetXaxis()->SetRangeUser(0,500);
0166     //    histHG[0][24]->Draw();
0167     TLegend *legend_dat = new TLegend(0.60,0.60,0.80,0.90);
0168     legend_dat->AddEntry(nhistHG[0][24],"Noise","l");
0169     legend_dat->AddEntry(shistHG[0][24],"Noise subtracted Triggered","l");
0170     legend_dat->SetBorderSize(0);
0171     nhistHG[0][24]->SetLineColor(kBlue);
0172     nhistHG[0][24]->GetXaxis()->SetTitle("ADC");
0173     nhistHG[0][24]->GetYaxis()->SetTitle("Entries");
0174     nhistHG[0][24]->Draw();
0175     shistHG[0][24]->SetLineColor(kGreen+4);
0176     shistHG[0][24]->Draw("SAME");
0177     legend_dat->Draw("SAME");
0178
0179 //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
0180     //Write this to a new file
0181     //Get_run_name from data file
0182
0183     TFile* fileOutput = new TFile(Form("%s/%s",outputDir.Data(),datafile),"RECREATE");
0184     for (Int_t j = 0; j < gMaxBoard; j++)
0185     {
0186         for (Int_t i = 0; i < gMaxChannels; i++)
0187         {
0188             shistHG[j][i]->Write();
0189             shistLG[j][i]->Write();
0190
0191         }
0192     }
0193
0194 }
0195
0196
0197