timing/src/InttStreamingTiming.cc

0001 #include "InttStreamingTiming.h"
0002
0003 using namespace std;
0004
0005 //____________________________________________________________________________..
0006 InttStreamingTiming::InttStreamingTiming(const std::string &name) :
0007   SubsysReco(name)
0008 {
0009
0010 }
0011
0012 //____________________________________________________________________________..
0013 InttStreamingTiming::~InttStreamingTiming()
0014 {
0015
0016 }
0017
0018 /////////////////////////////////////////////////////////////////////////
0019 // private
0020 /////////////////////////////////////////////////////////////////////////
0021 int InttStreamingTiming::GetNodes(PHCompositeNode *topNode)
0022 {
0023
0024   /////////////////////////////////////////////////////////////////////////
0025   // GL1 raw hit node
0026   string gl1_raw_node_name = "GL1RAWHIT";
0027   gl1_ = findNode::getClass<Gl1Packetv2>( topNode, gl1_raw_node_name);
0028   //Gl1Packet* gl1  = findNode::getClass<Gl1Packet>(topNode, "GL1RAWHIT");
0029
0030   if( !gl1_ )
0031     {
0032
0033       cerr << PHWHERE << gl1_raw_node_name << " node is missing." << endl;
0034       return Fun4AllReturnCodes::ABORTEVENT;
0035     }
0036
0037   /////////////////////////////////////////////////////////////////////////
0038   // INTT raw hit
0039   string node_name_inttrawhit = "INTTRAWHIT";
0040   node_inttrawhit_map_ =
0041     findNode::getClass<InttRawHitContainer>(topNode, node_name_inttrawhit);
0042
0043   if (!node_inttrawhit_map_)
0044     {
0045       cerr << PHWHERE << node_name_inttrawhit << " node is missing." << endl;
0046       return Fun4AllReturnCodes::ABORTEVENT;
0047     }
0048
0049   /////////////////////////////////////////////////////////////////////////
0050   // Trkr hit
0051   // TRKR_HITSET (IO,TrkrHitSetContainerv1)
0052
0053   string node_name_trkr_hitset = "TRKR_HITSET";
0054   // node_inttrawhit_map_ =
0055   node_trkrhitset_map_ =
0056     findNode::getClass<TrkrHitSetContainer>(topNode, node_name_trkr_hitset);
0057
0058   if (!node_trkrhitset_map_)
0059     {
0060       cerr << PHWHERE << node_name_trkr_hitset << " node is missing." << endl;
0061       return Fun4AllReturnCodes::ABORTEVENT;
0062     }
0063
0064   return 0;
0065 }
0066
0067 vector < std::pair < uint16_t, int > > InttStreamingTiming::GetBcoEventCounter()
0068 {
0069   //cout << "vector < std::pair < uint16_t, int > > InttStreamingTiming::GetBcoEventCounter()" << endl;
0070   vector < pair < InttNameSpace::Online_s,  unsigned int > > online_hits = this->GetHits(  node_trkrhitset_map_->getHitSets() );
0071
0072   vector < InttRawHit* > raw_hits = this->GetRawHitsWithoutClone();
0073   vector < std::pair < uint16_t, int > > rtn;
0074
0075   int counter = 0;
0076   for( int i=0; i<online_hits.size(); i++ )
0077     {
0078       ////////////////////////////////////////////
0079       // struct Online_s                        //
0080       //   int lyr, int ldr,  int arm           //
0081       //   int chp, int chn                     //
0082       ////////////////////////////////////////////
0083       auto online_hit = online_hits[i].first;
0084       auto dac = online_hits[i].second;
0085
0086       ////////////////////////////////////////////
0087       // struct RawData_s:                      //
0088       //   int felix_server, int felix_channel  //
0089       //   int chip,         int channel,       //
0090       ////////////////////////////////////////////
0091       // InttNameSpace::RawData_s
0092       auto raw_data  = InttNameSpace::ToRawData( online_hit );
0093
0094       ////////////////////////////////////////////
0095       // struct Offline_s                       //
0096       //   int layer,     int ladder_phi,       //
0097       //   int ladder_z,  int strip_x,          //
0098       //   int strip_y,                         //
0099       ////////////////////////////////////////////
0100       // InttNameSpace::Offline_s
0101       auto offline_hit = InttNameSpace::ToOffline( online_hit );
0102
0103       vector < int > used_index( raw_hits.size() );
0104
0105       for( int j=0; j<raw_hits.size(); j++ )
0106     {
0107
0108       auto raw_hit = raw_hits[j];
0109
0110       int chan_raw      = raw_hit->get_channel_id();
0111       // first selections
0112       if( chan_raw != raw_data.channel )
0113         continue;
0114
0115       // uint16_t InttRawHit::get_chip_id
0116       int chip_raw      = raw_hit->get_chip_id();
0117       if(chip_raw  > InttQa::kChip_num)
0118         chip_raw = chip_raw - InttQa::kChip_num;
0119
0120       chip_raw--; // to change from 1-base to 0-base
0121       if( chip_raw != raw_data.chip )
0122         continue;
0123
0124       int felix_raw     = raw_hit->get_packetid() - InttQa::kFirst_pid;
0125       int felix_ch_raw  = raw_hit->get_fee();
0126       bool is_same = (felix_raw == raw_data.felix_server
0127               || felix_ch_raw == raw_data.felix_channel );
0128               // || adc_raw == adc );
0129
0130       if( is_same == false )
0131         continue;
0132
0133       auto bco_raw      = raw_hit->get_FPHX_BCO();
0134       int event_counter_raw = raw_hit->get_event_counter(); // uint32_t IttRawHit::get_event_counter()
0135       pair < uint16_t, int > bco_event_counter( bco_raw, event_counter_raw );
0136       rtn.push_back( bco_event_counter );
0137       counter++;
0138       //used_index.push_back( j );
0139       //raw_hits.erase( raw_hits.begin() + j );
0140       //break;
0141     }
0142       //      cout << endl;
0143
0144     }
0145
0146
0147   // cout << counter << " / " << online_hits.size() << " found their corresponding raw hit while there are "
0148   //      << raw_hits.size() << " raw hits in this event."
0149   //      << endl;
0150   return rtn;
0151 }
0152
0153 std::vector < int > InttStreamingTiming::GetTriggerBits()
0154 {
0155
0156   //uint64_t trigger_vector = gl1_->getScaledVector();
0157   uint64_t trigger_vector = gl1_->getLiveVector();
0158
0159   vector < int > rtn;
0160   while( trigger_vector != 0 )
0161     {
0162       int this_bit = 0 ;
0163       this_bit = trigger_vector & 1;
0164       // cout << std::bitset<32>(trigger_vector) << "  "
0165       //       << this_bit << "\t";
0166
0167       trigger_vector = trigger_vector >> 1;
0168
0169       //cout << std::bitset<32>(trigger_vector) << endl;
0170
0171       rtn.push_back( this_bit );
0172     }
0173
0174
0175   // for( int i=0; i<rtn.size() ; i++ )
0176   //   {
0177   //     if( i == 6 || i == 7 )
0178   //    continue;
0179
0180   //     if( rtn[i] != 0 )
0181   //    {
0182   //      cout << trigger_names_[ i ] << ",\t";
0183   //    }
0184   //   }
0185   // cout << endl;
0186
0187   return rtn;
0188 }
0189
0190 bool InttStreamingTiming::IsSame( InttRawHit* hit1, InttRawHit* hit2 )
0191 {
0192
0193   if( hit1->get_word() != hit2->get_word() ) // only this might be enough
0194     return false;
0195   else if( hit1->get_bco() != hit2->get_bco() )
0196     return false;
0197   else if( hit1->get_FPHX_BCO() != hit2->get_FPHX_BCO() )
0198     return false;
0199   else if( hit1->get_channel_id() != hit2->get_channel_id() )
0200     return false;
0201   else if( hit1->get_chip_id() != hit2->get_chip_id() )
0202     return false;
0203   else if( hit1->get_fee() != hit2->get_fee() )
0204     return false;
0205   else if( hit1->get_packetid() != hit2->get_packetid() )
0206     return false;
0207   else if( hit1->get_adc() != hit2->get_adc() )
0208     return false;
0209   else if( hit1->get_full_FPHX() != hit2->get_full_FPHX() )
0210     return false;
0211   // else if( hit1->get_full_ROC() != hit2->get_full_ROC() ) // maybe no value is assigned?
0212   //   return false;
0213   // else if( hit1->get_amplitude()!= hit2->get_amplitude()) // maybe no value is assigned?
0214   //   return false;
0215
0216   // hits pass all comparison are the same
0217   return true;
0218 };
0219
0220 vector < InttRawHit* > InttStreamingTiming::GetRawHits()
0221 {
0222   vector < InttRawHit* > hits;
0223   auto raw_hit_num = node_inttrawhit_map_->get_nhits();
0224
0225   for (unsigned int i = 0; i < raw_hit_num; i++)
0226     {
0227       auto hit = node_inttrawhit_map_->get_hit(i);
0228
0229       hits.push_back( hit );
0230     }
0231
0232   return hits;
0233 }
0234
0235 vector < InttRawHit* > InttStreamingTiming::GetRawHitsWithoutClone()
0236 {
0237
0238   auto hits = this->GetRawHits();
0239   vector < InttRawHit* > rtn;
0240
0241   for( int i=0; i<hits.size(); i++ )
0242     {
0243
0244       bool is_same_found = false;
0245       for( int j=i+1; j<hits.size(); j++ )
0246         {
0247
0248           is_same_found = is_same_found || this->IsSame( hits[i], hits[j] );
0249           if( is_same_found == true )
0250             {
0251               break;
0252             }
0253
0254         } // end of for( hit_j )
0255
0256       // if this hit is still unique (same hit is not found), keep it
0257       if( is_same_found == false )
0258     {
0259       rtn.push_back( hits[i] );
0260     }
0261
0262     } // end of for( hit_i )
0263
0264   return rtn;
0265 }
0266
0267 vector < pair < InttNameSpace::Online_s,  unsigned int > >
0268 InttStreamingTiming::GetHits( TrkrHitSetContainer::ConstRange hitsets )
0269 {
0270   /*!
0271     @param ConstRange hitsets This is ConstRange = std::pair< ConstIterator, ConstIterator >
0272   */
0273
0274   vector < pair < InttNameSpace::Online_s, unsigned int > > hits_information;
0275   auto current = hitsets.first;
0276   while( current != hitsets.second )
0277     {
0278       auto hitset = (*current).second;
0279
0280       auto layer       = TrkrDefs::getLayer( (*current).first );
0281       // layer: 0, 1, 2   : MVTX
0282       //        3, 4, 5, 6: INTT
0283       //        6<        : TPC
0284
0285       // Take hits from only INTT
0286       if( !(3<= layer && layer <= 6) )
0287     continue;
0288
0289       auto detector_id = TrkrDefs::getTrkrId( (*current).first );
0290       auto phi         = InttDefs::getLadderPhiId( (*current).first ); //TrkrDefs::getPhiElement( (*current).first );
0291       auto z           = InttDefs::getLadderZId( (*current).first ); // TrkrDefs::getZElement( (*current).first );
0292       // cout << "-----> "
0293       //       << int(detector_id) << "\t" << int(layer) << "\t"
0294       //       << int(phi) << "\t" << int(z) ;
0295       //    //<< endl;
0296       // cout << " + " << (*current).first << "\t"
0297       //       << hitset->size()
0298       //       << endl;
0299
0300       auto hits = hitset->getHits(); // ConstRange = std::pair< ConstIterator, ConstIterator >
0301       auto current2 = hits.first; // ConstIterator = Map::const_iterator, here, Map = std::map < TrkrDefs::hitkey, TrkrHit* >, This Map is different from the Map for TrkrHitSetContainer!!!
0302       while( current2 != hits.second )
0303     {
0304       auto key = (*current2).first;
0305       auto hit = (*current2).second;
0306
0307       auto col = InttDefs::getCol( key );
0308       auto row = InttDefs::getRow( key );
0309
0310       InttNameSpace::Offline_s offline_hit;
0311       offline_hit.layer = layer;
0312       offline_hit.ladder_phi = phi;
0313       offline_hit.ladder_z = z;
0314       offline_hit.strip_x = row;
0315       offline_hit.strip_y = col;
0316
0317       InttNameSpace::Online_s online_hit = InttNameSpace::ToOnline( offline_hit );
0318
0319       auto hit_info = pair < InttNameSpace::Online_s, unsigned int >( online_hit, hit->getAdc() );
0320       hits_information.push_back( hit_info );
0321       // cout << "\t\tHit " << key << "\t"
0322       //      << col << "\t"
0323       //      << row << "\t"
0324       //      << hit->getAdc() << "\t|  "
0325       //      << online_hit.ldr << "\t"
0326       //      << online_hit.arm
0327       //      << endl;
0328
0329       current2++;
0330     }
0331
0332       current++;
0333     }
0334
0335   return hits_information;
0336 }
0337
0338 /////////////////////////////////////////////////////////////////////////
0339 // public
0340 /////////////////////////////////////////////////////////////////////////
0341
0342 void InttStreamingTiming::SetOutputDir( string dir )
0343 {
0344   if( dir != "" )
0345     {
0346       output_dir_ = dir;
0347     }
0348
0349   string run_num_str = string( 8 - to_string(run_num_).size(), '0' ) + to_string( run_num_ );
0350   output_root_ = output_dir_ + output_basename_ + run_num_str + ".root";
0351   output_txt_ = output_dir_ + output_basename_ + run_num_str + ".dat";
0352
0353 }
0354
0355 /////////////////////////////////////////////////////////////////////////
0356 // Fun4All stuff
0357 /////////////////////////////////////////////////////////////////////////
0358
0359 int InttStreamingTiming::Init(PHCompositeNode *topNode)
0360 {
0361   return Fun4AllReturnCodes::EVENT_OK;
0362 }
0363
0364 int InttStreamingTiming::InitRun(PHCompositeNode *topNode)
0365 {
0366
0367   // If an user sets the run number to recoConsts, it works. Can I assume it?
0368   recoConsts *rc = recoConsts::instance();
0369   run_num_ = rc->get_IntFlag( "RUNNUMBER" );
0370
0371   this->SetOutputDir();
0372   tf_output_ = new TFile( output_root_.c_str(), "RECREATE" );
0373
0374   // just a FPHX BCO distribution
0375   hist_fphx_bco_ = new TH1D( "fphx_bco", "FPHX BCO;INTT Local Clock [BCO];Counts",
0376                  120, 0, 120 );
0377
0378   hist_fphx_bco_raw_ = new TH1D( "fphx_bco_raw", "FPHX BCO (raw);INTT Local Clock [BCO];Counts",
0379                  120, 0, 120 );
0380
0381   hist_streaming_offset_ = new TH1D( "streaming_offset",
0382                      "BCO_{GL1} - BCO_{INTT GTM}; BCO_{GL1} - BCO_{INTT GTM} [BCO]; Counts",
0383                      256, -128, 128 );
0384
0385
0386   ofs_ = ofstream( output_txt_.c_str() ); // overwrite mode
0387   ofs_ << "#event" << "\t" // # means comment out
0388       << "GL1_BCO" << "\t"
0389       << "INTT_GTM_BCO" << "\t";
0390   for( int i=0; i<120; i++ )
0391     {
0392       ofs_ << "FPHX" << i << "\t";
0393     }
0394
0395   ofs_ << "\n"; // don't forget line break
0396
0397   return Fun4AllReturnCodes::EVENT_OK;
0398 }
0399
0400 int InttStreamingTiming::process_event(PHCompositeNode *topNode)
0401 {
0402
0403   // Get nodes to access data
0404   auto status = this->GetNodes(topNode);
0405
0406   // If a node is not available, skip this event
0407   if (status == Fun4AllReturnCodes::ABORTEVENT)
0408     return Fun4AllReturnCodes::ABORTEVENT;
0409
0410   // If no hit is found, skip it
0411   if( node_trkrhitset_map_->size() == 0 )
0412     {
0413       return Fun4AllReturnCodes::EVENT_OK;
0414     }
0415
0416   ////////////////////////////////////////////////////////////////////////
0417   // Get parameters of this event                                       //
0418   ////////////////////////////////////////////////////////////////////////
0419   event_counter_by_myself_++;
0420   auto bco_gl1 =  (gl1_->getBCO() & 0xFFFFFFFFFF);
0421   auto bco_intt = (node_inttrawhit_map_->get_hit( 0 )->get_bco());
0422
0423   auto trigger_bits = this->GetTriggerBits(); // 0, 1, 2, ...
0424   // if MBD N&S >= 1 triiger was not fired, skip it
0425   if( required_trigger_bit_ != -1 )
0426     if( trigger_bits[ required_trigger_bit_ ] == 0 )
0427       {
0428     return Fun4AllReturnCodes::EVENT_OK;
0429       }
0430
0431   event_counter_++;
0432
0433   ////////////////////////////////////////////////////////////////////////
0434   // Get TrkrHit on INTT                                                //
0435   ////////////////////////////////////////////////////////////////////////
0436   // get TrkrHitSet of INTT only, this is ConstRange = std::pair< ConstIterator, ConstIterator >
0437   auto hits = GetHits( node_trkrhitset_map_->getHitSets() );
0438   // FPHX BCO and event counter
0439   auto bco_event_counter_pair = this->GetBcoEventCounter(); // first: uint16_t (FPHX BCO), second: int (event counter)
0440
0441   ////////////////////////////////////////////////////////////////////////
0442   // loop over all hits                                                 //
0443   ////////////////////////////////////////////////////////////////////////
0444   int counts[ 128 ] = { 0 }; // BCO bin for this event
0445   for( int i=0; i<hits.size(); i++ )
0446     {
0447       auto adc = hits[i].second;
0448       int bco_fphx = bco_event_counter_pair[i].first;
0449
0450       //      counts[ bco_fphx ]++;
0451       counts[ bco_fphx ]++;
0452       hist_fphx_bco_raw_->Fill( bco_event_counter_pair[i].first ); // take all hits
0453       if( adc > 35 )
0454     {
0455       //counts[ bco_fphx ]++;
0456       hist_fphx_bco_->Fill( bco_event_counter_pair[i].first ); // take hits with ADC>1
0457
0458       auto bco_hit = bco_intt + bco_fphx;
0459       int diff = bco_hit - bco_gl1;
0460       hist_streaming_offset_->Fill( diff );
0461     } // end of if( adc > 35 )
0462
0463     } // end of for( int i=0; i<hits.size(); i++ )
0464
0465   auto raw_hits = this->GetRawHitsWithoutClone();
0466   for( auto& hit : raw_hits )
0467     {
0468       int bco_fphx = hit->get_FPHX_BCO();
0469       counts[ bco_fphx ]++;
0470     }
0471
0472   // write #hit for each BCO in this event to a text file
0473   if( event_counter_by_myself_ < event_max_ )
0474     {
0475       //      stringstream ss;
0476       ofs_ << event_counter_by_myself_ << "\t"
0477        << bco_gl1 << "\t"
0478        << bco_intt << "\t";
0479
0480       for( auto& val : counts )
0481     ofs_ << val << "\t";
0482
0483       ofs_ << "\n";
0484
0485       //ofs_ << ss.str();
0486     }
0487
0488   return Fun4AllReturnCodes::EVENT_OK;
0489 }
0490
0491 int InttStreamingTiming::ResetEvent(PHCompositeNode *topNode)
0492 {
0493
0494   return Fun4AllReturnCodes::EVENT_OK;
0495 }
0496
0497 int InttStreamingTiming::EndRun(const int runnumber)
0498 {
0499
0500   ofs_.close();
0501
0502   // Barrel correlation
0503   tf_output_->WriteTObject( hist_fphx_bco_, hist_fphx_bco_->GetName() );
0504   tf_output_->WriteTObject( hist_fphx_bco_raw_, hist_fphx_bco_raw_->GetName() );
0505   tf_output_->WriteTObject( hist_streaming_offset_, hist_streaming_offset_->GetName() );
0506
0507   // Close the ROOT file
0508   tf_output_->Close();
0509   return Fun4AllReturnCodes::EVENT_OK;
0510 }
0511
0512 int InttStreamingTiming::End(PHCompositeNode *topNode)
0513 {
0514
0515   return Fun4AllReturnCodes::EVENT_OK;
0516 }
0517
0518 int InttStreamingTiming::Reset(PHCompositeNode *topNode)
0519 {
0520
0521   return Fun4AllReturnCodes::EVENT_OK;
0522 }
0523
0524 void InttStreamingTiming::Print(const std::string &what) const
0525 {
0526
0527   int width = 100;
0528   cout << string( width, '-' ) << endl;
0529   cout << "InttStreamingTiming" << endl;
0530   cout << "  - Output (ROOT): " << output_root_ << endl;
0531   cout << string( width, '-' ) << endl;
0532 }