timing/macro/Analysis.cc

0001 #define Analysis_cc
0002 #include "Analysis.hh"
0003
0004 Analysis::Analysis( string data_no_suffix, bool is_preliminary ) :
0005   is_preliminary_( is_preliminary )
0006 {
0007   this->Init( data_no_suffix );
0008 }
0009
0010 void Analysis::Init( string data_no_suffix )
0011 {
0012   // init paths
0013   root_path_ = data_no_suffix + ".root";
0014   text_path_ = data_no_suffix + ".dat";
0015
0016   if( is_preliminary_ == true )
0017     output_ += "_preliminary.pdf";
0018   else
0019     output_ += "_internal.pdf";
0020
0021   // Init data an outptut ROOT file
0022   this->InitRoot();
0023
0024   // Init data from a text file
0025   this->InitText();
0026
0027   // normalize FPHX BCO dist
0028   hist_fphx_bco_->Scale( 1.0 / event_num_ );
0029
0030   if( is_logy_ == true )
0031     {
0032       ymin_ = 0.1;
0033       ymax_ = 5e5;
0034     }
0035   else
0036     {
0037       ymin_ = 0.0;
0038       ymax_ = 600;
0039     }
0040
0041   return;
0042 }
0043
0044 void Analysis::InitRoot()
0045 {
0046   /*!
0047     @brief Init for data from a ROOT file
0048     @details A histogram of FPHX BCO distribution of all analyzed events
0049     offset of INTT streaming readout, which is normally 23
0050     are obtained.
0051    */
0052   TFile* tf = new TFile( root_path_.c_str(), "READ" );
0053
0054   string title = ";Beam clock [BCO];Number of INTT hits";
0055
0056   hist_fphx_bco_raw_ = (TH1D*)tf->Get( "fphx_bco_raw" );
0057   hist_fphx_bco_raw_->SetTitle( title.c_str() );
0058
0059   hist_fphx_bco_ = (TH1D*)tf->Get( "fphx_bco" );
0060   this->ConfigureHist( hist_fphx_bco_ );
0061   //hist_fphx_bco_->SetLineColor( kGray + 1 );
0062   hist_fphx_bco_->SetLineColor( kBlue + 1 );
0063   //hist_fphx_bco_->SetLineStyle( 2 );
0064
0065   TH1D* hist_offset = (TH1D*)tf->Get( "streaming_offset" );
0066   int the_position = 0;
0067   int largest = 0;
0068   for( int i=1; i<hist_offset->GetNbinsX()+1; i++ )
0069     {
0070       int content = hist_offset->GetBinContent( i );
0071       if( largest < content )
0072     {
0073       largest = content;
0074       the_position = i;
0075     }
0076
0077     }
0078
0079   streaming_offset_ = hist_offset->GetBinLowEdge( the_position );
0080 }
0081
0082 void Analysis::InitText()
0083 {
0084   /*!
0085     @brief Init data from a text file
0086     @details BCO of GL1 and INTT and #hit for each FPHX BCO are taken from the text file event by event.
0087    */
0088   ifstream ifs( text_path_.c_str() );
0089   if( ifs.fail() )
0090     {
0091       cerr << text_path_ << " is not found" << endl;
0092       return;
0093     }
0094
0095   string line;
0096   while( getline( ifs, line ) )
0097     {
0098       if( line[0] == '#' )
0099     continue;
0100
0101       Event* ev = new Event( line, streaming_offset_ );
0102       if( ev->GetInttGtmBco() == 396957258755 )
0103     {
0104       //ev->Print();
0105       // auto hist = ev->GetHist();
0106       // for( int i=1; i<hist->GetNbinsX()+1; i++ )
0107         // cout << setw(3) << i << " "
0108         //   << hist->GetBinContent( i ) << " "
0109         //   << ev->GetGl1Timing() << " "
0110         //   << ev->GetGl1Bin()
0111         //   << endl;
0112     }
0113       events_.push_back( ev );
0114     }
0115
0116   ifs.close();
0117
0118   event_num_ = events_[ events_.size()-1 ]->GetEventNum();
0119   return;
0120 }
0121
0122 ////////////////////////////////////////////////////////
0123 // private functions
0124 ////////////////////////////////////////////////////////
0125
0126 string Analysis::GetDate()
0127 {
0128   return "9/13/2024";
0129
0130   TDatime dt;
0131   int year  = dt.GetYear();
0132   int month = dt.GetMonth();
0133   int day   = dt.GetDay();
0134
0135   // format: mm/dd/yyyy
0136   std::stringstream ss;
0137   ss << month << "/" << day << "/" << year;
0138
0139   return ss.str();
0140 }
0141
0142 void Analysis::DrawWords( Event* event )
0143 {
0144
0145   ///////////////////////////////////////////////////////////////////////////////
0146   // Writting words in the canvas                                              //
0147   ///////////////////////////////////////////////////////////////////////////////
0148   TLatex* tex = new TLatex();
0149
0150   double line_height = 0.05;
0151   double first_margin = 0.005;
0152   double pos_y = 1.0 - top_margin_ + first_margin;// - line_height;
0153
0154   // Date
0155   tex->DrawLatexNDC( 0.7, pos_y,
0156              string("#it{" + GetDate() + "}").c_str() );
0157
0158   // sPHENIX Internal or sPHENIX Prelimnary
0159   pos_y -= line_height - first_margin + 0.025;
0160   double pos_x = 0.2;
0161   if( is_preliminary_ == false )
0162     {
0163       tex->DrawLatexNDC( pos_x, pos_y, "#it{#bf{sPHENIX}} Internal" );
0164     }
0165   else
0166     {
0167       tex->DrawLatexNDC( pos_x, pos_y, "#it{#bf{sPHENIX}} Preliminary" );
0168     }
0169
0170   // p+p 200 GeV
0171   pos_y -= line_height;
0172   tex->DrawLatexNDC( pos_x, pos_y, "Run-24 #it{p+p} 200 GeV, Run 50889" );
0173
0174   pos_y -= line_height;
0175   tex->DrawLatexNDC( pos_x, pos_y, "INTT Streaming readout " );
0176
0177   pos_y -= line_height;
0178   stringstream ss;
0179   ss << "BCO " << event->GetInttGtmBco();
0180   tex->DrawLatexNDC( pos_x, pos_y, ss.str().c_str() );
0181
0182   // pos_y -= line_height;
0183   // if( event->IsGl1Correspondings() )
0184   //   tex->DrawLatexNDC( pos_x, pos_y, "GL1 match!" );
0185   // else
0186   //   tex->DrawLatexNDC( pos_x, pos_y, "GL1 not match" );
0187
0188 }
0189
0190 TH1D* Analysis::GetShiftedHist( TH1D* hist_arg )
0191 {
0192   TH1D* hist = (TH1D*)hist_arg->Clone();
0193
0194   for( int i=1; i<hist->GetNbinsX()+1; i++ )
0195     {
0196       auto content = hist_arg->GetBinContent( i + shift_value_ );
0197       // cout << hist->GetNbinsX() << "\t"
0198       //       << setw(3) << i << "\t"
0199       //       << setw(3) << i + shift_value_ << "\t"
0200       //       <<  (hist->GetNbinsX()+1 - shift_value_ ) << "\t"
0201       //       << i + shift_value_ - hist->GetNbinsX() << "\t"
0202       //       << (
0203       //           (i + shift_value_ ) > hist->GetNbinsX() ?
0204       //           i + shift_value_ - hist->GetNbinsX() :
0205       //           i + shift_value_
0206       //        )
0207       //       << endl;
0208
0209       if( (i + shift_value_ ) > hist->GetNbinsX() )
0210     content = hist_arg->GetBinContent( i + shift_value_ - hist->GetNbinsX() );
0211
0212       hist->SetBinContent( i, content );
0213     }
0214
0215   return hist;
0216 }
0217
0218 void Analysis::DrawGl1Timing( Event* event )
0219 {
0220   double timing = event->GetGl1Timing() + 0.5; // absolute value but not bin index
0221
0222   TLine* line = new TLine();
0223   line->SetLineColor( kAzure + 1 );
0224   line->SetLineWidth( 2 );
0225   line->SetLineStyle( 7 );
0226   line->DrawLine( timing, ymin_, timing, ymax_ );
0227 }
0228
0229 ////////////////////////////////////////////////////////
0230 // public functions
0231 ////////////////////////////////////////////////////////
0232 void Analysis::Draw()
0233 {
0234
0235   SetsPhenixStyle();
0236
0237   TCanvas* c = new TCanvas( output_.c_str(), "title", 800, 800 );
0238   // c->SetTopMargin( top_margin_ );
0239   // c->SetRightMargin( right_margin_ );
0240   c->Print( ((string)c->GetName() + "[").c_str() );
0241
0242   ///////////////////////////////////////////////////////////////////////////////
0243   // Drawing hists                                                             //
0244   ///////////////////////////////////////////////////////////////////////////////
0245   int counter = 0;
0246   for( auto& event : events_ )
0247     {
0248       counter++;
0249
0250       // bool condition = event->GetBinMoreThan( 50 ) > 8;
0251       // bool condition = event->GetBinMoreThan( 10 ) > 12;
0252       // bool condition = event->GetBinMoreThan( 0 ) > 30; // standard cut
0253
0254       // dedicated cut for better plots
0255       // bool condition = event->GetBinMoreThan( 0 ) > 30 // for low background
0256       //    //|| event->GetBinMoreThan( 50 ) < 3  // for some high peaks
0257       //    || (event->IsGl1Correspondings() == false)  // for GL1 maching event
0258       //    || (event->IsHitInAbortGap() == true )
0259       //    ;
0260
0261       int ignore_smaller_than = 7;
0262       bool condition = true
0263     // && event->GetBinMoreThan( 10 ) < 30 // for low background
0264     // || event->GetBinMoreThan( 50 ) < 3  // for some high peaks
0265     && ( event->IsGl1Correspondings() == true )     // for GL1 maching event
0266     //&& ( event->IsHitInAbortGap() == false )        // no hit in abort gap
0267     && ( event->GetHitNumMachingGl1() > 30 )        // high GL1 peak
0268     && ( event->GetFatPeakNum( ignore_smaller_than ) < 4 )            // few fat peaks
0269     && ( event->GetHitNumMaxEvent() / event->GetHitNumMachingGl1() < 3 ) // suppression of the highest peak
0270     ;
0271
0272       condition = true
0273     // && event->GetBinMoreThan( 10 ) < 30 // for low background
0274     // || event->GetBinMoreThan( 50 ) < 3  // for some high peaks
0275     && ( event->IsGl1Correspondings() == true )     // for GL1 maching event
0276     //&& ( event->IsHitInAbortGap() == false )        // no hit in abort gap
0277     && ( event->GetHitNumMachingGl1() > 30 )        // high GL1 peak
0278     && ( event->GetFatPeakNum( ignore_smaller_than ) > 10  )            // few fat peaks
0279     //  && ( event->GetHitNumMaxEvent() / event->GetHitNumMachingGl1() < 3 ) // suppression of the highest peak
0280     ;
0281
0282     if( event->GetInttGtmBco() == 396912729035 )
0283     //if( event->GetInttGtmBco() == 396903901595 )
0284     {
0285       cout << "The event!!!" << endl;
0286       condition = true;
0287     }
0288       else
0289     continue;
0290
0291
0292       if( !condition )
0293         continue;
0294
0295       //cout << event->GetFatPeakNum( ignore_smaller_than ) << endl;
0296       // if( event->GetInttGtmBco() != 396957258755 )
0297       //    continue;
0298
0299       // cout << event->GetInttGtmBco() << "\t"
0300       //       << event->GetHitNumMachingGl1() << "\t"
0301       //       << event->GetHitNumMaxEvent() << "\t"
0302       //       << event->IsGl1Correspondings() << "\t"
0303       //       << event->GetHist()->GetBinContent( event->GetGl1Bin() )
0304       //       << endl;
0305
0306       hist_fphx_bco_->GetYaxis()->SetRangeUser( ymin_, ymax_ );
0307       hist_fphx_bco_->Draw( "HISTE" );
0308       gPad->SetLogy( is_logy_ );
0309
0310       event->GetHist()->Draw( "same" );
0311       event->GetHistGl1Match()->Draw( "same" );
0312
0313       TLegend* leg = new TLegend( 0.2, 0.55, 0.35, 0.7);
0314       leg->AddEntry( hist_fphx_bco_, "#INTT hits (100k strobes average)", "l" );
0315       leg->AddEntry( event->GetHist(), "#INTT hits (Single strobe)", "l" );
0316       leg->AddEntry( event->GetHistGl1Match(), "#INTT hits (Matching trigger timing)", "l" );
0317       leg->SetTextSize( 0.04 );
0318
0319       this->DrawWords( event );
0320       // this->DrawGl1Timing( event );
0321       leg->Draw();
0322       c->Print( c->GetName() );
0323
0324       if( counter > process_num_ )
0325     break;
0326
0327     }
0328
0329   counter = 0;
0330   int counter_gl1_match = 0;
0331   for( auto& event : events_ )
0332     {
0333       counter++;
0334       counter_gl1_match += ( event->IsGl1Correspondings() ? 1 : 0 );
0335     }
0336
0337   cout << "GL1 matching ratio: "
0338        << counter_gl1_match << " / " << counter << " = "
0339        << double(counter_gl1_match) / counter
0340        << endl;
0341
0342   c->Print( ((string)c->GetName() + "]").c_str() );
0343 }