ROOT36: 波形分析代码2

#include "iostream"

using namespace std;

void ana2(){

  // 1. read the WFM.root file, and get data

  TFile *ifile = new TFile( "../root35/WFM.root", "read" );

  TTree *itree = (TTree*)ifile->Get( "WFM" );

  double timestamp =0.0;

  vector<double> *time = 0;

  vector<double> *channel1 = 0;

  vector<double> *channel2 = 0;

  itree->SetBranchAddress( "TimeStamp", &timestamp );

  itree->SetBranchAddress( "Time", &time );

  itree->SetBranchAddress( "Channel1", &channel1 );

  itree->SetBranchAddress( "Channel2", &channel2 );

  

  // 2. analyze the waveform and get some charateristic values, ex. timestamp, peakenergy, and flag

  // energy = peakheight *511 / 0.31 [keV]

  // time = 5*timeindex [ns]

  double TimeStamp;

  double PeakHeight1;

  double PeakEnergy1;

  double PeakHeight2;

  double PeakEnergy2;

  double PeakTime1;

  double PeakTime2;

  bool Flag_Coin;

  const int N = itree->GetEntries();

  TFile *ofile = new TFile( "WFM_Analysis.root", "recreate" );

  TTree *otree = new TTree( "WFM", "analyzed WFM" );

  otree->Branch( "TimeStamp", &TimeStamp, "TimeStamp/D" );

  otree->Branch( "PeakEnergy1", &PeakEnergy1, "PeakEnergy1/D"  );

  otree->Branch( "PeakEnergy2", &PeakEnergy2, "PeakEnergy2/D"  );

  otree->Branch( "PeakTime1", &PeakTime1, "PeakTime1/D"  );

  otree->Branch( "PeakTime2", &PeakTime2, "PeakTime2/D"  );

  otree->Branch( "Flag_Coin", &Flag_Coin, "Flag_Coin/O" );

  //3 . save the analyzed data into a new file WFM_analysis.root

  for( int j=0; j<N; j++ ){

    // j is the number-th of event

    cout << "Analyzing " << j << "-th event" <<"\r" << flush;

    

    itree->GetEntry(j);

    TimeStamp = timestamp;

    auto pk1 = std::min_element( channel1->begin(), channel1->end() );

    PeakHeight1 = *pk1;

    PeakEnergy1 = -PeakHeight1*(511/0.31);

    auto PeakIndex1 = std::distance( channel1->begin(), pk1 );

    PeakTime1 = 5*PeakIndex1;

    auto pk2 = std::min_element( channel2->begin(), channel2->end() );

    PeakHeight2 = *pk2;

    PeakEnergy2 = -PeakHeight2*(511/0.31);

    auto PeakIndex2 = std::distance( channel2->begin(), pk2 );

    PeakTime2 = 5*PeakIndex2;

    if( (450<PeakEnergy2&&PeakEnergy2<560)&&(400<PeakTime2&&PeakTime2<700) ){

      Flag_Coin = 1;

    }else{

      Flag_Coin = 0;

    }

    

    otree->Fill();

  }

  cout << endl;

  otree->Write();

  ofile->Close();

}