mbd/WAVEFORMS/cal_tail.C

0001 //
0002 // To calibrate the tails from events in previous crossings,
0003 //
0004 //    1. first create the mdbsigXX.txt files (uncomment pileupfile in MbdSig.cc)
0005 //
0006 //    2. get_tails() : extracts the mean of the tails from all time and charge ch's from all events
0007 //                     can be used to play with functional fits
0008 //                     (once fit fcn determined, not necessary and can skip to 2)
0009 //
0010 //    3. calc_residuals() : fit the waveforms from events where there was no event in crossing,
0011 //                          but there was an event in prev. crossing
0012 //                          writes out fit parameters to tfit and qfit tree
0013 //                          makes plots of residuals to fits
0014 //                          only the charge fit makes sense here, since the time can be done by simple
0015 //                          correlation
0016 //
0017 //    4. fit_time_pileupcorr() : find slope of s8 vs s0 in time channels, write out mbd_pileup.calib file
0018 //                               this is already done at end of calc_residuals
0019 //
0020 #include <TGraphErrors.h>
0021 #include <TF1.h>
0022
0023 #include "/phenix/u/chiu/macros/find_th2ridge.C"
0024
0025
0026 const int NFEECH = 256;
0027 //const int NFEECH = 10;
0028 const int NSAMPLES = 16;
0029
0030 TH2 *h2_tail[NFEECH];
0031 TH2 *h2_residuals[NFEECH];
0032 TH2 *h2_s8s0[NFEECH]{nullptr};
0033 TGraphErrors *g_tail[NFEECH];
0034 TF1 *fit[NFEECH] = {nullptr};
0035 TF1 *fit_s8s0[NFEECH] = {nullptr};
0036
0037 // Tree variables
0038 TTree *_tfit_tree{nullptr};
0039 Short_t _tch{-1};
0040 Float_t _tchi2ndf{-1.};
0041 Float_t _tfitpar[10]{0.};
0042 TTree *_qfit_tree{nullptr};
0043 Short_t _qch{-1};
0044 Float_t _qchi2ndf{-1.};
0045 Float_t _qfitpar[10]{0.};
0046
0047 Double_t powerlaw(Double_t *x, Double_t *par)
0048 {
0049   Float_t xx =x[0];
0050
0051   Double_t f = par[0]*pow((par[1]/(par[1]+xx)),par[2]);
0052
0053   return f;
0054 }
0055
0056 //
0057 // Do test fits of the g_tail, which is the mean tail shape from all events
0058 // Used to develop the fit function
0059 //
0060 void fit_test(const int feech = 8, const std::string &fname = "mbdsig_tails.root")
0061 {
0062   TString name;
0063
0064   static TFile *tfile{nullptr};
0065
0066   if ( tfile == nullptr )
0067   {
0068     cout << "Opening " << fname << endl;
0069     tfile = new TFile(fname.c_str(),"READ");
0070   }
0071
0072   name = "g_tail"; name += feech;
0073   TGraphErrors *g = (TGraphErrors*)tfile->Get( name );
0074
0075   TF1 *fit2mean{nullptr};
0076
0077   if ( ((feech/8)%2)==0 ) // time channel
0078   {
0079     //fit2mean = new TF1("fit2mean","gaus+[3]*x+[4]",0,16);
0080     fit2mean = new TF1("fit2mean","gaus+pol2(3)",0,16);
0081     fit2mean->SetParameters(131.5,-18,5.87,-0.002,-0.002/16.,1e-3);
0082     //fit2mean->SetParameters(49.18,-14.544,5.2014,0.0139,-0.00147,4.064e-5);
0083
0084     g->Draw("ap");
0085     g->Fit(fit2mean,"R");
0086
0087     TF1 *gaussian = new TF1("gaussian","gaus",0,16);
0088     gaussian->SetParameters(fit2mean->GetParameter(0),fit2mean->GetParameter(1),fit2mean->GetParameter(2));
0089     gaussian->SetLineColor(4);
0090     gaussian->Draw("same");
0091   }
0092   else
0093   {
0094     //fit2mean = new TF1("fit2mean","gaus+pol2(3)",0,16);
0095     //fit2mean->SetParameters(131.5,-18,5.87,-0.002,-0.002/16.,1e-3);
0096
0097     fit2mean = new TF1("fit2mean","expo + pol2(2)",0,16);
0098     fit2mean->SetParameters(1,-1/16.,1e-2,1e-3,1e-4);
0099
0100     g->Draw("ap");
0101     g->Fit(fit2mean,"R");
0102   }
0103
0104   gPad->SetLogy(1);
0105 }
0106
0107 void fit_tail(TGraph *g, const int feech)
0108 {
0109   int verbose = 0;
0110
0111   TString name;
0112   if ( fit[feech] == nullptr )
0113   {
0114     if ( ((feech/8)%2)==0 ) // time channel
0115     {
0116       name = "fit"; name += feech;
0117       fit[feech] = new TF1(name,"gaus+pol2(3)",-0.1,16);
0118       fit[feech]->SetLineColor(2);
0119       fit[feech]->SetParameters(5.5*g->GetPointY(0),-4.8,2.6,-0.002,-0.002/16.,1e-3);
0120       //fit[feech]->SetParameters(49.18,-14.544,5.2014,0.0139,-0.00147,4.064e-5);
0121     }
0122     else
0123     {
0124       name = "fit"; name += feech;
0125       //fit[feech] = new TF1(name,"expo+pol2(2)",0,16);
0126       //fit[feech]->SetParameters(1,-1/16.,1e-2,1e-3,1e-4);
0127
0128       //fit[feech] = new TF1(name,"expo",0,4);
0129       //fit[feech]->SetParameters(log(g->GetPointY(0)),-1/16);
0130
0131       fit[feech] = new TF1(name,"gaus",-0.1,4.1);
0132       fit[feech]->SetParameters(1.43*g->GetPointY(0),-3.2,3.8);
0133
0134       fit[feech]->SetLineColor(2);
0135     }
0136   }
0137
0138   if ( ((feech/8)%2)==0 ) // time channel
0139   {
0140     if ( verbose )
0141     {
0142       g->Draw("ap");
0143       fit[feech]->SetRange(-0.1,16);
0144       g->Fit(fit[feech],"R");
0145       gPad->SetLogy(1);
0146       Double_t chi2 = fit[feech]->GetChisquare();
0147       Double_t ndf = fit[feech]->GetNDF();
0148       cout << "chi2 ndf " << chi2 << "\t" << ndf << endl;
0149     }
0150     else
0151     {
0152       fit[feech]->SetRange(-0.1,16);
0153       g->Fit(fit[feech],"RNQ");
0154
0155       // save to tree
0156       _tch = feech;
0157       for (int ipar=0; ipar<6; ipar++)
0158       {
0159         _tfitpar[ipar] = static_cast<Float_t>( fit[feech]->GetParameter(ipar) );
0160       }
0161
0162       // copy samp 0 and samp 8 values. samp 8 should be timing maxsamp
0163       _tfitpar[7] = static_cast<Float_t>( g->GetPointY(0) );
0164       _tfitpar[8] = static_cast<Float_t>( g->GetPointY(1) );
0165       _tfitpar[9] = static_cast<Float_t>( g->GetPointY(8) );
0166
0167       h2_s8s0[feech]->Fill( g->GetPointY(0), g->GetPointY(8) );
0168
0169       Double_t chi2 = fit[feech]->GetChisquare();
0170       Double_t ndf = fit[feech]->GetNDF();
0171       _tchi2ndf = chi2/ndf;
0172       _tfit_tree->Fill();
0173     }
0174
0175     if ( verbose )
0176     {
0177       /*
0178       TF1 *gaussian = new TF1("gaussian","gaus",0,16);
0179       gaussian->SetParameters(fit[feech]->GetParameter(0),fit[feech]->GetParameter(1),fit[feech]->GetParameter(2));
0180       gaussian->SetLineColor(4);
0181       gaussian->Draw("same");
0182       */
0183     }
0184
0185   }
0186   else    // charge channel
0187   {
0188     if ( verbose )
0189     {
0190       g->Draw("ap");
0191       //gPad->SetLogy(1);
0192       fit[feech]->SetRange(-0.1,4.1);
0193       //fit[feech]->SetParameters(1.43*g->GetPointY(0),-3.2,3.8);
0194       fit[feech]->SetParameters(200.*g->GetPointY(0),-32.,9.8);
0195       g->Fit(fit[feech],"R");
0196
0197       Double_t chi2 = fit[feech]->GetChisquare();
0198       Double_t ndf = fit[feech]->GetNDF();
0199       cout << "chi2 ndf " << chi2 << "\t" << ndf << endl;
0200
0201     }
0202     else
0203     {
0204       fit[feech]->SetRange(-0.1,4.1);
0205       //fit[feech]->SetParameters(1.43*g->GetPointY(0),-3.2,3.8);
0206       fit[feech]->SetParameters(200.*g->GetPointY(0),-32.,9.8);
0207       g->Fit(fit[feech],"RNQ");
0208
0209       // save to tree
0210       _qch = feech;
0211       for (int ipar=0; ipar<3; ipar++)
0212       {
0213         _qfitpar[ipar] = static_cast<Float_t>( fit[feech]->GetParameter(ipar) );
0214       }
0215       Double_t chi2 = fit[feech]->GetChisquare();
0216       Double_t ndf = fit[feech]->GetNDF();
0217       _qchi2ndf = chi2/ndf;
0218
0219       _qfit_tree->Fill();
0220
0221     }
0222
0223     /*
0224     TF1 *expon = new TF1("expon","expo",0,16);
0225     expon->SetParameters(fit[feech]->GetParameter(0),fit[feech]->GetParameter(1));
0226     expon->SetLineColor(4);
0227     expon->Draw("same");
0228     */
0229
0230   }
0231
0232   if ( verbose )
0233   {
0234     fit[feech]->SetRange(0,16);
0235     fit[feech]->Draw("same");
0236
0237     gPad->Modified();
0238     gPad->Update();
0239     string junk;
0240     cin >> junk;
0241   }
0242 }
0243
0244
0245 //
0246 // get the events with only tails in them to use for fitting
0247 //
0248 void get_tails()
0249 {
0250   ifstream infile;
0251
0252   TString name;
0253
0254   name = "mbdsig_tails.root";
0255   TFile *savefile = new TFile(name,"RECREATE");
0256
0257   for (int ifeech=0; ifeech<NFEECH; ifeech++)
0258   {
0259     name = "h2_tail"; name += ifeech;
0260     h2_tail[ifeech]  = new TH2F(name,name,NSAMPLES,-0.5,NSAMPLES-0.5,2200,-0.1,1.1);
0261   }
0262
0263   // set up the sample number array
0264   double xsamp[NSAMPLES];
0265   double yval[NSAMPLES];
0266   for (int isamp=0; isamp<NSAMPLES; isamp++)
0267   {
0268     xsamp[isamp] = isamp;
0269   }
0270
0271
0272   string junk1, junk2;
0273   int ch;
0274   double mean;
0275   double pedsigma = 4.0;
0276
0277   for (int ifeech=0; ifeech<NFEECH; ifeech++)
0278   {
0279     name = "mbdsig"; name += ifeech; name += ".txt";
0280     infile.open( name.Data() );
0281     int evt = 1;
0282     while ( infile >> junk1 >> ch >> junk2 >> mean )
0283     {
0284       if ( evt%1000 == 1 )
0285       {
0286         cout << "evt " << evt << endl;
0287       }
0288
0289       for (int isamp=0; isamp<NSAMPLES; isamp++)
0290       {
0291         infile >> yval[isamp];
0292       }
0293
0294       // find those tail only events with high enough amplitude
0295       if ( (fabs(yval[12]-mean) < pedsigma*3) && ((yval[0]-mean) > 200) )
0296       {
0297         // fill histogram
0298         for (int isamp=0; isamp<NSAMPLES; isamp++)
0299         {
0300           h2_tail[ifeech]->Fill( isamp, (yval[isamp]-mean)/(yval[0]-mean) );
0301         }
0302       }
0303
0304       /*
0305       if ( fabs(yval[12]-mean) < pedsigma*3 )
0306       {
0307         g_subpulse->Draw("ap");
0308         gPad->Modified();
0309         gPad->Update();
0310         cin >> junk2;
0311       }
0312       */
0313
0314       evt++;
0315     }
0316
0317     infile.close();
0318   }
0319
0320   for (int ifeech=0; ifeech<NFEECH; ifeech++)
0321   {
0322     g_tail[ifeech] = find_th2ridge(h2_tail[ifeech]);
0323     //g_tail[ifeech]->InsertPointBefore(0,0,1);
0324     g_tail[ifeech]->SetPointError(0,0,g_tail[ifeech]->GetErrorY(1));
0325     name = "g_tail"; name += ifeech;
0326     g_tail[ifeech]->SetName(name);
0327
0328     h2_tail[ifeech]->Draw("colz");
0329     g_tail[ifeech]->Draw("cp");
0330
0331     name = "tail_"; name += ifeech; name += ".png";
0332     gPad->Print( name );
0333
0334     gPad->SetLogz(1);
0335     gPad->Modified();
0336     gPad->Update();
0337     /*
0338     string junk;
0339     cin >> junk;
0340     */
0341
0342     g_tail[ifeech]->Write();
0343   }
0344
0345   savefile->Write();
0346 }
0347
0348 //
0349 // fits s8 vs s0 in time channels to get time pileup correction
0350 // then writes out mbd_pileup.calib file
0351 //
0352 void fit_time_pileupcorr()
0353 {
0354   ofstream calibfile("mbd_pileup.calib");
0355
0356   TString name;
0357
0358   for (int ifeech=0; ifeech<NFEECH; ifeech++)
0359   {
0360     if ( ((ifeech/8)%2) == 0 )  // time ch
0361     {
0362       name = "f_s8s0_"; name += ifeech;
0363       fit_s8s0[ifeech] = new TF1(name,"[0]*x",0,4500);
0364       fit_s8s0[ifeech]->SetParameters(0.01);
0365
0366       h2_s8s0[ifeech]->Fit(fit_s8s0[ifeech],"R");
0367
0368       double p0 = fit_s8s0[ifeech]->GetParameter(0);
0369       double p0err = fit_s8s0[ifeech]->GetParError(0);
0370       double chi2ndf = fit_s8s0[ifeech]->GetChisquare()/fit_s8s0[ifeech]->GetNDF();
0371
0372       calibfile << ifeech
0373                 << "\t" << p0 << "\t" << p0err
0374                 << "\t" << 0. << "\t" << 0.
0375                 << "\t" << 0. << "\t" << 0.
0376                 << "\t" << chi2ndf << std::endl;
0377     }
0378     else  // charge ch
0379     {
0380       calibfile << ifeech
0381                 << "\t" << 200. << "\t" << 0.
0382                 << "\t" << -32. << "\t" << 0.
0383                 << "\t" << 9.8 << "\t" << 0.
0384                 << "\t" << 1. << std::endl;
0385     }
0386   }
0387
0388   calibfile.close();
0389 }
0390
0391
0392 //
0393 // get the residuals event by event
0394 //
0395 void calc_residuals()
0396 {
0397   ifstream infile;
0398
0399   TString name;
0400
0401   name = "mbdsig_tailresiduals.root";
0402   TFile *savefile = new TFile(name,"RECREATE");
0403
0404   for (int ifeech=0; ifeech<NFEECH; ifeech++)
0405   {
0406     name = "h2_residuals"; name += ifeech;
0407     h2_residuals[ifeech]  = new TH2F(name,name,NSAMPLES,-0.5,NSAMPLES-0.5,200,-100,100);
0408     h2_residuals[ifeech]->SetXTitle("sample");
0409
0410     name = "h2_s8s0_"; name += ifeech;
0411     h2_s8s0[ifeech]  = new TH2F(name,name,900,-0.5,4500-0.5,100,-20,120);
0412     h2_s8s0[ifeech]->SetXTitle("s0");
0413     h2_s8s0[ifeech]->SetYTitle("s8");
0414   }
0415
0416   _tfit_tree = new TTree("tfit","time ch fits");
0417   _tfit_tree->Branch("ch",&_tch,"ch/S");
0418   _tfit_tree->Branch("chi2ndf",&_tchi2ndf,"chi2ndf/F");
0419   _tfit_tree->Branch("ampl",&_tfitpar[0],"ampl/F");
0420   _tfit_tree->Branch("mean",&_tfitpar[1],"mean/F");
0421   _tfit_tree->Branch("sig",&_tfitpar[2],"sig/F");
0422   _tfit_tree->Branch("p0",&_tfitpar[3],"p0/F");
0423   _tfit_tree->Branch("p1",&_tfitpar[4],"p1/F");
0424   _tfit_tree->Branch("p2",&_tfitpar[5],"p2/F");
0425   _tfit_tree->Branch("s0",&_tfitpar[7],"s0/F");   // samp0
0426   _tfit_tree->Branch("s1",&_tfitpar[8],"s1/F");   // samp1
0427   _tfit_tree->Branch("s8",&_tfitpar[9],"s8/F");   // samp8 (used for timing)
0428
0429   _qfit_tree = new TTree("qfit","charge ch fits");
0430   _qfit_tree->Branch("ch",&_qch,"ch/S");
0431   _qfit_tree->Branch("chi2ndf",&_qchi2ndf,"chi2ndf/F");
0432   _qfit_tree->Branch("ampl",&_qfitpar[0],"ampl/F");
0433   _qfit_tree->Branch("mean",&_qfitpar[1],"mean/F");
0434   _qfit_tree->Branch("sig",&_qfitpar[2],"sig/F");
0435
0436   TGraphErrors *g_subpulse = new TGraphErrors(NSAMPLES);
0437
0438   // set up the sample number array
0439   double xsamp[NSAMPLES];
0440   double yval[NSAMPLES];
0441   for (int isamp=0; isamp<NSAMPLES; isamp++)
0442   {
0443     xsamp[isamp] = isamp;
0444   }
0445
0446   string junk1, junk2;
0447   int ch;
0448   double mean;
0449   double pedsigma = 4.0;
0450
0451   for (int ifeech=0; ifeech<NFEECH; ifeech++)
0452   {
0453     name = "mbdsig"; name += ifeech; name += ".txt";
0454     cout << "Processing " << name << endl;
0455     infile.open( name.Data() );
0456     int evt = 1;
0457     while ( infile >> junk1 >> ch >> junk2 >> mean )
0458     {
0459       if ( evt%1000 == 1 )
0460       {
0461         cout << "evt " << evt << "\tch " << ifeech << endl;
0462       }
0463
0464       for (int isamp=0; isamp<NSAMPLES; isamp++)
0465       {
0466         infile >> yval[isamp];
0467         g_subpulse->SetPoint( isamp, isamp, yval[isamp] - mean);
0468         g_subpulse->SetPointError( isamp, 0., 5.);
0469       }
0470
0471       // find those tail only events with high enough amplitude
0472       if ( (fabs(yval[12]-mean) < pedsigma*3) && ((yval[0]-mean) > 200) )
0473       {
0474         fit_tail( g_subpulse, ifeech );
0475
0476         // fill histogram
0477         for (int isamp=0; isamp<NSAMPLES; isamp++)
0478         {
0479           h2_residuals[ifeech]->Fill( isamp, (yval[isamp]-mean) - fit[ifeech]->Eval(isamp) );
0480           /*
0481           if ( fabs((yval[isamp]-mean) - fit[ifeech]->Eval(isamp)) > 100 )
0482           {
0483             cout << "xxx " << isamp << "\t" << yval[isamp]-mean << "\t" << fit[ifeech]->Eval(isamp) << endl;
0484             isbad = 1;
0485           }
0486           */
0487         }
0488       }
0489
0490       evt++;
0491     }
0492
0493     infile.close();
0494   }
0495
0496   fit_time_pileupcorr();
0497
0498   savefile->Write();
0499 }
0500
0501