epic-lfhcal-tbana/OldStructure/MultiMultiphoton.C

0001 #include "TFile.h"
0002 #include "TTree.h"
0003 #include "TH1D.h"
0004 #include "TSpectrum.h"
0005 #include "TCanvas.h"
0006 #include "TLatex.h"
0007 #include "TMath.h"
0008 #include "TLeaf.h"
0009 #include "TF1.h"
0010 #include "TGraph.h"
0011 #include "TMultiGraph.h"
0012 #include "TVirtualFitter.h"
0013 #include <iostream>
0014
0015 typedef double Double_t;
0016
0017 using namespace std;
0018 using std::ofstream;
0019
0020 Int_t npeaks = 0;
0021
0022 const int gNumSipm = 4;
0023 const int gNumOv = 1;
0024
0025 Double_t fpeaks(Double_t * x, Double_t * par)
0026 {
0027   Double_t result = par[0] + par[1] * x[0];
0028   for (Int_t p = 0; p < npeaks; p++)
0029     {
0030       Double_t norm = par[3 * p +2]; // "height" or "area"
0031       Double_t mean = par[3 * p + 3];
0032       Double_t sigma = par[3 * p + 4];
0033 #if defined(__PEAKS_C_FIT_AREAS__)
0034       norm /= sigma * (TMath::Sqrt(TMath::TwoPi())); // "area"
0035 #endif                                               /* defined(__PEAKS_C_FIT_AREAS__) */
0036       result += norm * TMath::Gaus(x[0], mean, sigma);
0037     }
0038   return result;
0039 }
0040
0041
0042 void MultiMultiphoton(){
0043
0044     const int ch[4] = {00, 01, 32, 33};
0045     TH1D *h[4];
0046     TH1 *bg[4];
0047     TSpectrum *s[4];
0048
0049     TCanvas *c = new TCanvas("c", "c",515,81,837,763);
0050     c->Divide(4, 2); // The pad is 1 coloumn and two lines
0051
0052     std::string filenames[gNumOv];
0053     filenames[0] = "Truncated_Dummy.txt";
0054 //    filenames[1] = "/Users/ar2545/Downloads/S14160_1315_Scans/Trunc_Run1_list_45.txt";
0055 //    filenames[2] = "/Users/ar2545/Downloads/S14160_1315_Scans/Trunc_Run1_list_44.txt";
0056     double Vop[gNumOv];
0057     double yAvg[gNumSipm][gNumOv];
0058
0059     for(int i=0; i < gNumOv; i++)
0060     {
0061         //        if(i>0)break;
0062         TTree *t = new TTree(Form("ntuple_%d",i),"data from ascii file"); // A temporary tree
0063
0064         t->ReadFile(filenames[i].c_str(), "a:b:c"); // Read the data file
0065
0066
0067         for (Int_t k = 0; k < gNumSipm; k++) {
0068
0069             h[k] = new TH1D(Form("h_%d_%d",k,i), "Multiphoton Spectrum; ADC ; Counts", t->GetEntries(), 0., (Double_t)t->GetEntries()); // A new histogram
0070             h[k]->Sumw2(kFALSE);
0071             h[k]->ResetStats();
0072
0073             // Fill histogram with values from column c where column a is equal to 0
0074             for (Int_t i = 0; i < t->GetEntries(); ++i) {
0075                 t->GetEntry(i);
0076                 if (t->GetLeaf("a")->GetValue() == ch[k]) {
0077                     //                      cout<< "check 1 ======>"<<ch[k]<< endl;
0078                     //                    cout<< "check 2 ======>" << t->GetLeaf("c")->GetValue()<<endl;
0079
0080
0081                     h[k]->Fill(t->GetLeaf("c")->GetValue());
0082                 }
0083             }
0084             cout<< "check 3 ---------->"<<endl;
0085
0086             //delete t; // No longer needed
0087
0088             s[k] = new TSpectrum(); // Use TSpectrum to find the peak candidates
0089
0090             bg[k] = s[k]->Background(h[k], 20, "smoothing"); // Estimate the background
0091
0092             bg[k]->SetLineColor(kRed); // Set the color
0093
0094             c->cd(k+5);        // Move to line one
0095             bg[k]->Draw("SAME"); // Draw the bg in the first line
0096
0097             //h->SetTitleFont(82);
0098
0099             //h->GetYaxis()->SetTitleFont(82);
0100             h[k]->GetYaxis()->SetTitleSize(0.05);
0101             h[k]->GetXaxis()->SetTitleSize(0.05);
0102             h[k]->GetYaxis()->SetTitleOffset(1);
0103             h[k]->GetXaxis()->SetTitleOffset(1);
0104             h[k]->GetXaxis()->SetLabelSize(0.05);
0105             h[k]->GetYaxis()->SetLabelSize(0.05);
0106             h[k]->GetXaxis()->SetRange(0,2000);
0107             h[k]->SetStats(0);
0108             h[k]->SetTitle("S1416-1315 @  V_{op} = 7 V (High-Gain)");
0109             h[k]->DrawCopy("");
0110             TLatex *   tex1 = new TLatex();
0111
0112 //            bg[k]->Draw("SAME"); // Draw the bg in the first line
0113
0114             //if (gPad)
0115             //gPad->SetGrid(1, 1);
0116
0117             c->cd(k+1); // Move to line two
0118
0119             h[k]->Add(bg[k], -1.); // Subtract the estimated background
0120
0121             Int_t nfound = s[k]->Search(h[k], 7., "", 0.05); // Search for peaks
0122
0123             printf("We found %d candidate peaks to fit.\n", nfound);
0124
0125             Double_t *xpeaks;
0126             Double_t par[3000];
0127             xpeaks = s[k]->GetPositionX();
0128
0129             std::sort(xpeaks, xpeaks + nfound, [&](Double_t a, Double_t b) {
0130                 return a > b;
0131             });
0132
0133             std::cout << "Differences between consecutive entries:" << std::endl;
0134             Double_t sum_diff = 0.0;
0135             for (size_t i = 1; i < (Int_t)nfound; ++i) {
0136                 Double_t diff = xpeaks[i - 1] - xpeaks[i];
0137                 std::cout << diff << " ";
0138                 sum_diff += diff;
0139             }
0140             std::cout << std::endl;
0141
0142             Double_t average_diff = sum_diff / (nfound - 1);
0143             std::cout << "Average of differences: " << average_diff << std::endl;
0144             yAvg[k][i] = average_diff;
0145
0146             for (Int_t p = 0; p < (Int_t)nfound; p++)
0147             {
0148                 Double_t xp = xpeaks[p];
0149                 Int_t bin = h[k]->GetXaxis()->FindBin(xp);
0150                 Double_t yp = h[k]->GetBinContent(bin);
0151
0152                 par[3 * npeaks + 2] = yp;     // "height"
0153                 par[3 * npeaks + 3] = xp; // "mean"
0154                 par[3 * npeaks + 4] = 3;  // "sigma"
0155
0156
0157                 cout<<xp<<endl;
0158                 npeaks++;
0159
0160             }
0161
0162
0163             printf("Howdy, now fitting: Be patient.\n");
0164
0165             TF1 *fit = new TF1("fit", fpeaks, 0, 4096, 3 * npeaks);
0166
0167             TVirtualFitter::Fitter(h[k], 3 * npeaks);
0168             fit->SetParameters(par);
0169             fit->SetNpx(1000);
0170
0171             h[k]->Fit("fit");
0172
0173             //h->GetYaxis()->SetTitleFont(82);
0174             //h->GetXaxis()->SetTitleFont(82);
0175
0176             h[k]->Draw("E"); // Draw the plot without bg
0177
0178
0179             //if (gPad)
0180             //gPad->SetGrid(1, 1);
0181
0182             /* c->cd(2); */
0183
0184             /* TLatex *   tex = new TLatex(); */
0185             /* tex->DrawLatex(1400,2225.444,"Background Subtracted"); */
0186             /* tex->DrawLatex(1030.066,1693.444,Form("Average peak-to-peak value = %0.2f",average_diff)); */
0187             /* tex->SetTextFont(82); */
0188             /* tex->Draw("SAME"); */
0189
0190             printf("There, done! You happy?\n");
0191
0192
0193         }
0194     }
0195     c->SaveAs("fittedPeak.pdf");
0196 //    cout<<yAvg[0][0]<<endl;
0197 //    cout<<yAvg[3][2]<<endl;
0198
0199     double y0[gNumOv];
0200 //    double y1[gNumOv];
0201 //    double y2[gNumOv];
0202 //    double y3[gNumOv];
0203
0204     for (int i = 0; i < gNumOv; ++i) {
0205         y0[i] = yAvg[0][i];
0206         cout<<y0[i]<<endl;
0207 //        y1[i] = yAvg[1][i];
0208 //        cout<<y1[i]<<endl;
0209 //        y2[i] = yAvg[2][i];
0210 //        cout<<y2[i]<<endl;
0211 //        y3[i] = yAvg[3][i];
0212
0213     }
0214
0215     TGraph* g = new TGraph(gNumOv,y0,Vop);
0216 //    TGraph* g1 = new TGraph(gNumOv,y1,Vop);
0217 //    TGraph* g2 = new TGraph(gNumOv,y2,Vop);
0218 //    TGraph* g3 = new TGraph(gNumOv,y3,Vop);
0219 //
0220     TMultiGraph *multiGraph = new TMultiGraph();
0221     multiGraph->Add(g);
0222 //    multiGraph->Add(g1);
0223
0224
0225     TCanvas *canvas = new TCanvas("canvas", "Multiple Graphs", 800, 600);
0226     multiGraph->Draw("APL"); // "APL" stands for axes, points, and line
0227 }