/************************************************************** * * 2015/07/08 Kei Moriya * * Program that will take in trees created by TPOL_waveform * plugin and * - make plots of waveforms * - make average subtraction of pedestal from waveforms * - make trees for further processing * **************************************************************/ #include #include #include #include #include #include #include #include "unistd.h" // to use optarg #include // ROOT header files #include "TROOT.h" #include "TFile.h" #include "TH1F.h" #include "TH2F.h" #include "TF1.h" #include "TProfile.h" #include "TTree.h" #include "TMath.h" #include "TLatex.h" #include "TLine.h" #include "TLegend.h" #include "TPaveText.h" #include "TCanvas.h" #include "TPad.h" #include "TLorentzVector.h" #include "TProfile.h" #include "TGraphErrors.h" #include "TStyle.h" // Header file that includes fit function to signal, // calibration for each channel, sectorToPhi, // utilities to set and get Tree branches #include "tpol.h" using namespace std; void helpMessage(void); int main(int argc, char **argv){ gStyle->SetOptStat(0); const Int_t colors[] = {kBlack, kRed, kYellow+2, kGreen+2, kBlue, kMagenta, kCyan, 32}; /////////////////////////////////////////////////////////////// // // // Set up configuration parameters // // // /////////////////////////////////////////////////////////////// Int_t MAXEVENTS = 100 * 1000 * 1000; std::string infilename = ""; char outfilename[400]; Bool_t outfile_specified = false; sprintf(outfilename,"out_plot_waveform.root"); // default Bool_t verbose = false; // Number of fitted waveforms to keep in output file. Int_t NMAXKEEP = 50; // Difference in max fADC counts in waveform and pedestal // to draw and save waveforms for an event. Double_t DRAW_THRESHOLD = 100; // Difference of maximum and minimum fADC counts in waveform // for event to be fit with signal function Double_t FIT_THRESHOLD = 30; // Difference of maximum and minimum fADC counts in waveform // for event to be counted in multiplicity studies Double_t MULTIPLICITY_THRESHOLD = 30; // Maximum difference of maxcounts and pedestal to be // used for average pedestal calculation Double_t AVERAGE_PEDESTAL_THRESHOLD = 30; // Range of pedestal histogram Double_t PEDESTAL_MAX = 800; // Range of maxcounts Double_t MAXCOUNTS_MAX = 3000; // Range of signal size Double_t SIGNAL_MAX = 300; // Readout threshold (depends on DAQ setting) Double_t READOUT_THRESHOLD = 1120; char *progName = argv[0]; extern char* optarg; // Check command line arguments int c; while((c = getopt(argc,argv,"hi:o:vM:N:t:d:m:p:q:r:s:T:")) != -1){ switch(c){ case 'h': helpMessage(); exit(-1); break; case 'i': // Specify infile name. infilename = optarg; cout << "infile name: " << infilename << endl; break; case 'o': // Specify outfile name sprintf(outfilename,optarg); outfile_specified = true; cout << "outfile: " << outfilename << endl; break; case 'M': // Specify max # of events to process MAXEVENTS = atoi(optarg); cout << "max # of events: " << MAXEVENTS << endl; break; case 'v': verbose = true; cout << "vebose mode" << endl; break; case 'N': // Maximum number of waveform histograms to save to output file NMAXKEEP = atoi(optarg); cout << "NMAXKEEP = " << NMAXKEEP << endl; break; case 't': // Threshold to save images for event DRAW_THRESHOLD = atof(optarg); cout << "draw threshold: " << DRAW_THRESHOLD << " ADC counts" << endl; break; case 'd': // Threshold to do signal fit FIT_THRESHOLD = atof(optarg); cout << "FIT_THRESHOLD = " << FIT_THRESHOLD << endl; break; case 'm': // Threshold to be counted in multiplicity MULTIPLICITY_THRESHOLD = atof(optarg); cout << "MULTIPLICITY_THRESHOLD = " << MULTIPLICITY_THRESHOLD << endl; break; case 'p': // Threshold to be counted in average waveform (no signal) AVERAGE_PEDESTAL_THRESHOLD = atof(optarg); cout << "AVERAGE_PEDESTAL_THRESHOLD = " << AVERAGE_PEDESTAL_THRESHOLD << endl; break; case 'q': // Maximum range of pedestal histogram PEDESTAL_MAX = atof(optarg); cout << "PEDESTAL_MAX = " << PEDESTAL_MAX << endl; break; case 'r': // Maximum range of maxcounts histogram MAXCOUNTS_MAX = atof(optarg); cout << "MAXCOUNTS_MAX = " << MAXCOUNTS_MAX << endl; break; case 's': // Maximum range of signal size SIGNAL_MAX = atof(optarg); cout << "SIGNAL_MAX = " << SIGNAL_MAX << endl; break; case 'T': // Readout threshold READOUT_THRESHOLD = atof(optarg); cout << "READOUT_THRESHOLD = " << READOUT_THRESHOLD << endl; break; default: break; } } if(infilename == ""){ helpMessage(); abort(); } //___________________________________________________________________________________________ /////////////////////////////////////////////////////////////// // // // Set up input file, tree // // // /////////////////////////////////////////////////////////////// TFile *infile = new TFile(infilename.c_str()); TTree *intree = (TTree*)infile->Get("TPOL_waveformTree"); const Int_t NMAX = 40; const Int_t NSECTORS = 32; // Value that pedestal will be adjusted to const Double_t PEDESTAL = 100.; // Utility variables to retrieve graph points Double_t x,y; // TPOL waveform variables UInt_t nwaveforms; Int_t eventnum; Int_t slot[NMAX]; Int_t sector[NMAX]; Double_t pedestal[NMAX]; UShort_t maxcounts[NMAX]; UShort_t waveform[NMAX][100]; // PSC stuff UInt_t nhits_PSC; Int_t iLeft_PSC[NMAX]; Int_t iRight_PSC[NMAX]; Float_t t_iLeft_PSC[NMAX]; Float_t t_iRight_PSC[NMAX]; // PS stuff UInt_t nhits_PS; Int_t iLeft_PS[NMAX]; Int_t iRight_PS[NMAX]; Float_t e_iLeft_PS[NMAX]; Float_t e_iRight_PS[NMAX]; // ADC hit info intree->SetBranchAddress("nwaveforms",&nwaveforms); intree->SetBranchAddress("eventnum",&eventnum); intree->SetBranchAddress("slot",slot); intree->SetBranchAddress("sector",sector); intree->SetBranchAddress("pedestal",pedestal); intree->SetBranchAddress("maxcounts",maxcounts); intree->SetBranchAddress("waveform",waveform); // PSC, PS stuff intree->SetBranchAddress("nhits_PSC",&nhits_PSC); intree->SetBranchAddress("iLeft_PSC",iLeft_PSC); intree->SetBranchAddress("iRight_PSC",iRight_PSC); intree->SetBranchAddress("t_iLeft_PSC",t_iLeft_PSC); intree->SetBranchAddress("t_iRight_PSC",t_iRight_PSC); intree->SetBranchAddress("nhits_PS",&nhits_PS); intree->SetBranchAddress("iLeft_PS",iLeft_PS); intree->SetBranchAddress("iRight_PS",iRight_PS); intree->SetBranchAddress("e_iLeft_PS",e_iLeft_PS); intree->SetBranchAddress("e_iRight_PS",e_iRight_PS); // Get file name only with no directory info // to set up outfile name char command[200]; sprintf(command,"basename %s > ___tmp.txt",infilename.c_str()); system(command); ifstream IN___tmp("___tmp.txt"); std::string basename; IN___tmp >> basename; system("rm -f ___tmp.txt"); // default outfile name if(outfile_specified == false){ sprintf(outfilename,"results/results.%s.root",basename.c_str()); } TFile *outfile = new TFile(outfilename,"recreate"); /////////////////////////////////////////////////////////////// // // // Set up fit function, output file, tree // // // /////////////////////////////////////////////////////////////// //Fit function (in tpol.h file) MySigFun * fptr = new MySigFun(); int nPar = 5; TF1 *sigFun = new TF1("sigFun",fptr,&MySigFun::Evaluate,0,1,nPar, "MySigFun","Evaluate"); TF1 *fsigFun_copy[NSECTORS]; // Save the fit parameters for all hits within an event // into a vector that is a pair of sector numbers and // fit parameters. std::vector > > vsector_fitparams; //Setup the output tree TTree *outTree = new TTree("tpolTree1","tpol1"); tpol1_t tpol1; // struct tpol1_t defined in tpol.h setBranchesT1(outTree,&tpol1); // setBranchesT1 defined in tpol.h /////////////////////////////////////////////////////////////// // // // Graphs for waveforms // // 1. raw waveforms // // 2. adjusted waveforms - pedestals moved to PEDESTAL // // 3. normalized waveforms - average pedestal is subtracted // // // /////////////////////////////////////////////////////////////// TGraph *gwaveform[NSECTORS]; TGraph *gwaveform_pedestalAdjusted[NSECTORS]; TGraph *gwaveform_normalized[NSECTORS]; TGraph *gwaveform_average = new TGraph(100); char gname[200]; // Format waveforms for drawing gwaveform_average->SetName("gwaveform_average"); gwaveform_average->SetMarkerColor(kCyan); gwaveform_average->SetMarkerStyle(20); gwaveform_average->SetMarkerSize(0.6); gwaveform_average->SetLineColor(kCyan); gwaveform_average->SetLineStyle(1); gwaveform_average->SetLineWidth(1); for(Int_t sec=0;secSetName(gname); gwaveform[sec]->SetTitle(""); gwaveform[sec]->GetXaxis()->SetTitle("sample #"); gwaveform[sec]->GetXaxis()->CenterTitle(); gwaveform[sec]->GetYaxis()->SetTitle("ADC counts"); gwaveform[sec]->GetYaxis()->CenterTitle(); gwaveform[sec]->SetMarkerColor(colors[sec%6]); gwaveform[sec]->SetMarkerStyle(20 + sec/6); gwaveform[sec]->SetMarkerSize(0.4); gwaveform[sec]->SetLineColor(colors[sec%6]); gwaveform[sec]->SetLineStyle(1); gwaveform[sec]->SetLineWidth(1); gwaveform_pedestalAdjusted[sec] = new TGraph(100); sprintf(gname,"gwaveform_pedestalAdjusted%2.2d",sec+1); gwaveform_pedestalAdjusted[sec]->SetName(gname); gwaveform_pedestalAdjusted[sec]->SetTitle(""); gwaveform_pedestalAdjusted[sec]->GetXaxis()->SetTitle("sample #"); gwaveform_pedestalAdjusted[sec]->GetXaxis()->CenterTitle(); gwaveform_pedestalAdjusted[sec]->GetYaxis()->SetTitle("ADC counts"); gwaveform_pedestalAdjusted[sec]->GetYaxis()->CenterTitle(); gwaveform_pedestalAdjusted[sec]->SetMarkerColor(colors[sec%6]); gwaveform_pedestalAdjusted[sec]->SetMarkerStyle(20 + sec/6); gwaveform_pedestalAdjusted[sec]->SetMarkerSize(0.4); gwaveform_pedestalAdjusted[sec]->SetLineColor(colors[sec%6]); gwaveform_pedestalAdjusted[sec]->SetLineStyle(1); gwaveform_pedestalAdjusted[sec]->SetLineWidth(1); gwaveform_normalized[sec] = new TGraph(100); sprintf(gname,"gwaveform_normalized%2.2d",sec+1); gwaveform_normalized[sec]->SetName(gname); gwaveform_normalized[sec]->SetTitle(""); gwaveform_normalized[sec]->GetXaxis()->SetTitle("sample #"); gwaveform_normalized[sec]->GetXaxis()->CenterTitle(); gwaveform_normalized[sec]->GetYaxis()->SetTitle("ADC counts"); gwaveform_normalized[sec]->GetYaxis()->CenterTitle(); gwaveform_normalized[sec]->SetMarkerColor(colors[sec%6]); gwaveform_normalized[sec]->SetMarkerStyle(20 + sec/6); gwaveform_normalized[sec]->SetMarkerSize(0.4); gwaveform_normalized[sec]->SetLineColor(colors[sec%6]); gwaveform_normalized[sec]->SetLineStyle(1); gwaveform_normalized[sec]->SetLineWidth(1); } TCanvas *canvasAll = new TCanvas("canvasAll","canvasAll",1200,600); canvasAll->SetRightMargin(0.02); canvasAll->SetTopMargin(0.02); canvasAll->SetLeftMargin(0.10); canvasAll->SetBottomMargin(0.10); canvasAll->Draw(); TLegend *legend = new TLegend(0.15,0.70,0.95,0.95); legend->SetTextColor(kBlack); legend->SetTextAlign(12); legend->SetTextFont(132); legend->SetTextSize(0.040); legend->SetBorderSize(0); legend->SetFillStyle(0); legend->SetNColumns(6); char text[200]; TLatex *latex = new TLatex(); latex->SetTextColor(kBlack); latex->SetTextAlign(12); latex->SetTextFont(132); latex->SetTextSize(0.040); latex->SetNDC(0); TLine *line_time0[NSECTORS]; for(Int_t sec=0;secSetLineColor(colors[sec%6]); line_time0[sec]->SetLineStyle(4); line_time0[sec]->SetLineWidth(1); } /////////////////////////////////////////////////////////////// // // // Histograms for waveforms // // // /////////////////////////////////////////////////////////////// char hname[400]; TH1F *hnvchannels_for_average = new TH1F("hnvchannels_for_average",";# channels used for average pedestal calculation",NSECTORS,0.5,NSECTORS+0.5); TH2F *hnvchannels_for_average_nwaveforms = new TH2F("hnvchannels_for_average_nwaveforms",";# waveforms; # channels used for average pedestal calculation",NSECTORS,0.5,NSECTORS+0.5,NSECTORS,0.5,NSECTORS+0.5); // Threshold for using sector in average waveform calculation // vs // number of sectors used in average waveform calculation TH2F *hthreshold_nvchannels_for_average = new TH2F("hthreshold_nvchannels_for_average",";# channels used for average pedestal calculation;threshold",NSECTORS+1,-0.5,NSECTORS+0.5,19,5,100); // Event numbers when #waveforms is 0 TH1F *heventnum_nwaveforms_0 = new TH1F("heventnum_nwaveforms_0",";event num",2000,0,20*1000*1000); // Event numbers when #waveforms is not 32 TH1F *heventnum_nwaveforms_not_32 = new TH1F("heventnum_nwaveforms_not_32",";event num",2000,0,20*1000*1000); // Event numbers when #sectors used for average is 24 TH1F *heventnum_nvchannels_for_average_less_than_24 = new TH1F("heventnum_nvchannels_for_average_less_than_24",";event num",2000,0,20*1000*1000); // Event numbers when #sectors used for average is 32 TH1F *heventnum_nvchannels_for_average_is_32 = new TH1F("heventnum_nvchannels_for_average_is_32",";event num",2000,0,20*1000*1000); // Event numbers when #sectors used for average is 0 TH1F *heventnum_nvchannels_for_average_is_0 = new TH1F("heventnum_nvchannels_for_average_is_0",";event num",2000,0,20*1000*1000); // Event numbers when #signals > 5 TH1F *heventnum_nSignals_greater_than_5 = new TH1F("heventnum_nSignals_greater_than_5",";event num",2000,0,20*1000*1000); // Check which sector(s) were above readout threshold TH1F *hsector_with_pedestal_above_threshold = new TH1F("hsector_with_pedestal_above_threshold",";sector with pedestal above threshold;",NSECTORS,0.5,NSECTORS+0.5); TH1F *hsector_with_pedestal_above_threshold_nvchannels_for_average_is_32 = new TH1F("hsector_with_pedestal_above_threshold_nvchannels_for_average_is_32",";sector with pedestal above threshold;",NSECTORS,0.5,NSECTORS+0.5); // Pedestal, maxcounts for raw waveforms TH1F *hpedestal[NSECTORS]; TH1F *hmaxcounts[NSECTORS]; TH1F *hdiffmaxcounts_pedestal[NSECTORS]; TH2F *hmaxcounts_pedestal[NSECTORS]; // Diff of maxcounts and pedestal when no tail TH1F *hdiffmaxcounts_pedestal_notail[NSECTORS]; // Pedestal, maxcounts for normalized waveforms // TH1F *hpedestal_normalized[NSECTORS]; TH1F *hwaveform_normalized[NSECTORS]; TH1F *hdiffmaxcounts_pedestal_normalized[NSECTORS]; TH1F *hdiffmaxcounts_pedestal_normalized_oneHit[NSECTORS]; TH1F *hnormalized_integral[NSECTORS]; TH1F *hhittime[NSECTORS]; TH2F *hsector_hittime = new TH2F("hsector_hittime",";hit time (ns); hit sector",400,0,400,NSECTORS,0.5,NSECTORS+0.5); // When #PS pairs, #PS pairs >= 1 TH1F *hdiffmaxcounts_pedestal_normalized_PSHIT[NSECTORS]; TH1F *hnormalized_integral_PSHIT[NSECTORS]; TH2F *hsector_hittime_PSHIT = new TH2F("hsector_hittime_PSHIT",";hit time (ns); hit sector",400,0,400,NSECTORS,0.5,NSECTORS+0.5); // When within coherent peak of PS TH1F *hdiffmaxcounts_pedestal_normalized_COHPEAK[NSECTORS]; TH1F *hnormalized_integral_COHPEAK[NSECTORS]; TH2F *hsector_hittime_COHPEAK = new TH2F("hsector_hittime_COHPEAK",";hit time (ns); hit sector",400,0,400,NSECTORS,0.5,NSECTORS+0.5); // Comparison of normalized (maxcounts - pedestal) and hit amplitude from fit TH2F *hdiffmaxcounts_pedestal_fitAmp[NSECTORS]; // Comparison of normalized (maxcounts - pedestal) and normalized integral TH2F *hdiffmaxcounts_pedestal_normalized_integral[NSECTORS]; // Comparison of normalized integral and hit amplitude from fit TH2F *hnormalized_integral_fitAmp[NSECTORS]; for(Int_t i=0;i= 1 sprintf(hname,"hdiffmaxcounts_pedestal_normalized_PSHIT%2.2d",i+1); hdiffmaxcounts_pedestal_normalized_PSHIT[i] = new TH1F(hname,";maxcounts - pedestal (normalized) (ADC counts);",8000,-20,SIGNAL_MAX-20); sprintf(hname,"hnormalized_integral_PSHIT%2.2d",i+1); hnormalized_integral_PSHIT[i] = new TH1F(hname,";normalized integral;",5000,-1500,MAXCOUNTS_MAX*65-1500); // When within coherent peak of PS sprintf(hname,"hdiffmaxcounts_pedestal_normalized_COHPEAK%2.2d",i+1); hdiffmaxcounts_pedestal_normalized_COHPEAK[i] = new TH1F(hname,";maxcounts - pedestal (normalized) (ADC counts);",8000,-20,SIGNAL_MAX-20); sprintf(hname,"hnormalized_integral_COHPEAK%2.2d",i+1); hnormalized_integral_COHPEAK[i] = new TH1F(hname,";normalized integral;",5000,-1500,MAXCOUNTS_MAX*65-1500); // Comparison of normalized (maxcounts - pedestal) and hit amplitude from fit sprintf(hname,"hdiffmaxcounts_pedestal_fitAmp%2.2d",i+1); hdiffmaxcounts_pedestal_fitAmp[i] = new TH2F(hname,";fit amplitude;maxcounts - pedestal",250,0,SIGNAL_MAX,250,0,SIGNAL_MAX); // Comparison of normalized (maxcounts - pedestal) and normalized integral sprintf(hname,"hdiffmaxcounts_pedestal_normalized_integral%2.2d",i+1); hdiffmaxcounts_pedestal_normalized_integral[i] = new TH2F(hname,";normalized integral;maxcounts - pedestal;",250,-1500,MAXCOUNTS_MAX*65-1500,250,0,SIGNAL_MAX); // Comparison of normalized integral and hit amplitude from fit sprintf(hname,"hnormalized_integral_fitAmp%2.2d",i+1); hnormalized_integral_fitAmp[i] = new TH2F(hname,";fit amplitude;normalized integral",250,0,MAXCOUNTS_MAX,250,-1500,MAXCOUNTS_MAX*65-1500); } // Multiplicity of signals // Signals are defined as having // [normalized maxcounts] - [normalized pedestal] // larger than MULTIPLICITY_THRESHOLD (set with -m) TH1F *hnwaveforms = new TH1F("hnwaveforms",";#waveforms;",NMAX,-0.5,NMAX-0.5); TH1F *hnSignals = new TH1F("hnSignals",";# signals;",NMAX,-0.5,NMAX-0.5); TH2F *hmultiplicity = new TH2F("hmultiplicity",";sector #1;sector #2",NSECTORS,0.5,NSECTORS+0.5,NSECTORS,0.5,NSECTORS+0.5); TH2F *hnvchannels_for_average_nSignals = new TH2F("hnvchannels_for_average_nSignals",";# signals;# channels used for average pedestal calculation",NSECTORS,0.5,NSECTORS+0.5,NSECTORS,0.5,NSECTORS+0.5); // PSC TH1F *hnPSCPair = new TH1F("hnPSCPair",";# of PSC pairs;",20,-0.5,19.5); TH2F *hPSCLvsR = new TH2F("hPSCLvsR",";PSC R;PSC L",8,0.5,8.5,8,0.5,8.5); TH2F *hPSCtLvstR = new TH2F("hPSCtLvstR",";t_{R} (ns);t_{L} (ns)",200,-10,10,200,-10,10); TH1F *hPSCtdiff = new TH1F("hPSCtdiff",";PSC #Deltat_{L-R} (ns)",200,-10,10); // PS TH1F *hnPSPair = new TH1F("hnPSPair",";# of PS pairs;",20,-0.5,19.5); TH2F *hPSlvsr = new TH2F("hPSlvsr",";PS r;PS l",145,0.5,145.5,145,0.5,145.5); TH2F *hPSElvsEr = new TH2F("hPSElvsEr",";E_{R};E_{L}",200,1.0,3.0,200,1.0,3.0); TH1F *hPSEsum = new TH1F("hPSEsum",";PS E sum (GeV);",400,2.0,6.0); // PSC vs PS TH2F *hnPSCPair_nPSPair = new TH2F("hnPSCPair_nPSPair",";# PS pair;# PSC pair",20,-0.5,19.5,20,-0.5,19.5); TH2F *hPSCL_PSl = new TH2F("hPSCL_PSl",";PS l;PSC L",145,0.5,145.5,8,0.5,8.5); TH2F *hPSCR_PSr = new TH2F("hPSCR_PSr",";PS r;PSC R",145,0.5,145.5,8,0.5,8.5); TH2F *hPSCtdiff_PSEsum = new TH2F("hPSCtdiff_PSEsum",";PS E sum (GeV);PSC #Deltat_{L-R} (ns)",400,2.0,6.0,200,-10,10); // Histograms for Mike's analysis TH1F *hTemp1 = new TH1F("hTemp1","",100,0.0,400.0); TH1F *hVarX = new TH1F("hVarX",";variance",1000,0.0,100.0); TH2F *hSecVsVarX = new TH2F("hSecVsVarX",";variance;sector",1000,0.0,100.0,32,0.5,32.5); TH1F *hfitAmpAll = new TH1F("hfitAmpAll",";fitAmp",1000,0,MAXCOUNTS_MAX); TH2F *hSectorVsfitAmp = new TH2F("hSectorVsfitAmp",";fitAmp;sector",200,0,MAXCOUNTS_MAX,NSECTORS,0.5,NSECTORS+0.5); TH1D *hPrimeTime = new TH1D("hPrimeTime",";t_{0} (ns)",1600,0.0,400.0); TH2D *hSectorVsPrimeTime = new TH2D("hSectorVsPrimeTime",";t_{0} (ns);sector",1600,0.0,400.0,NSECTORS,0.5,NSECTORS+0.5); TH1D *hGammaRise = new TH1D("hGammaRise",";#Gamma_{r} (ns^{-1})",200,0,0.04); TH2D *hSectorVsGammaRise = new TH2D("hSectorVsGammaRise",";#Gamma_{r} (ns^{-1});sector",200,0,0.04,NSECTORS,0.5,NSECTORS+0.5); TH1D *hGammaFall = new TH1D("hGammaFall",";#Gamma_{f} (ns^{-1})",200,0,0.0002); // fixed TH1D *hcVal = new TH1D("hcVal",";constant",200,-5,5); TH2D *hSectorVscVal = new TH2D("hSectorVscVal",";constant;sector",200,-5,5,NSECTORS,0.5,NSECTORS+0.5); TH1F *hfitAmp[NSECTORS]; for(Int_t i=0;iGetEntries() && n < MAXEVENTS;n++){ intree->GetEntry(n); if(verbose) cout << "entry = " << n << endl; // This is for Mike's analysis tree tpol1.nSectorHits = 0; // Clear vector of sectors and fit params vsector_fitparams.clear(); hnwaveforms->Fill(nwaveforms); if(nwaveforms==0){ cout << "no hits for event # " << setw(8) << eventnum << endl; continue; } if(n % 10000 == 0) cout << "------------------- entry " << n << " / " << (intree->GetEntries() < MAXEVENTS ? intree->GetEntries() : MAXEVENTS) << endl; if(nwaveforms == 0){ heventnum_nwaveforms_0->Fill(eventnum); } if(nwaveforms != NSECTORS){ OUT_STRANGE << "event " << setw(8) << eventnum << " : nwaveforms was not " << NSECTORS << " but " << nwaveforms << endl; heventnum_nwaveforms_not_32->Fill(eventnum); if(verbose){ cout << "nwaveforms was not " << NSECTORS << " but " << nwaveforms << endl; } } //////////////////////////////////////////////// // // // Fill PS, PSC histograms // // // //////////////////////////////////////////////// // PSC hnPSCPair->Fill(nhits_PSC); for(Int_t nPSC=0;nPSCFill(iRight_PSC[nPSC],iLeft_PSC[nPSC]); hPSCtLvstR->Fill(t_iRight_PSC[nPSC],t_iLeft_PSC[nPSC]); hPSCtdiff->Fill(t_iRight_PSC[nPSC] - t_iLeft_PSC[nPSC]); } // PS hnPSPair->Fill(nhits_PS); for(Int_t nPS=0;nPSFill(iRight_PS[nPS],iLeft_PS[nPS]); hPSElvsEr->Fill(e_iRight_PS[nPS],e_iLeft_PS[nPS]); hPSEsum->Fill(e_iRight_PS[nPS] + e_iLeft_PS[nPS]); } hnPSCPair_nPSPair->Fill(nhits_PS, nhits_PSC); // PSC vs PS for(Int_t nPSC=0;nPSCFill(iLeft_PS[nPS],iLeft_PSC[nPSC]); hPSCR_PSr->Fill(iRight_PS[nPS],iRight_PSC[nPSC]); hPSCtdiff_PSEsum->Fill(e_iRight_PS[nPS] + e_iLeft_PS[nPS], t_iRight_PSC[nPSC] - t_iLeft_PSC[nPSC]); } } // Set up flags for PS hits bool IS_PSHIT = false; bool IS_COHPEAK = false; // Was there both a PS and PSC pair? if(nhits_PSC >=1 && nhits_PS >=1){ IS_PSHIT = true; // Was there a PS pair where the sum of PS energies was in the coherent peak? for(Int_t iPS=0;iPSSetPoint(sample,0,0); gwaveform_pedestalAdjusted[sec]->SetPoint(sample,0,0); gwaveform_normalized[sec]->SetPoint(sample,0,0); } gwaveform_average->SetPoint(sample,0,0); } // Waveforms are drawn and saved as pdf files only if // the difference of [maximum fADC counts] - [pedestal] // is larger than DRAW_THRESHOLD (specify with -t) Bool_t draw_waveforms = false; // Vector holding channels that will be used to // calculate average waveform std::vector vchannels_for_average; // Initialize array that holds # sectors used for average waveform for a given threshold for(Int_t nthresholds=0;nthresholds 100) hadSignal = true; // If there was a hit and the final waveform is close to the // pedestal, we had a signal with no tail if(sample == 99){ if(hadSignal && fabs(waveform[i][sample] - pedestal[i]) < 10){ notailSignal = true; } } } // end of loop over samples of waveform if(notailSignal==true){ hdiffmaxcounts_pedestal_notail[sec-1]->Fill(maxcounts[i] - pedestal[i],sec); // draw_waveforms = true; } if(fabs(maxcounts[i] - pedestal[i]) > DRAW_THRESHOLD){ draw_waveforms = true; if(verbose) cout << "writing out event " << eventnum << endl; } } // end of loop over hits within entry // If # of sectors used for average is low, draw waveforms. if(vchannels_for_average.size() < 24){ OUT_STRANGE << "event " << setw(8) << eventnum << " : #sectors used for average = " << vchannels_for_average.size() << endl; cout << "event " << eventnum << " had #sectors used for average = " << vchannels_for_average.size() << endl; draw_waveforms = true; heventnum_nvchannels_for_average_less_than_24->Fill(eventnum); } for(Int_t i=0;i READOUT_THRESHOLD) hsector_with_pedestal_above_threshold->Fill(sector[i]); } if(vchannels_for_average.size() == NSECTORS){ // See which events had 32 sectors used for average (no signal, fluctuation of baseline) OUT_STRANGE << "event " << setw(8) << eventnum << " : #sectors used for average = " << vchannels_for_average.size() << endl; // cout << "event " << eventnum << " had #sectors used for average = " << vchannels_for_average.size() << endl; // Too many // draw_waveforms = true; heventnum_nvchannels_for_average_is_32->Fill(eventnum); // Make histogram of sectors that went above readout threshold for(Int_t i=0;i READOUT_THRESHOLD) hsector_with_pedestal_above_threshold_nvchannels_for_average_is_32->Fill(sector[i]); } } hnvchannels_for_average->Fill(vchannels_for_average.size()); hnvchannels_for_average_nwaveforms->Fill(nwaveforms,vchannels_for_average.size()); for(Int_t nthresholds=0;nthresholdsGetBinContent(xbin,ybin); hthreshold_nvchannels_for_average->SetBinContent(xbin,ybin,currentbincontent + 1); } if(verbose) cout << "# of channels used for event " << eventnum << " = " << vchannels_for_average.size() << endl; if(vchannels_for_average.size() == 0){ OUT_STRANGE << "event " << setw(8) << eventnum << " : no channels used to calculate average" << endl; cout << "No channels used to calculate average for event " << eventnum << "...." << endl; heventnum_nvchannels_for_average_is_0->Fill(eventnum); } Double_t average_waveform[100]; for(Int_t i=0;i<100;i++){ average_waveform[i] = 0; } ///////////////////////////////////////////////////////////// // // // Adjust pedestals to 100 for all channels // // // ///////////////////////////////////////////////////////////// // Adjust all channels to have pedestal of PEDESTAL (= 100) for(Int_t i=0;iGetPoint(sample,x,y); gwaveform_pedestalAdjusted[sec-1]->SetPoint(sample,x,y - pedestal[i] + PEDESTAL); } } ///////////////////////////////////////////////////////////// // // // Calculate average waveform // // Use pedestal-adjusted waveforms // // // ///////////////////////////////////////////////////////////// for(Int_t sample=0;sample<100;sample++){ // Get the channels to be included in calculation of average // vchannels_for_average contains the sector numbers to be used for(std::vector::iterator it = vchannels_for_average.begin() ; it != vchannels_for_average.end(); ++it){ Int_t sec = *it; // SINCE WAVEFORM GRAPHS HAVE INDEX OF [SECTOR # - 1], // WE CAN USE THE SECTOR NUMBERES SAVED IN vchannels_for_average // TO MATCH UP THE CORRECT WAVEFORMS gwaveform_pedestalAdjusted[sec-1]->GetPoint(sample,x,y); average_waveform[sample] += y; } // end of loop over nwaveforms // cout << "vchannels_for_average.size() = " << vchannels_for_average.size() << endl; if(vchannels_for_average.size() == 0){ average_waveform[sample] = 0; }else{ average_waveform[sample] /= vchannels_for_average.size(); } gwaveform_average->SetPoint(sample,4. * (sample + 1),average_waveform[sample]); if(average_waveform[sample] != average_waveform[sample]) cout << sample << " " << average_waveform[sample] << endl; if(verbose) cout << "average " << 4. * (sample + 1) << " " << average_waveform[sample] << endl; } // end of loop over samples to calculate average // Set all fsigFun_copy to 0. // These are set only if there was a fit to the signal. for(Int_t sec=0;secGetPoint(sample,x,y); gwaveform_average->GetPoint(sample,average_x,average_y); gwaveform_normalized[sec-1]->SetPoint(sample,4. * (sample + 1), y - average_y); // Fill histogram for normalized waveform hwaveform_normalized[sec-1]->Fill(y - average_y); // Get normalized maxcounts if(y - average_y > maxcounts_normalized[sec-1]) maxcounts_normalized[sec-1] = y - average_y; // Calculate normalized integrals normalized_integral += y - average_y; // Calculate pedestal if(sample<5){ pedestal_normalized[sec-1] += y - average_y; } // Find if there was a hit if(hadHit == false && y - average_y > FIT_THRESHOLD){ // This ensures that we only fill this for // the first sample above FIT_THRESHOLD hadHit = true; hhittime[sec-1]->Fill(4. * (sample + 1)); hsector_hittime->Fill(4. * (sample + 1),sec); ////////////////////////////////////////////////// // When #PS pairs, #PS pairs >= 1 // ////////////////////////////////////////////////// if(IS_PSHIT) hsector_hittime_PSHIT->Fill(4. * (sample + 1),sec); ////////////////////////////////////////////////// // When PS E is in coh. peak // ////////////////////////////////////////////////// if(IS_COHPEAK) hsector_hittime_COHPEAK->Fill(4. * (sample + 1),sec); } // end of had hit } // end of loop over samples // Normalized pedestal pedestal_normalized[sec-1] /= 5.; // Fill difference of normalized pedestal and normalized maxcounts hdiffmaxcounts_pedestal_normalized[sec-1]->Fill(maxcounts_normalized[sec-1] - pedestal_normalized[sec-1]); // Fill difference of normalized pedestal and normalized maxcounts // when there was only one hit. // This is determined by how many channels were used in // the pedestal average // (vchannels_for_average.size() should be NSECTORS-1) if(vchannels_for_average.size() == NSECTORS-1) hdiffmaxcounts_pedestal_normalized_oneHit[sec-1]->Fill(maxcounts_normalized[sec-1] - pedestal_normalized[sec-1]); // Subtract pedestal to get normalized integrals normalized_integral -= pedestal_normalized[sec-1] * 100.; hnormalized_integral[sec-1]->Fill(normalized_integral); if(verbose){ if(normalized_integral > 40000.) cout << normalized_integral << endl; } ////////////////////////////////////////////////// // When #PS pairs, #PS pairs >= 1 // ////////////////////////////////////////////////// if(IS_PSHIT){ hdiffmaxcounts_pedestal_normalized_PSHIT[sec-1]->Fill(maxcounts_normalized[sec-1] - pedestal_normalized[sec-1]); hnormalized_integral_PSHIT[sec-1]->Fill(normalized_integral); } ////////////////////////////////////////////////// // When PS E is in coh. peak // ////////////////////////////////////////////////// if(IS_COHPEAK){ hdiffmaxcounts_pedestal_normalized_COHPEAK[sec-1]->Fill(maxcounts_normalized[sec-1] - pedestal_normalized[sec-1]); hnormalized_integral_COHPEAK[sec-1]->Fill(normalized_integral); } /////////////////////////////////////////////////////////// // // // Start of Mike's fit analysis. // // // /////////////////////////////////////////////////////////// // The normalized waveforms are contained in the // TGraphs gwaveform_normalized[sec-1] // The pedestals are extremely close to 0. for (Int_t tNum = 1; tNum<=100; tNum++){ gwaveform_normalized[sec-1]->GetPoint(tNum-1,x,y); hTemp1->SetBinContent(tNum,y); } double diffTest = hTemp1->GetMaximum() - hTemp1->GetMinimum(); //Get Estimate of variance for no signal if(diffTest < FIT_THRESHOLD){ double fadcSum = 0.0; for (Int_t tNum = 1; tNum<=100; tNum++){ gwaveform_normalized[sec-1]->GetPoint(tNum-1,x,y); fadcSum += y; } double fadcAve = fadcSum/100.0; //Let xVal = fadc - fadcAve double xSumSq = 0; for (Int_t tNum = 1; tNum<=100; tNum++){ gwaveform_normalized[sec-1]->GetPoint(tNum-1,x,y); xSumSq += pow(y - fadcAve,2.0); } double varX = xSumSq/100.0; hVarX->Fill(varX); hSecVsVarX->Fill(varX,sec*1.0); } if(diffTest > FIT_THRESHOLD){ nCount++; nCountFits++; if(verbose) cout << "nCountFits = " << nCountFits << " , diffTest = " << diffTest << endl; if (nCount > NMAXKEEP) nCount = 1; if (nCount <= NMAXKEEP){ tpol1.nSectorHits++; tpol1.sector[tpol1.nSectorHits-1] = sec; tpol1.phi[tpol1.nSectorHits-1] = sectorToPhi(sec); double difVal = 0.0; for (Int_t tNum = 1; tNum<=100; tNum++){ if (tNum > 1) difVal = hTemp1->GetBinContent(tNum) - hTemp1->GetBinContent(tNum-1); hT[nCount-1]->SetBinContent(tNum,difVal*tpolEnergyCal(sec)); gwaveform_normalized[sec-1]->GetPoint(tNum-1,x,y); hTR[nCount-1]->SetBinContent(tNum,y); //hTR[nCount-1]->SetBinError(tNum,sqrt(8.0)); } //-->START FIT double sigIntLow40 = hTR[nCount-1]->Integral(1,10)/10.0; int binTimeSeed = hT[nCount-1]->GetMaximumBin(); double startTimeSeed = hT[nCount-1]->GetBinCenter(binTimeSeed) - 23.0; //Make parameter limits sigFun->SetParLimits(0,0.0,10000); //fADC counts // sigFun->SetParLimits(1,1.0e-2,4.0e-2); // sigFun->SetParLimits(3,sigIntLow40*0.5,sigIntLow40*1.5); // sigFun->SetParLimits(4,-50.0,325.0); if(verbose){ cout << "maxcounts_normalized[sec-1] = " << maxcounts_normalized[sec-1] << endl; cout << "sigIntLow40 = " << sigIntLow40 << endl; cout << "startTimeSeed = " << startTimeSeed << endl; } //Seed parameters sigFun->SetParameter(0,maxcounts_normalized[sec-1] + 50); //fADC counts sigFun->SetParameter(1,0.028); //GammaRise sigFun->FixParameter(2,0.000167); //GammaFall (1/6us) sigFun->SetParameter(3,sigIntLow40); //offset sigFun->SetParameter(4,startTimeSeed); //startTime //Fit the function hTR[nCount-1]->Fit("sigFun","BREQ","",0.0,400.0); //Load up the fit values into the tree tpol1.energy[tpol1.nSectorHits-1] = sigFun->GetParameter(0)*tpolEnergyCal(sec); tpol1.energyErr[tpol1.nSectorHits-1] = sigFun->GetParError(0)*tpolEnergyCal(sec); tpol1.riseTime[tpol1.nSectorHits-1] = sigFun->GetParameter(1); tpol1.riseTime[tpol1.nSectorHits-1] = sigFun->GetParError(1); tpol1.energyOffset[tpol1.nSectorHits-1] = sigFun->GetParameter(3)*tpolEnergyCal(sec); tpol1.energyOffsetErr[tpol1.nSectorHits-1] = sigFun->GetParError(3)*tpolEnergyCal(sec); tpol1.startTime[tpol1.nSectorHits-1] = sigFun->GetParameter(4); tpol1.startTimeErr[tpol1.nSectorHits-1] = sigFun->GetParError(4); if(verbose){ // Check difference of initial seeds and fitted values cout << "fit time initial seed: " << startTimeSeed << endl << "fitted time : " << sigFun->GetParameter(4) << endl << "diff : " << startTimeSeed - sigFun->GetParameter(4) << endl; cout << "fit amp initial seed : " << maxcounts_normalized[sec-1] + 50 << endl << "fit amp : " << sigFun->GetParameter(0) << endl << "diff : " << maxcounts_normalized[sec-1] + 50 - sigFun->GetParameter(0) << endl; } if (tpol1.energy[tpol1.nSectorHits-1] > 150.0){ if(verbose) cout << sigFun->GetParameter(0) << "\t" << sigFun->GetParameter(1) << "\t" << sigFun->GetParameter(2) << "\t" << sigFun->GetParameter(3) << "\t" << sigFun->GetParameter(4) << endl; hPrimeTime->Fill(tpol1.startTime[tpol1.nSectorHits-1]); hSectorVsPrimeTime->Fill(tpol1.startTime[tpol1.nSectorHits-1],sec); hfitAmpAll->Fill(sigFun->GetParameter(0)); hSectorVsfitAmp->Fill(sigFun->GetParameter(0),sec); hGammaRise->Fill(sigFun->GetParameter(1)); hSectorVsGammaRise->Fill(sigFun->GetParameter(1),sec); hGammaFall->Fill(sigFun->GetParameter(2)); hcVal->Fill(sigFun->GetParameter(3)); hSectorVscVal->Fill(sigFun->GetParameter(3),sec); } if(verbose) cout<<"startTime = "<GetChisquare(); int ndfVal = sigFun->GetNDF(); tpol1.confidenceLevel[tpol1.nSectorHits-1] = sigFun->GetProb(); tpol1.chi2Reduced[tpol1.nSectorHits-1] = chi2Val/ndfVal; if(verbose){ cout << " fit to sector " << sec << " chi2 = " << chi2Val << " ndf = " << ndfVal << endl; } // Fill in fitted sector hfitFrequency->Fill(sec); // Fill in fitAmp hfitAmp[sec-1]->Fill(sigFun->GetParameter(0)); // Fill in maxcounts - pedestal vs fit amplitude if(verbose) cout << tpol1.energy[tpol1.nSectorHits-1] << " " << maxcounts_normalized[sec-1] - pedestal_normalized[sec-1] << endl; hdiffmaxcounts_pedestal_fitAmp[sec-1]->Fill(sigFun->GetParameter(0), // tpol1.energy[tpol1.nSectorHits-1], maxcounts_normalized[sec-1] - pedestal_normalized[sec-1]); // Fill in maxcounts - pedestal vs normalized integral hdiffmaxcounts_pedestal_normalized_integral[sec-1]->Fill(normalized_integral, maxcounts_normalized[sec-1] - pedestal_normalized[sec-1]); // Fill in normalized integral vs fit amplitude hnormalized_integral_fitAmp[sec-1]->Fill(sigFun->GetParameter(0),normalized_integral); if(verbose){ cout << "sec = " << sec << " energy = " << tpol1.energy[tpol1.nSectorHits-1] << " maxcounts - ped = " << maxcounts_normalized[sec-1] - pedestal_normalized[sec-1] << endl; } // !!! THIS DOES NOT WORK !!! // FOR SOME REASON WE CANNOT COPY sigFun DIRECTLY // AND DRAW IT. SAVE FIT PARAMETERS INSTEAD. // Copy sigFun to fsigFun_copy[sec-1] so that we can draw it. // sprintf(hname,"fsigFun_copy%8.8d_%2.2d",eventnum,sec); // fsigFun_copy[sec-1] = (TF1*)sigFun->Clone(hname); vector params; for(int p=0;p<5;p++) params.push_back(sigFun->GetParameter(p)); // Create pair of sector number and params pair > thispair = std::make_pair(sec,params); // Add into vector of sector numbers and params vsector_fitparams.push_back(thispair); } // end of nCounts <= NMAXKEEP } // end of diffTest < FIT_THRESHOLD } // end of loop over nwaveforms in an event if (tpol1.nSectorHits > 0) { tpol1.nhits_PS = nhits_PS; tpol1.nhits_PSC = nhits_PSC; for (int pscIndex = 0; pscIndex < tpol1.nhits_PSC; pscIndex++){ tpol1.iLeft_PSC[pscIndex] = iLeft_PSC[pscIndex]; tpol1.iRight_PSC[pscIndex] = iRight_PSC[pscIndex]; tpol1.t_iLeft_PSC[pscIndex] = t_iLeft_PSC[pscIndex]; tpol1.t_iRight_PSC[pscIndex] = t_iRight_PSC[pscIndex]; } for (int psIndex = 0; psIndex < tpol1.nhits_PS; psIndex++){ tpol1.iLeft_PS[psIndex] = iLeft_PS[psIndex]; tpol1.iRight_PS[psIndex] = iRight_PS[psIndex]; tpol1.e_iLeft_PS[psIndex] = e_iLeft_PS[psIndex]; tpol1.e_iRight_PS[psIndex] = e_iRight_PS[psIndex]; } outTree->Fill(); } // Still looping over events in tree ////////////////////////////////////////////////////////////// // // // Check how many sectors had maxcount - pedestal // // above some threshold value (multiplicity) // // // ////////////////////////////////////////////////////////////// bool hadHit[NSECTORS]; for(int i=0;i MULTIPLICITY_THRESHOLD){ hadHit[sec-1] = true; nSignals++; } } // end of loop over hits // Plot multiplicity hnSignals->Fill(nSignals); hnvchannels_for_average_nSignals->Fill(nSignals,vchannels_for_average.size()); // If we have more than 5 signals, draw the waveforms. // For 1M 241Am events, there are ~1600 such events if(nSignals > 5){ OUT_STRANGE << "event " << setw(8) << eventnum << " : #signals = " << nSignals << endl; cout << "event " << eventnum << " had #signals = " << nSignals << endl; draw_waveforms = true; heventnum_nSignals_greater_than_5->Fill(eventnum); } // Plot sectors that shared hits if multiplicity = 2 if(nSignals==2){ Int_t sec1 = -999; Int_t sec2 = -999; for(Int_t i=0;iFill(sec1,sec2); if(sec1 == sec2){ cout << "sec1 = " << sec1 << " sec2 = " << sec2 << endl; for(int num=0;numSetHeader(text); for(Int_t i=0;iSetMinimum(0); gwaveform[sec-1]->SetMaximum(MAXCOUNTS_MAX); gwaveform[sec-1]->GetXaxis()->SetLimits(0,400); if(i==0) gwaveform[sec-1]->Draw("ALP"); else gwaveform[sec-1]->Draw("LP"); sprintf(text,"sector %2.2d",sector[sec-1]); legend->AddEntry(gwaveform[sec-1],text,"LP"); } gwaveform_average->Draw("LP"); legend->AddEntry(gwaveform_average,"average","LP"); legend->Draw("same"); sprintf(text,"mkdir -p figures/%s/waveforms",basename.c_str()); system(text); sprintf(text,"figures/%s/waveforms/%8.8d___01___rawwaveform.pdf",basename.c_str(),eventnum); canvasAll->SaveAs(text); //_________________________________________________________________________________ // Draw all pedestal-adjusted waveforms top of each other canvasAll->cd(); legend->Clear(); sprintf(text,"pedestal subtracted waveforms event%8.8d",eventnum); legend->SetHeader(text); for(Int_t i=0;iGetXaxis()->SetLimits(0,400); gwaveform_pedestalAdjusted[sec-1]->SetMinimum(0); gwaveform_pedestalAdjusted[sec-1]->SetMaximum(MAXCOUNTS_MAX); if(i==0) gwaveform_pedestalAdjusted[sec-1]->Draw("ALP"); else gwaveform_pedestalAdjusted[sec-1]->Draw("LP"); sprintf(text,"sector %2.2d",sec); legend->AddEntry(gwaveform_pedestalAdjusted[sec-1],text,"LP"); } gwaveform_average->Draw("LP"); legend->AddEntry(gwaveform_average,"average","LP"); legend->Draw("same"); legend->Draw("same"); sprintf(text,"mkdir -p figures/%s/waveforms",basename.c_str()); system(text); sprintf(text,"figures/%s/waveforms/%8.8d___02___pedestalAdjusted.pdf",basename.c_str(),eventnum); canvasAll->SaveAs(text); //_________________________________________________________________________________ // Draw all waveforms with average subtracted top of each other canvasAll->cd(); legend->Clear(); sprintf(text,"normalized waveforms event%8.8d",eventnum); legend->SetHeader(text); bool timeSector27Less66 = false; bool timeSector27Greater66 = false; for(Int_t i=0;iGetXaxis()->SetLimits(0,400); gwaveform_normalized[sec-1]->SetMinimum(-20); gwaveform_normalized[sec-1]->SetMaximum(SIGNAL_MAX); if(i==0) gwaveform_normalized[sec-1]->Draw("ALP"); else gwaveform_normalized[sec-1]->Draw("LP"); sprintf(text,"sector %2.2d",sec); legend->AddEntry(gwaveform_normalized[sec-1],text,"LP"); // If fit was performed, draw fit function. // vsector_fitparams contains vector of sectors and fitparams vector thisparams; bool foundFit = false; for (std::vector > >::iterator it = vsector_fitparams.begin() ; it != vsector_fitparams.end(); ++it){ std::pair > thispair = *it; if(thispair.first == sec){ foundFit = true; for(int p=0;pSetParameter(p,thisparams[p]); } fsigFun_copy[sec-1]->SetLineColor(colors[(sec-1)%6]); fsigFun_copy[sec-1]->SetLineStyle(2); fsigFun_copy[sec-1]->SetLineWidth(1); fsigFun_copy[sec-1]->Draw("same"); // fsigFun_copy->Write(); // Draw time sprintf(text,"t_{0} = %5.3f",fsigFun_copy[sec-1]->GetParameter(4)); latex->SetTextColor(colors[(sec-1)%6]); latex->DrawLatex(fsigFun_copy[sec-1]->GetParameter(4) + 1.,gwaveform_normalized[sec-1]->GetMaximum() * sec/NSECTORS,text); if(sec == 27 && fsigFun_copy[sec-1]->GetParameter(4) < 66.) timeSector27Less66 = true; if(sec == 27 && fsigFun_copy[sec-1]->GetParameter(4) > 66.) timeSector27Greater66 = true; // Draw line at t0 line_time0[sec-1]->SetX1(fsigFun_copy[sec-1]->GetParameter(4)); line_time0[sec-1]->SetX2(fsigFun_copy[sec-1]->GetParameter(4)); line_time0[sec-1]->SetY2(gwaveform_normalized[sec-1]->GetYaxis()->GetXmax()); line_time0[sec-1]->Draw("same"); } // end of foundFit } // end of loop over waveforms gwaveform_average->Draw("LP"); legend->AddEntry(gwaveform_average,"average","LP"); legend->Draw("same"); sprintf(text,"mkdir -p figures/%s/waveforms",basename.c_str()); system(text); sprintf(text,"figures/%s/waveforms/%8.8d___03___normalized.pdf",basename.c_str(),eventnum); canvasAll->SaveAs(text); // Save a copy depending on fitted time if(timeSector27Less66){ sprintf(text,"mkdir -p figures/%s/waveforms/smalltime",basename.c_str()); system(text); sprintf(text,"figures/%s/waveforms/smalltime/%8.8d___03___normalized.pdf",basename.c_str(),eventnum); canvasAll->SaveAs(text); } if(timeSector27Greater66){ sprintf(text,"mkdir -p figures/%s/waveforms/largetime",basename.c_str()); system(text); sprintf(text,"figures/%s/waveforms/largetime/%8.8d___03___normalized.pdf",basename.c_str(),eventnum); canvasAll->SaveAs(text); } //_________________________________________________________________________________ } // end of write is true } // end of loop over tree entries cout << "total of " << nCountFits << " waveforms had hits larger than FIT_THRESHOLD = " << FIT_THRESHOLD << endl; ////////////////////////////////////////////////////// // Save // ////////////////////////////////////////////////////// outfile->cd(); outfile->Write(); for(Int_t i=0;iWrite(); gwaveform_pedestalAdjusted[i]->Write(); gwaveform_normalized[i]->Write(); } gwaveform_average->Write(); } void helpMessage(){ cout << "plot_waveform: " << endl; cout << "Options:" << endl; cout << "\t-i infile" << endl; cout << "\t-o outfile" << endl; cout << "\t-v verbose mode" << endl; cout << "\t-I invert polarity for 2015 Spring runs" << endl; cout << "\t-M max # of events to process" << endl; cout << "\t-N maximum number of waveforms to save" << endl; cout << "\t-t threshold to draw and save difference" << endl; cout << "\t-d minimum difference to attempt fit" << endl; cout << "\t-m threshold for multiplicity" << endl; cout << "\t-p maximum difference of maxcounts and pedestal to be used for average pedestal" << endl; cout << "\t-T readout threshold of fADCs" << endl; }