Initial commit

ascii2rootWF.cc

+/*
+Valerio, Insulab cosmic, ASACUSA bars
+Jan 2021
+
+
+=== 1st row ===
+4 floating (XY * (2 beam chambers))
+16 integers (CAEN V1730 digitizer: 8 wavef max + 8 max times)
+3 integers (iev1 + unix time + iev2)
+-------------------
+TOT 23 columns
+
+=== 2nd row ===
+512 ADC (positive)
+
+=== 3rd row ===
+512 ADC (positive)
+*/
+#include <string>
+#include <sstream>
+#include "Riostream.h"
+#include "TClassTable.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TGraph.h"
+
+TGraph * Filter(const TGraph * g, double alpha=0.2)
+{
+  TGraph * gbis = new TGraph();
+  double x,y,yprec;
+  g->GetPoint(0,x,y);
+  g->GetPoint(0,x,yprec);
+  gbis->SetPoint(0,x,alpha*y);
+  for(int i=1;i<g->GetN();i++) {
+    g->GetPoint(i,x,y);
+    gbis->SetPoint(i,x,alpha*y+(1-alpha)*yprec);
+    gbis->GetPoint(i,x,yprec);
+  }
+  return gbis;
+}
+TGraph * Deriv(const TGraph * g,int nsamples=3)
+{
+
+  TGraph * gbis = new TGraph();
+  double x,y,sum,average_delta;
+  for(int i=0;i<g->GetN();i++) {
+    sum=0;average_delta=0;
+    g->GetPoint(i,x,y);
+    if (i==0) {
+      average_delta=y;
+    }
+    
+    if (i>0 && i<nsamples) {
+      for (int j=0;j<i;j++) {sum += g->GetY()[i] - g->GetY()[i-j];}
+      average_delta=sum/i;
+    }
+    if (i>=nsamples) {
+      for (int j=0;j<nsamples;j++) {sum += g->GetY()[i] - g->GetY()[i-j];}
+      average_delta=sum/nsamples;
+    }
+    gbis->SetPoint(i,x,average_delta);
+  }
+  return gbis;
+}
+void SearchPeakstemp(vector <double> * vtime, vector <double> * vampli, TGraph * wave, double filter=0.9, int derivsamp=5)
+{
+  for(int i=0; i<wave->GetN(); i++) {
+    if ( wave->GetY()[i]>3000 ) {
+      vtime->push_back(i);
+      vampli->push_back(i);
+    }
+  }
+  return;
+}
+
+void SearchPeaks(vector <double> * vtime, vector <double> * vampli, vector <double> * vdeltat, TGraph * wave, double filter=0.9, int derivsamp=5)
+{
+  // graph deriv filt
+  TGraph * gf= new TGraph();
+  gf  = Filter(wave,filter);
+  
+  TGraph * gdf= new TGraph();
+  gdf = Deriv(gf,derivsamp);
+
+  float thr1 = filter * derivsamp*5;
+  float thr2 = 0;
+  float sum_over_threshold=0;
+  int isABOVE=0;
+  int peak_count=0;
+  for(int i=0; i<gdf->GetN(); i++) {
+    
+    // transition ABOVE thr1
+    if (i>0 && gdf->GetY()[i] > thr1 && gdf->GetY()[i-1]<=thr1 && !isABOVE) {
+      sum_over_threshold+=gdf->GetY()[i];
+      isABOVE=1;
+      // store transition time
+      vtime->push_back(i);
+      if (vtime->size()==1) {
+	vdeltat->push_back(0);
+      }
+      else {
+	int last = vtime->size();
+	vdeltat->push_back(vtime->at(last-1)-vtime->at(last-2));
+      }
+    }
+    // sum if it stays above thr
+    if (isABOVE) sum_over_threshold += gdf->GetY()[i];
+    // transition BELOW thr2
+    if (isABOVE && gdf->GetY()[i]<=thr2) {
+      vampli->push_back(sum_over_threshold);
+      sum_over_threshold=0;
+      isABOVE=0;
+      peak_count++;
+
+    }
+
+    if (isABOVE && i==gdf->GetN()-1) vampli->push_back(sum_over_threshold);
+
+    //if (peak_count>0) cout << "i "<<i<<" y deriv = "<<wave->GetY()[i]<<","<<gdf->GetY()[i]<<" => isABOVE "<<isABOVE<<" size T,A "<<vtime.size()<<","<<vampli.size()<<" (ampli "<<vampli[peak_count-1]<<")"<<endl;
+  }
+
+  
+  delete gf;
+  delete gdf;
+  return;
+}
+
+
+
+
+
+
+
+void ascii2rootWF(const int ndigi=1, const int nsili=4, string inFileName="dummy", string outFileName="out") {
+
+    
+//  string outFileName = "root_"+removeExtension(inFileName)+".root";
+  std::cout<<"---> Converting file "<<inFileName<<std::endl;
+  std::cout<<"--->       into file "<<outFileName<<std::endl;
+  ifstream myfile;
+  string line;
+  Int_t lineNumber=0;
+  myfile.open(inFileName.c_str()); 
+
+  // OUTPUT TREE DEFINITION  
+  TFile* outFile = new TFile(outFileName.c_str(),"RECREATE"); 
+  TTree* t = new TTree("t","Insulab tree from ascii file 2019");
+  //gROOT->cd();
+  // silicon trackers
+  Float_t hitPos[nsili];
+  t->Branch("hitPos",hitPos,Form("hitPos[%d]/F",nsili));
+  // digitizer
+  Float_t digiPh[8*ndigi]; 
+  Float_t digiTime[8*ndigi];
+  t->Branch("digiPh",digiPh,Form("digiPh[%d]/F",8*ndigi));
+  t->Branch("digiTime",digiTime,Form("digiTime[%d]/F",8*ndigi));
+  Int_t iev1, unixTime, iev2;
+  t->Branch("iev1",&iev1,"iev1/I");
+  t->Branch("unixTime",&unixTime,"unixTime/I");
+  t->Branch("iev2",&iev2,"iev2/I");
+
+  // 2 waveforms
+  Int_t wf2[512],wf3[512],x[512];
+  for (Int_t i=0;i<512;i++) x[i]=i;
+  t->Branch("wf2",wf2,Form("wf2[%d]/I",512));
+  t->Branch("wf3",wf3,Form("wf3[%d]/I",512));
+  // Float_t wf2Base=0, wf3Base=0;
+  // t->Branch("wf2Base",&wf2Base,"wf2Base/I");
+  // t->Branch("wf3Base",&wf3Base,"wf3Base/I");
+  // Int_t wf2Area[3],wf3Area[3],wf2Integral=0, wf3Integral=0;
+  // t->Branch("wf2Area",wf2Area,Form("wf2Area[%d]/I",3));
+  // t->Branch("wf3Area",wf3Area,Form("wf3Area[%d]/I",3));
+
+
+  // output vectors for peak search (time and amplitude)
+
+  vector <double> FoundPeak_time2, FoundPeak_time3;
+  vector <double> FoundPeak_ampl2, FoundPeak_ampl3;
+   vector <double> FoundPeak_deltat2, FoundPeak_deltat3;
+  t->Branch("FoundPeak_time2",&FoundPeak_time2);
+  t->Branch("FoundPeak_time3",&FoundPeak_time3);
+  t->Branch("FoundPeak_ampl2",&FoundPeak_ampl2);
+  t->Branch("FoundPeak_ampl3",&FoundPeak_ampl3);
+  t->Branch("FoundPeak_deltat2",&FoundPeak_deltat2);
+  t->Branch("FoundPeak_deltat3",&FoundPeak_deltat3);
+  
+  // std::fill(wf2, wf2+512, 1022);
+  // std::fill(wf3, wf3+512, 1033);
+  // end of OUTPUT TREE DEF
+
+  float filter=0.7;
+  int derivsamp=3;
+  if (myfile.is_open()) {
+    
+    //    while (!myfile.eof()) {
+    while(getline(myfile,line)){
+      std::stringstream ss(line);
+     
+      // 1st line
+      if (lineNumber%3==0) {
+     	for (Int_t i=0;i<nsili;i++) ss >> hitPos[i];
+	for (Int_t idigi=0;idigi<8*ndigi;idigi++) ss >> digiPh[idigi] ;
+	for (Int_t idigi=0;idigi<8*ndigi;idigi++) ss >> digiTime[idigi] ;
+	ss >> iev1 >> unixTime >> iev2;
+      }
+
+      if (lineNumber%3==1) {
+	//	wf2Base=0;wf2Integral=0;
+	for (Int_t i=0;i<512;i++) {
+	  ss >> wf2[i];
+	}
+	
+	// waveform analysys
+	FoundPeak_time2.clear();
+	FoundPeak_ampl2.clear();
+	FoundPeak_deltat2.clear();
+	TGraph * g2 = new TGraph(512,x,wf2);
+	SearchPeaks(&FoundPeak_time2,&FoundPeak_ampl2,&FoundPeak_deltat2,g2,filter,derivsamp);
+	delete g2;
+	
+      }
+      
+      if (lineNumber%3==2) {
+	//	wf3Base=0;wf3Integral=0;
+	for (Int_t i=0;i<512;i++) {
+	  ss >> wf3[i];
+	}
+	// waveform analysys
+	FoundPeak_time3.clear();
+	FoundPeak_ampl3.clear();
+	FoundPeak_deltat3.clear();
+	TGraph * g3 = new TGraph(512,x,wf3);
+	SearchPeaks(&FoundPeak_time3,&FoundPeak_ampl3,&FoundPeak_deltat3,g3,filter,derivsamp);
+	delete g3;
+
+	
+
+	// before going to the next line group, save it!
+	if (abs(hitPos[0]+hitPos[1]+hitPos[2]+hitPos[3])<100) {
+	  t->Fill();
+	}
+	
+      }
+      
+    
+      //std::cout <<  lineNumber<< std::endl;
+      lineNumber++;
+    }
+    myfile.close();
+
+  }
+
+  else {
+    std::cout<<"Can't open file!!"<<std::endl;
+  }
+  
+  outFile->cd();
+  t->Write();
+  
+  outFile->Close();
+    
+  return;
+}