ANASEN_analysis/script.C

// #include "/home/splitpole/FSUDAQ_Qt6/Aux/fsuReader.h"
// #include "/home/splitpole/FSUDAQ_Qt6/Aux/fsutsReader.h"
#include <TGraph.h>
#include <TFile.h>
#include <TTree.h>
#include <TROOT.h>
#include <TString.h>
#include <TMath.h>
#include <TRandom.h>

#include "mapping.h"

#include "Armory/ClassAnasen.h"

class PulserChecker {
public:
  PulserChecker(int sn) : SN(sn){
    t0 = 0;
    t1 = 0;
    dt = 0;
    tStart = -1;
    has1stData = false;
    n = 0;
    mean = 0;
    m2 = 0;
  }

  void addTime(unsigned long long time){
    t0 = t1;
    t1 = time;
    if( has1stData ) {
      dt = t1/1e6 - t0/1e6;
      int mod = 1;
      if( dt >= 2* mean && mean > 0 ){ 
        mod = std::round(dt/mean); 
      }
      if( mod > 1 ) return;
      double new_dt = dt/mod;
      
      n ++;
      double delta = new_dt - mean;
      mean += delta / n;

      double delta2 = new_dt - mean;
      m2 += delta * delta2;
      //if( mod > 1 ) printf("%6d | %llu, %llu (dt %.6f, m : %d) | %.6f\n", SN, t1, t0, dt, mod, new_dt);
    }else{
      tStart = time;
      has1stData = true; 
    }
    // printf("%6d | %13llu, %13llu | %.6f | %13llu\n", SN, t1, t0, dt, tStart);
  }

  unsigned long long getTime0() const {return tStart;}
  int                getSN()    const {return SN;}
  double             getDT()    const {return has1stData ? dt : TMath::QuietNaN() ;}
  double             getMean()  const {return mean;}
  double             getSTD()   const { return n < 2 ? 0.0 :sqrt(m2 / (n - 1)); }

  void Print() { printf("%6d (%2d) | %16llu| n : %6d, mean : %10.6f msec, sd : %10.6f msec\n", 
             SN, SN & 0x1F, tStart, n, getMean(), getSTD());} 

private:
  const int SN;
  unsigned long long tStart, t0, t1;
  double dt;
  bool has1stData;

  int n;
  double mean;
  double m2;

};

