BBCStudy/Mwibanda/MBDStudy.cc

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