analysis/BBCStudy/BBCStudy.cc

0001 #include "BBCStudy.h"
0002
0003 #include <phool/phool.h>
0004 #include <phool/getClass.h>
0005 #include <g4main/PHG4HitContainer.h>
0006 #include <g4main/PHG4TruthInfoContainer.h>
0007 #include <g4main/PHG4Particle.h>
0008 #include <g4main/PHG4Hit.h>
0009 #include <g4main/PHG4VtxPoint.h>
0010 #include <fun4all/PHTFileServer.h>
0011 #include <fun4all/Fun4AllServer.h>
0012 #include <ffaobjects/EventHeaderv1.h>
0013 #include <mbd/MbdOut.h>
0014 #include <mbd/MbdPmtContainer.h>
0015 #include <mbd/MbdPmtHit.h>
0016 #include <mbd/MbdGeom.h>
0017 #include <mbd/MbdDefs.h>
0018
0019 #include <TFile.h>
0020 #include <TTree.h>
0021 #include <TH1F.h>
0022 #include <TH2F.h>
0023 #include <TF1.h>
0024 #include <TCanvas.h>
0025 #include <TString.h>
0026 #include <TLorentzVector.h>
0027 //#include <TMath.h>
0028 #include <TDatabasePDG.h>
0029 #include <TRandom3.h>
0030 #include <TSystem.h>
0031
0032 #include <iostream>
0033 #include <cmath>
0034
0035 using namespace std;
0036 using namespace MbdDefs;
0037
0038 //const Double_t C = 29.9792458; // cm/ns
0039
0040 //____________________________________
0041 BBCStudy::BBCStudy(const string &name) : SubsysReco(name),
0042     _tree( 0 ),
0043     f_evt( 0 ),
0044     _pdg( 0 ),
0045     _rndm( 0 ),
0046     _tres( 0.05 ),
0047     _savefname( "bbcstudy.root" ),
0048     _savefile( 0 )
0049 {
0050
0051 }
0052
0053 //___________________________________
0054 int BBCStudy::Init(PHCompositeNode *topNode)
0055 {
0056   cout << PHWHERE << " Saving to file " << _savefname << endl;
0057   //PHTFileServer::get().open( _outfile, "RECREATE");
0058   _savefile = new TFile( _savefname.c_str(), "RECREATE" );
0059
0060   // Global BBC variables
0061   _tree = new TTree("t","BBC Study");
0062   _tree->Branch("evt",&f_evt,"evt/I");
0063   _tree->Branch("bimp",&f_bimp,"bimp/F");
0064   _tree->Branch("ncoll",&f_ncoll,"ncoll/I");
0065   _tree->Branch("npart",&f_npart,"npart/I");
0066   _tree->Branch("vx",&f_vx,"vx/F");
0067   _tree->Branch("vy",&f_vy,"vy/F");
0068   _tree->Branch("vz",&f_vz,"vz/F");
0069   _tree->Branch("bns",&f_bbcn[0],"bns/S");
0070   _tree->Branch("bnn",&f_bbcn[1],"bnn/S");
0071   _tree->Branch("bqs",&f_bbcq[0],"bqs/F");
0072   _tree->Branch("bqn",&f_bbcq[1],"bqn/F");
0073   _tree->Branch("bts",&f_bbct[0],"bts/F");
0074   _tree->Branch("btn",&f_bbct[1],"btn/F");
0075   _tree->Branch("btes",&f_bbcte[0],"btes/F");
0076   _tree->Branch("bten",&f_bbcte[1],"bten/F");
0077   _tree->Branch("bz",&f_bbcz,"bz/F");
0078   _tree->Branch("bt0",&f_bbct0,"bt0/F");
0079
0080   _pdg = TDatabasePDG::Instance();  // database of PDG info on particles
0081   _rndm = new TRandom3(0);
0082
0083   TString name, title;
0084   for (int ipmt=0; ipmt<MbdDefs::MBD_N_PMT; ipmt++)
0085   {
0086     name = "h_bbcq"; name += ipmt;
0087     title = "bbc charge, ch "; title += ipmt;
0088     h_bbcq[ipmt] = new TH1F(name,title,1200,0,120*90);
0089
0090     /*
0091     name = "h_tdiff_ch"; name += ipmt;
0092     title = "tdiff, ch "; title += ipmt;
0093     h_tdiff_ch[ipmt] = new TH1F(name,title,600,-3,3);
0094     */
0095   }
0096
0097   for (int iarm=0; iarm<2; iarm++)
0098   {
0099     name = "h_bbcqtot"; name += iarm;
0100     title = "bbc charge, arm "; title += iarm;
0101     h_bbcqtot[iarm] = new TH1F(name,title,1200,0,120*60); // npe for 1 mip = 120, and up to 30 mips are possible
0102
0103     //
0104     name = "hevt_bbct"; name += iarm;
0105     title = "bbc times, arm "; title += iarm;
0106     hevt_bbct[iarm] = new TH1F(name,title,200,7.5,11.5);
0107     hevt_bbct[iarm]->SetLineColor(4);
0108   }
0109   h2_bbcqtot = new TH2F("h2_bbcqtot","north BBCQ vs South BBCQ",300,0,120*900,300,0,120*900);
0110
0111   h_ztrue = new TH1F("h_ztrue","true z-vtx",600,-30,30);
0112   h_tdiff = new TH1F("h_tdiff","dt (measured - true_time)",6000,-3,3);
0113   h2_tdiff_ch = new TH2F("h2_tdiff_ch","dt (measured - true time) vs ch",MBD_N_PMT,-0.5,MBD_N_PMT-0.5,200,-2,2);
0114
0115   gaussian = new TF1("gaussian","gaus",0,20);
0116   gaussian->FixParameter(2,0.05);   // set sigma to 50 ps
0117
0118   c_bbct = new TCanvas("c_bbct","BBC Times",1200,800);
0119   c_bbct->Divide(1,2);
0120
0121   return 0;
0122 }
0123
0124 //___________________________________
0125 int BBCStudy::InitRun(PHCompositeNode *topNode)
0126 {
0127   _bbcgeom = new MbdGeom();
0128   GetNodes(topNode);
0129
0130   return 0;
0131 }
0132
0133 //__________________________________
0134 //Call user instructions for every event
0135 int BBCStudy::process_event(PHCompositeNode *topNode)
0136 {
0137   nprocessed++;
0138
0139   //GetNodes(topNode);
0140
0141   f_evt = _evtheader->get_EvtSequence();
0142   if(f_evt%1==0) cout << PHWHERE << "Event " << f_evt << "\t" << ", nprocessed = " << nprocessed << endl;
0143
0144   // Initialize Variables
0145   f_bbcn[0] = 0;
0146   f_bbcn[1] = 0;
0147   f_bbcq[0] = 0.;
0148   f_bbcq[1] = 0.;
0149   f_bbct[0] = -9999.;
0150   f_bbct[1] = -9999.;
0151   f_bbcte[0] = -9999.;
0152   f_bbcte[1] = -9999.;
0153   f_bbcz = NAN;
0154   f_bbct0 = NAN;
0155   hevt_bbct[0]->Reset();
0156   hevt_bbct[1]->Reset();
0157
0158   // Get Truth Centrality Info
0159   f_bimp = _evtheader->get_ImpactParameter();
0160   f_ncoll = _evtheader->get_ncoll();
0161   f_npart = _evtheader->get_npart();
0162
0163   // Get Primaries AND Secondaries
0164   /*
0165   PHG4TruthInfoContainer::ConstRange primary_range = _truth_container->GetPrimaryParticleRange();
0166   int nprimaries = 0;
0167
0168   for (auto part_iter = primary_range.first; part_iter != primary_range.second; ++part_iter)
0169   {
0170     PHG4Particle *particle = part_iter->second;
0171     Float_t e = particle->get_e();
0172     Float_t px = particle->get_px();
0173     Float_t py = particle->get_py();
0174     Float_t pz = particle->get_pz();
0175     Float_t eta = 0.5*log((particle->get_e()+particle->get_pz())/ (particle->get_e()-particle->get_pz()));
0176     Float_t pt = sqrt( px*px + py*py );
0177     Float_t phi = atan2(particle->get_py(), particle->get_px());
0178     Float_t pid = particle->get_pid();
0179
0180     nprimaries++;
0181   }
0182   if ( f_evt < 20 )
0183   {
0184     cout << "Num primaries = " << nprimaries << "\t" << _truth_container->GetNumPrimaryVertexParticles() << endl;
0185   }
0186   */
0187
0188
0189   // Get True Vertex
0190   // NB: Currently PrimaryVertexIndex is always 1, need to figure out how to handle collision pile-up
0191   PHG4VtxPoint *vtxp = _truth_container->GetPrimaryVtx( _truth_container->GetPrimaryVertexIndex() );
0192   if ( vtxp != 0 )
0193   {
0194     f_vx = vtxp->get_x();
0195     f_vy = vtxp->get_y();
0196     f_vz = vtxp->get_z();
0197     f_vt = vtxp->get_t();
0198
0199     if ( f_evt<20 )
0200     {
0201       cout << "VTXP " << "\t" << f_vx << "\t" << f_vy << "\t" << f_vz << "\t" << f_vt << endl;
0202     }
0203
0204     h_ztrue->Fill( f_vz );
0205
0206   }
0207
0208
0209   // Get All Vertices
0210   PHG4TruthInfoContainer::ConstVtxRange vtx_range = _truth_container->GetVtxRange();
0211   unsigned int nvtx = 0;
0212
0213   /*
0214   for (auto vtx_iter = vtx_range.first; vtx_iter != vtx_range.second; ++vtx_iter)
0215   {
0216     PHG4VtxPoint *vtx = vtx_iter->second;
0217     double vx = vtx->get_x();
0218     double vy = vtx->get_y();
0219     double vz = vtx->get_z();
0220     double vt = vtx->get_t();
0221
0222     // look at first few tracks
0223     // What are negative track id's?
0224     if ( abs(vtx->get_id()) < 8 && f_evt<20 )
0225     {
0226       cout << "vtx " << nvtx << "\t" << vtx->get_id() << "\t" << vt << "\t" << vx << "\t" << vy << "\t" << vz << endl;
0227     }
0228     nvtx++;
0229   }
0230   */
0231
0232   if ( f_evt<20 )
0233   {
0234     cout << "Num Vertices = " << nvtx << "\t" << _truth_container->GetNumVertices() << endl;
0235     //cout << "Primary Vertex = " << _truth_container->GetPrimaryVertexIndex() << endl;
0236   }
0237
0238   // Go through BBC hits to see what they look like
0239
0240   float len[MbdDefs::MBD_N_PMT] = {0.};
0241   float edep[MbdDefs::MBD_N_PMT] = {0.};
0242   float first_time[MbdDefs::MBD_N_PMT];    // First hit time for each tube
0243   std::fill_n(first_time, MbdDefs::MBD_N_PMT, 1e12);
0244
0245
0246   unsigned int nhits = 0;
0247
0248   TLorentzVector v4;
0249   //TLorentzVector hitpos;
0250
0251   PHG4HitContainer::ConstRange bbc_hit_range = _bbchits->getHits();
0252   for (auto hit_iter = bbc_hit_range.first; hit_iter != bbc_hit_range.second; hit_iter++)
0253   {
0254     PHG4Hit *this_hit = hit_iter->second;
0255
0256     unsigned int ipmt = this_hit->get_layer();  // pmt channel
0257     int arm = ipmt/64;;                         // south=0, north=1
0258
0259     int trkid = this_hit->get_trkid();
0260     if ( trkid>0 && f_evt<20 ) cout << "TRKID " << trkid << endl;
0261
0262     PHG4Particle *part = _truth_container->GetParticle( trkid );
0263     v4.SetPxPyPzE( part->get_px(), part->get_py(), part->get_pz(), part->get_e() );
0264
0265     int pid = part->get_pid();
0266     TParticlePDG *partinfo = _pdg->GetParticle( pid );
0267     Double_t charge = -9999.;
0268     if ( partinfo )
0269     {
0270       charge = partinfo->Charge()/3;  // PDG gives charge in 1/3 e
0271     }
0272     else if ( part->isIon() )
0273     {
0274       charge = part->get_IonCharge();
0275     }
0276
0277     Double_t beta = v4.Beta();
0278
0279     // Determine time
0280     //hitpos.SetXYZT( this_hit->get_x(), this_hit->get_y(), this_hit->get_z(), this_hit->get_t() );
0281
0282     float xsign = 1.;
0283     float zsign = -1.;
0284     if ( arm == 1 )
0285     {
0286       xsign = -1.;
0287       zsign = 1.;
0288     }
0289
0290     float tube_x = _bbcgeom->get_x(ipmt);
0291     float tube_y = _bbcgeom->get_y(ipmt);
0292     float tube_z = zsign*253.;
0293     float flight_z = fabs(tube_z - vtxp->get_z());
0294
0295     float flight_time = sqrt( tube_x*tube_x + tube_y*tube_y + flight_z*flight_z )/C;
0296     float tdiff = this_hit->get_t(1) - ( vtxp->get_t() + flight_time );
0297
0298     // get the first time
0299     if ( this_hit->get_t(1) < first_time[ipmt] )
0300     {
0301       if ( fabs( this_hit->get_t(1) ) < 106.5 )
0302       {
0303         first_time[ipmt] = this_hit->get_t(1) - vtxp->get_t();
0304         Float_t dt = static_cast<float>( _rndm->Gaus( 0, _tres ) ); // get fluctuation in time
0305         first_time[ipmt] += dt;
0306       }
0307       else
0308       {
0309         cout << "BAD " << ipmt << "\t" << this_hit->get_t(1) << endl;
0310       }
0311     }
0312
0313     if ( f_evt<10 )
0314     {
0315       cout << "hit " << ipmt << "\t" << trkid << "\t" << pid
0316         //<< "\t" << v4.M()
0317         << "\t" << beta
0318         << "\t" << this_hit->get_path_length()
0319         << "\t" << this_hit->get_edep()
0320         //<< "\t" << v4.Eta()
0321         //<< "\t" << v4.Pt()
0322         //<< "\t" << v4.P()
0323         << "\t" << this_hit->get_x(1)
0324         << "\t" << this_hit->get_y(1)
0325         << "\t" << this_hit->get_z(1)
0326         << "\t" << this_hit->get_t(1)
0327         << "\t" << tdiff
0328         << endl;
0329
0330       //cout << "WWW " << first_time[ipmt] << endl;
0331     }
0332
0333     edep[ipmt] += this_hit->get_edep();
0334
0335     // get summed path length for particles that can create CKOV light
0336     // n.p.e. is determined from path length
0337     if ( beta > MbdDefs::v_ckov && charge != 0. )
0338     {
0339       len[ipmt] += this_hit->get_path_length();
0340
0341       if ( trkid>0 )
0342       {
0343         h_tdiff->Fill( tdiff );
0344         h2_tdiff_ch->Fill( ipmt, tdiff );
0345       }
0346
0347       _pids[pid] += 1;
0348     }
0349
0350     // sanity check
0351     if ( part->get_track_id() != trkid )
0352     {
0353       cout << "ERROR " << endl;
0354     }
0355
0356     nhits++;
0357   }
0358
0359   if ( f_evt<20 )
0360   {
0361     cout << "******** " << f_evt << " **************" << endl;
0362   }
0363
0364
0365   // process the data from each channel
0366   for (int iarm=0; iarm<2; iarm++)
0367   {
0368     f_bbct[iarm] = 0.;
0369   }
0370
0371   vector<float> hit_times[2];   // times of the hits in each [arm]
0372
0373   for (int ipmt=0; ipmt<MbdDefs::MBD_N_PMT; ipmt++)
0374   {
0375     //cout << "ZZZ " << ipmt << "\t" << first_time[ipmt] << "\t" << edep[ipmt] << "\t" << len[ipmt] << endl;
0376
0377     int arm = _bbcgeom->get_arm(ipmt); // ch 0-63 = south, ch 64-127 = north
0378
0379     // Fill charge and time info
0380     if ( len[ipmt]>0. )
0381     {
0382       if ( f_evt<20 )
0383       {
0384         cout << "ipmt " << ipmt << "\t" << len[ipmt] << "\t" << edep[ipmt] << endl;
0385       }
0386
0387       // Get charge in BBC tube
0388       float npe = len[ipmt]*(120/3.0);  // we get 120 p.e. per 3 cm
0389       float dnpe = static_cast<float>( _rndm->Gaus( 0, sqrt(npe) ) ); // get fluctuation in npe
0390
0391       npe += dnpe;  // apply the fluctuations in npe
0392
0393       f_bbcq[arm] += npe;
0394
0395       h_bbcq[ipmt]->Fill( npe );
0396       h_bbcqtot[arm]->Fill( npe );
0397
0398       // Now time
0399       if ( first_time[ipmt] < 9999. )
0400       {
0401         // Fill evt histogram
0402         hevt_bbct[arm]->Fill( first_time[ipmt] );
0403         hit_times[arm].push_back( first_time[ipmt] );
0404
0405         f_bbct[arm] += first_time[ipmt];
0406         //cout << "XXX " << ipmt << "\t" << f_bbct[arm] << "\t" << first_time[ipmt] << endl;
0407
0408       }
0409       else  // should never happen
0410       {
0411         cout << "ERROR, have hit but no time" << endl;
0412       }
0413
0414       // threshold should be > 0.
0415       ++f_bbcn[arm];
0416     }
0417   }
0418
0419   // Get best t
0420   if ( f_bbcn[0]>0 && f_bbcn[1]> 0 )
0421   {
0422     for (int iarm=0; iarm<2; iarm++)
0423     {
0424       c_bbct->cd(iarm+1);
0425
0426       std::sort( hit_times[iarm].begin(), hit_times[iarm].end() );
0427       float earliest = hit_times[iarm][0];
0428
0429       gaussian->SetParameter(0,5);
0430       gaussian->SetParameter(1, earliest);
0431       gaussian->SetRange(6, earliest+ 5*0.05);
0432       //gaussian->SetParameter(1,hevt_bbct[iarm]->GetMean());
0433       //gaussian->SetRange(6,hevt_bbct[iarm]->GetMean()+0.125);
0434
0435       hevt_bbct[iarm]->Fit(gaussian,"BLR");
0436       hevt_bbct[iarm]->Draw();
0437
0438       if ( f_bbcn[iarm] > 0 )
0439       {
0440         //f_bbct[iarm] = f_bbct[iarm] / f_bbcn[iarm];
0441         f_bbct[iarm] = gaussian->GetParameter(1);
0442         f_bbcte[iarm] = earliest;
0443       }
0444     }
0445
0446     // Now calculate zvtx, t0 from best times
0447     f_bbcz = (f_bbct[0] - f_bbct[1])*C/2.0;
0448     f_bbct0 = (f_bbct[0] + f_bbct[1])/2.0;
0449
0450   }
0451
0452   c_bbct->Modified();
0453   c_bbct->Update();
0454
0455   if ( fabs(f_bbcz-f_vz)>5.0 && f_bbcn[0]>1 && f_bbcn[1]>1 )
0456   {
0457     string response;
0458     cout << "Event " << f_evt << " ? ";
0459     //cin >> response;
0460     //if ( response[0] == 'q' )
0461     TString name = "evt_"; name += f_evt; name += ".png";
0462     c_bbct->SaveAs( name );
0463   }
0464
0465   h2_bbcqtot->Fill( f_bbcq[0], f_bbcq[1] );
0466
0467   _tree->Fill();
0468
0469   //CheckDST(topNode);
0470
0471   return 0;
0472 }
0473
0474 //___________________________________
0475 void BBCStudy::GetNodes(PHCompositeNode *topNode)
0476 {
0477   // Get the DST objects
0478
0479   // Truth container
0480   _truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0481   if(!_truth_container && f_evt<10) cout << PHWHERE << " PHG4TruthInfoContainer node not found on node tree" << endl;
0482
0483   // BBC hit container
0484   _bbchits = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_BBC");
0485   if(!_bbchits && f_evt<10) cout << PHWHERE << " G4HIT_BBC node not found on node tree" << endl;
0486
0487   _evtheader = findNode::getClass<EventHeaderv1>(topNode, "EventHeader");
0488   if (!_evtheader && f_evt<10) cout << PHWHERE << " G4HIT_BBC node not found on node tree" << endl;
0489
0490 }
0491
0492 //______________________________________
0493 int BBCStudy::End(PHCompositeNode *topNode)
0494 {
0495   _savefile->cd();
0496   _savefile->Write();
0497   _savefile->Close();
0498
0499   // print out list of pids that hit BBC
0500   cout << "PIDs of Particles that hit BBC" << endl;
0501   unsigned int ipid = 0;
0502   double npart = 0;
0503   for (auto & pid : _pids)
0504   {
0505     npart += pid.second;
0506   }
0507   for (auto & pid : _pids)
0508   {
0509     cout << pid.first << "\t" << pid.second << "\t" << pid.second/npart << endl;
0510     ipid++;
0511   }
0512   cout << "There were " << ipid << " different particle types" << endl;
0513
0514   return 0;
0515 }
0516
0517 void BBCStudy::CheckDST(PHCompositeNode *topNode)
0518 {
0519   // BbcOut
0520   _bbcout = findNode::getClass<MbdOut>(topNode, "MbdOut");
0521   if(!_bbcout && f_evt<4) cout << PHWHERE << " MbdOut node not found on node tree" << endl;
0522
0523   // BbcOut
0524   _bbcpmts = findNode::getClass<MbdPmtContainer>(topNode, "MbdPmtContainer");
0525   if(!_bbcpmts && f_evt<4) cout << PHWHERE << " MbdPmtContainer node not found on node tree" << endl;
0526
0527   Float_t bbcz = _bbcout->get_zvtx();
0528   if ( f_bbcz != bbcz )
0529   {
0530     cout << "ERROR, f_bbcz != bbcz, " << f_bbcz << "\t" << bbcz << endl;
0531   }
0532
0533   for (int iarm=0; iarm<2; iarm++)
0534   {
0535     if ( f_bbcq[iarm] != _bbcout->get_q(iarm) )
0536     {
0537       cout << "ERROR, f_bbcq != bbcq, arm " << iarm << "\t" << f_bbcq[iarm] << "\t" << _bbcout->get_q(iarm) << endl;
0538     }
0539   }
0540
0541 }