void script(TString fileName = "", int maxEvent = -1){

  /*
  //+++++++++++++++++++++++++++++++++++++++++++
  FSUReader * reader = new FSUReader("pulsertest_013_15528_QDC_16_000.fsu", 100);

  reader->ScanNumBlock(1, 1);
  reader->PrintHit(10);

  // for( int i = 0; i < 24 ; i++){
  //   printf("#################### agg : %d \n", i);

  //   reader->ReadNextBlock(false, 0, 1);

  //   reader->GetData()->PrintAllData();

  // }


  // TGraph * graph = new TGraph();
  // int count = 0;
  // unsigned long long tree0 = 0;
  // unsigned long numHit = reader->GetHitCount();
  // for( unsigned long i = 0; i < numHit; i++){
  //   //printf("-------------%6lu \n\033[A\r", i);
  //   unsigned long long t1 = reader->GetHit(i).timestamp;
  //   unsigned short ch = reader->GetHit(i).ch;
  //   unsigned short ee = reader->GetHit(i).energy;
  //   //printf("%6lu | %2u, %6u, %16llu\n", i, ch, ee, t1);
  //   count ++;
  //   graph->SetPoint(i, i, t1/1e9);
  // } 
  // graph->Draw("APL");

  // for( int i = 40; i < 60; i++){
  //   reader->GetHit(i).Print();
  // }
  */

  /*
  //+++++++++++++++++++++++++++++++++++++++++++

  TFile * file0 = new TFile("test_002_1000_noTrace.root"); 
  TFile * f1 = new TFile("test_002_1000.root"); 

  TTree * tree0 = (TTree*) file0->Get("tree");
  TTree * t1 = (TTree*) f1->Get("tree");

  unsigned int m0, m1;
  unsigned short ch0[100], ch1[100];
  unsigned long long ts0[100], ts1[100];


  tree0->SetBranchAddress("multi", &m0);
  tree0->SetBranchAddress("ch", ch0);
  tree0->SetBranchAddress("e_t", ts0);

  t1->SetBranchAddress("multi", &m1);
  t1->SetBranchAddress("ch", ch1);
  t1->SetBranchAddress("e_t", ts1);

  unsigned long long nEntries0 = tree0->GetEntries();
  unsigned long long nEntries1 = t1->GetEntries();

  for( unsigned long long ev = 0; ev < nEntries0; ev++){
    tree0->GetEntry(ev);
    t1->GetEntry(ev);


    if( ch0[0] != ch1[0] || ts0[0] != ts1[0] ) {
      printf("======================= %llu \n", ev);
      printf("ch : %u   %u\n", ch0[0], ch1[0]);
      printf("ts : %llu   %llu\n", ts0[0], ts1[0]);
    }
  }

  file0->Close();
  f1->Close();
  */
  /*
  //+++++++++++++++++++++++++++++++++++++++++++
  printf("######### file : %s \n", fileName.Data());
  TFile * file0 = new TFile(fileName); 

  TTree * tree0 = (TTree*) file0->Get("tree");

  const int MAX_MULTI = 1000;
  unsigned long long                evID = 0;
  unsigned int                     multi = 0;
  unsigned short           sn[MAX_MULTI] = {0};
  unsigned short           ch[MAX_MULTI] = {0};
  unsigned short            e[MAX_MULTI] = {0};
//  unsigned short           e2[MAX_MULTI] = {0};
  unsigned long long      e_t[MAX_MULTI] = {0};
//  unsigned short          e_f[MAX_MULTI] = {0};

  tree0->SetBranchAddress("evID", &evID);
  tree0->SetBranchAddress("multi", &multi);
  tree0->SetBranchAddress("sn", sn);
  tree0->SetBranchAddress("ch", ch);
  tree0->SetBranchAddress("e", e);
  tree0->SetBranchAddress("e_t", e_t);

  unsigned long long nEntries = tree0->GetEntries();

  const std::map<int, int> board = {
    {0, 17122},
    {1, 17123},
    {2, 22320},
    {3, 22130},
    {4, 22129},
    {5, 15529},
    {6, 15528},
    {7,   379},
    {8,   409},
    {9,   405},
  };
  const int nBd = board.size();

  PulserChecker * stat[9];
  for( int i = 0; i < nBd; i++){ stat[i] = new PulserChecker(board.at(i)); }

  printf("Number of Entries : %llu\n", nEntries);

  unsigned long long evEnd = nEntries;
  if( maxEvent >= 0 ) evEnd = maxEvent;

  for( unsigned long long ev = 0; ev < evEnd; ev++){
    tree0->GetEntry(ev);

    // printf("sn : %5d | %14llu \n", sn[0], e_t[0]);

    for( int i = 0; i< nBd; i++){
      for( unsigned int j = 0; j < multi; j++){
        if( sn[j] == stat[i]->getSN() ) {
          stat[i]->addTime(e_t[j]);
          printf(" sn : %5d | %16llu \n", sn[j], e_t[j]);
          break;
        }
      }
    }

    if( ev > 1){
      unsigned long long ts1 = e_t[0];
      tree0->GetEntry(ev-1);
      unsigned long long ts0 = e_t[0];

      if( ts1 <= ts0 )printf("--------------- ev : %llu\n", ev);
    }
  }
  file0->Close();

  printf("=========================\n");

  for( int i = 0; i< nBd; i++){ stat[i]->Print(); }

  unsigned long long time0 = stat[0]->getTime0(); 
  for( int i = 1; i< nBd; i++){ 
    if( stat[i]->getTime0() < time0 ) time0 = stat[i]->getTime0();
  }

  printf("=========================\n");
  printf(" S/N      time Offset ns \n");
  for( int i = 0; i< nBd; i++){ 
    printf(" %5d | %16llu \n", stat[i]->getSN(), stat[i]->getTime0() - time0);
  }
  */
  //+++++++++++++++++++++++++++++++++++++++++++
  // TCanvas * c1 = new TCanvas();
  // ANASEN * kaka = new ANASEN();
  // kaka->DrawAnasen();

  TCanvas * c2 = new TCanvas("c2", "ANASEN Simulation", 800, 800);
  c2->cd();


  ANASEN * haha = new ANASEN();

  haha->SetUncertainties(0, 0, 0, 0);

  PW * pw = haha->GetPW();
  SX3 * sx3 = haha->GetSX3();

  // haha->DrawAnasen(0, 1, 0, 1, -1, false);

  gRandom->SetSeed(0);

  int xRan = gRandom->Integer(20) - 10;
  int yRan = gRandom->Integer(20) - 10;
  int zRan = gRandom->Integer(20) - 10;
  int pRan = gRandom->Integer(60) -120; // phi deg

  int tRan = 0 ;
  do{
    tRan = gRandom->Integer(100) + 40 ; // theta deg
  }while( 75 < tRan && tRan < 105); 

  TVector3 pos (xRan, yRan, zRan);
  TVector3 dir (1, 0, 0);

  dir.SetPhi( pRan * TMath::DegToRad());
  dir.SetTheta( tRan * TMath::DegToRad());

  pos.Print();
  dir.Print();

  haha->DrawTrack(pos, dir, true);

  pw->FindWireID(pos, dir, true);
  sx3->FindSX3Pos(pos, dir, true);


  // haha->CalTrack(sx3.hitPos, wireID.first, wireID.second, true);

  //haha->DrawDeducedTrack(sx3.hitPos, wireID.first, wireID.second);


}