BBCStudy/MBDLaser/MbdLaser.cc

0001 #include "MbdLaser.h"
0002
0003 #include <phool/phool.h>
0004 #include <phool/getClass.h>
0005 #include <phool/recoConsts.h>
0006 #include <fun4all/Fun4AllServer.h>
0007 #include <fun4all/Fun4AllReturnCodes.h>
0008 #include <ffaobjects/EventHeader.h>
0009
0010 #include <mbd/MbdDefs.h>
0011 #include <mbd/MbdOut.h>
0012 #include <mbd/MbdPmtContainer.h>
0013 #include <mbd/MbdPmtHit.h>
0014 #include <mbd/MbdGeom.h>
0015
0016 #include <calobase/TowerInfoContainer.h>
0017 #include <calobase/TowerInfo.h>
0018 #include <globalvertex/GlobalVertexMap.h>
0019 #include <globalvertex/GlobalVertex.h>
0020
0021 #include <TFile.h>
0022 #include <TTree.h>
0023 #include <TH1F.h>
0024 #include <TH2F.h>
0025 #include <TGraphErrors.h>
0026 #include <TF1.h>
0027 #include <TCanvas.h>
0028 #include <TString.h>
0029 #include <TLorentzVector.h>
0030 //#include <TMath.h>
0031 #include <TSystem.h>
0032
0033 #include <iostream>
0034 #include <cmath>
0035 #include <cstdio>
0036
0037 using namespace std;
0038 using namespace MbdDefs;
0039
0040 //---------------------------------------------------------------->
0041 MbdLaser::MbdLaser(const string &name) : SubsysReco(name)
0042 {
0043   _savefname = "mbdlaser.root";
0044 }
0045
0046 //----------------------------------------------------------------->
0047 int MbdLaser::Init(PHCompositeNode *topNode)
0048 {
0049   cout << PHWHERE << " Saving to file " << _savefname << endl;
0050   //PHTFileServer::get().open( _outfile, "RECREATE");
0051   _savefile = new TFile( _savefname.c_str(), "RECREATE" );
0052
0053   _tree = new TTree("t","MbdLaser");
0054   _tree->Branch("run",&f_run,"run/I");
0055   _tree->Branch("ch",&f_ch,"ch/I");
0056   _tree->Branch("qmean",&f_qmean,"qmean/F");
0057   _tree->Branch("qmerr",&f_qmerr,"qmerr/F");
0058   _tree->Branch("tmean",&f_tmean,"tmean/F");
0059   _tree->Branch("tmerr",&f_tmerr,"tmerr/F");
0060
0061   TString name, title;
0062
0063   //charge
0064   for (int ipmt=0; ipmt<128; ipmt++)
0065   {
0066     name = "h_mbdq"; name += ipmt;
0067     title = "mbd charge, ch "; title += ipmt;
0068     h_mbdq[ipmt] = new TH1F(name,title,1600,0,16000);
0069
0070
0071
0072     // TGraph to track mean of mbdq distribution
0073     name = "g_mbdq"; name += ipmt;
0074     title = "mbdq, ch "; title += ipmt;
0075     g_mbdq[ipmt] = new TGraphErrors();
0076     g_mbdq[ipmt]->SetName(name);
0077     g_mbdq[ipmt]->SetTitle(name);
0078
0079     /*
0080        name = "h_tdiff_ch"; name += ipmt;
0081        title = "tdiff, ch "; title += ipmt;
0082        h_tdiff_ch[ipmt] = new TH1F(name,title,600,-3,3);
0083        */
0084   }
0085   //Time
0086   for (int ipmt=0; ipmt<128; ipmt++)
0087   {
0088     //Time A for each ch filled with (t)
0089     name = "h_mbdt"; name += ipmt;
0090     title = "mbd Time, ch "; title += ipmt;
0091     h_mbdt[ipmt] = new TH1F(name,title,128,-0.5,127.5);
0092
0093
0094     // TGraph to track mean of mbdq distribution
0095     name = "g_mbdt"; name += ipmt;
0096     title = "mbdt, ch "; title += ipmt;
0097     g_mbdt[ipmt] = new TGraphErrors();
0098     g_mbdt[ipmt]->SetName(name);
0099     g_mbdt[ipmt]->SetTitle(name);
0100
0101
0102   }
0103
0104   for (int iarm=0; iarm<2; iarm++)
0105   {
0106     //
0107     name = "hevt_mbdt"; name += iarm;
0108     title = "mbd times, arm "; title += iarm;
0109     hevt_mbdt[iarm] = new TH1F(name,title,200,7.5,11.5);
0110     hevt_mbdt[iarm]->SetLineColor(4);
0111   }
0112
0113   // Time vs ch filled with (btt,pmt)
0114   h2_mbdbtt = new TH2F("h2_mbdtt","Time[ADC] vs Ch",128,-0.5,127.5,250,0.,25.);
0115   h2_mbdbtt->GetXaxis()->SetTitle("ch");
0116   h2_mbdbtt->GetYaxis()->SetTitle("Time");
0117
0118   //charge A vs ch D-2 filled with(bq,pmt)
0119   h2_mbdbq = new TH2F("h2_mbdbq","Charge[ADC] vs Ch", 128,-0.5,127.5,1600,0.,16000.);
0120   h2_mbdbq->GetXaxis()->SetTitle("ch");
0121   h2_mbdbq->GetYaxis()->SetTitle("Charge");
0122
0123
0124
0125
0126
0127   c_mbdt = new TCanvas("c_mbdt","MBD Times",1200,800);
0128   c_mbdt->Divide(1,2);
0129
0130   return 0;
0131 }
0132
0133 //------------------------------------------------------------------>
0134 int MbdLaser::InitRun(PHCompositeNode *topNode)
0135 {
0136   recoConsts *rc = recoConsts::instance();
0137   f_run = rc->get_IntFlag("RUNNUMBER");
0138
0139   GetNodes(topNode);
0140
0141   return 0;
0142 }
0143
0144 //-------------------------------------------------------------------->
0145 //Call user instructions for every event
0146 int MbdLaser::process_event(PHCompositeNode *topNode)
0147 {
0148   nprocessed++;
0149
0150   // skip the first 100 events
0151   // these are bad since the laser takes time to ramp up
0152   f_evt = _evtheader->get_EvtSequence();
0153   if ( f_evt<100 )
0154   {
0155     return Fun4AllReturnCodes::DISCARDEVENT;
0156   }
0157
0158   //if(f_evt%1==0) cout << PHWHERE << "Event " << f_evt << "\t" << ", nprocessed = " << nprocessed << endl;
0159
0160   // Initialize Variables
0161   f_mbdt[0] = -9999.;
0162   f_mbdt[1] = -9999.;
0163   f_mbdte[0] = -9999.;
0164   f_mbdte[1] = -9999.;
0165   f_mbdt0 = NAN;
0166   hevt_mbdt[0]->Reset();
0167   hevt_mbdt[1]->Reset();
0168
0169   for (int ipmt=0; ipmt<_mbdpmts->get_npmt(); ipmt++)
0170   {
0171     Float_t q = _mbdpmts->get_pmt(ipmt)->get_q();
0172     Float_t t = _mbdpmts->get_pmt(ipmt)->get_time();
0173     //Float_t phi = mbdgeom->get_phi(ipmt);   // get phi angle
0174
0175     h_mbdq[ipmt]->Fill( q );
0176     //cout << ipmt << ":\t" << q << "\t" << t << endl;
0177
0178     h_mbdt[ipmt]->Fill( t );
0179     // cout << ipmt << ":\t" << q << "\t" << t << endl;
0180
0181     // declear variables
0182     Float_t bq   = _mbdpmts->get_pmt(ipmt)->get_q();
0183     Float_t btt  = _mbdpmts->get_pmt(ipmt)->get_time();
0184     Short_t pmt =_mbdpmts->get_pmt(ipmt)->get_pmt();
0185     //cout << " bq = " << bq << ":\t" << "btt = " <<btt << ":\t" << "pmt = " <<pmt<<endl;
0186
0187
0188     //Fill 2D charge vs ch and Time vs ch
0189     h2_mbdbtt->Fill(pmt,btt);
0190     h2_mbdbq->Fill(pmt,bq);
0191
0192   }
0193
0194
0195   // LaserDST(topNode);
0196
0197   return 0;
0198 }
0199
0200 //-------------------------------------------------------------------->
0201 void MbdLaser::GetNodes(PHCompositeNode *topNode)
0202 {
0203   // Get the DST objects
0204
0205   _evtheader = findNode::getClass<EventHeader>(topNode, "EventHeader");
0206   if (!_evtheader && f_evt<10) cout << PHWHERE << " EventHeader node not found on node tree" << endl;
0207
0208   // MbdOut
0209   _mbdout = findNode::getClass<MbdOut>(topNode, "MbdOut");
0210   if(!_mbdout && f_evt<4) cout << PHWHERE << " MbdOut node not found on node tree" << endl;
0211
0212   // MbdPmt Info
0213   _mbdpmts = findNode::getClass<MbdPmtContainer>(topNode, "MbdPmtContainer");
0214   if(!_mbdpmts && f_evt<4) cout << PHWHERE << " MbdPmtContainer node not found on node tree" << endl;
0215
0216 }
0217
0218 //------------------------------------------------------------------->
0219 int MbdLaser::End(PHCompositeNode *topNode)
0220 {
0221   _savefile->cd();
0222
0223   //Filling g_mbdq and g_mbdt
0224   for (int ipmt=0; ipmt<MBD_N_PMT; ipmt++)
0225   {
0226
0227     // Fill info on q distribution and t
0228     f_ch = ipmt;
0229     f_qmean = h_mbdq[ipmt]->GetMean();
0230     f_qmerr = h_mbdq[ipmt]->GetMeanError();
0231     f_tmean = h_mbdt[ipmt]->GetMean();
0232     f_tmerr = h_mbdt[ipmt]->GetMeanError();
0233
0234     _tree->Fill();
0235
0236     //cout << "run =  "<< f_run << "\t" << " pmt =  "<< f_ch << "\t" <<" qmean =  "<< f_qmean << "\t" << " qmerr =  "<<f_qmerr << "\t"<< "tmean =  "<< f_tmean << "\t" << "tmerr = " << f_tmerr<< endl;;
0237
0238     g_mbdq[ipmt]->SetPoint(0,f_run,f_qmean);
0239     g_mbdq[ipmt]->SetPointError(0,0,f_qmerr);
0240     g_mbdq[ipmt]->Write();
0241
0242     g_mbdt[ipmt]->SetPoint(0,f_run,f_tmean);
0243     g_mbdt[ipmt]->SetPointError(0,0,f_tmerr);
0244     g_mbdt[ipmt]->Write();
0245
0246   }
0247
0248   _savefile->Write();
0249   _savefile->Close();
0250
0251   return 0;
0252 }
0253
0254 void MbdLaser::LaserDST(PHCompositeNode *topNode)
0255 {
0256
0257   //Float_t r = (4.4049/4.05882);
0258   Float_t r = 1.0;
0259
0260   // Fill info from each PMT
0261   for (int ipmt=0; ipmt<_mbdpmts->get_npmt(); ipmt++)
0262   {
0263     Float_t q = _mbdpmts->get_pmt(ipmt)->get_q();
0264     Float_t t = _mbdpmts->get_pmt(ipmt)->get_time();
0265     //Float_t phi = mbdgeom->get_phi(ipmt);   // get phi angle
0266
0267     if ( fabs(t) < 25. )
0268     {
0269       h_mbdq[ipmt]->Fill( q*r );
0270       h_mbdt[ipmt]->Fill(t);
0271       // cout << ipmt << ":\t" << q << "\t" << t << endl;
0272     }
0273
0274     // cout << ipmt << ":\t" << q << "\t" << t << endl;
0275   }
0276
0277 }
0278
0279 int MbdLaser::Getpeaktime(TH1 * h)
0280 {
0281   int getmaxtime, tcut = -1;
0282
0283   for(int bin = 1; bin < h->GetNbinsX()+1; bin++)
0284   {
0285     double c = h->GetBinContent(bin);
0286     double max = h->GetMaximum();
0287     int bincenter = h->GetBinCenter(bin);
0288     if(max == c)
0289     {
0290       getmaxtime = bincenter;
0291       if(getmaxtime != -1) tcut = getmaxtime;
0292     }
0293   }
0294
0295   return tcut;
0296
0297 }
0298
0299