/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/Microroc/Calibration/Calib_asu_microroc.cpp

Go to the documentation of this file.

#include <map>
#include <fstream>
#include "TKey.h"
#include <TROOT.h>
#include <TRint.h>
#include <TKey.h>
#include <TTree.h>
#include <TFile.h>
#include <TSystem.h>
#include <iostream>
#include <TF1.h>
#include <TH1F.h>
#include <TGraph.h>
#include <TCanvas.h>
#include <sstream>
#include "Log.hh"

#include "MicroException.hh"


#include "root/MTRun.hh"
#include "root/MTDetector.hh"
#include "root/MTChamber.hh"
#include "root/MTBoard.hh"
#include "root/MTChip.hh"
#include "root/MTEvent.hh"
#include "root/MTChannel.hh"
#include "root/MTMicrorocChip.hh"

using namespace std;


int main(int argc, char**argv){

  if (argc!=4){
    
    cout<<"/*****************************************************************/"<<endl;
    cout<<"/* ARGUMENTS TO THE PROGRAM **************************************/"<<endl;
    cout<<"/* 1 - DIRECTORY OF THE ROOT FILE PRODUCED BY THE RECONSTRUCTION */"<<endl;
    cout<<"/* 2 - NAME OF THE ROOT FILE *************************************/"<<endl;
    cout<<"/* 3 - THRESHOLD THAT IS MOVED DURING THE CALIBRATION ************/"<<endl;
    cout<<"/*    0 -> DAC_0 (LOW THRESHOLD) *********************************/"<<endl;
    cout<<"/*    1 -> DAC_1 (MEDIUM THRESHOLD) ******************************/"<<endl;
    cout<<"/*    2 -> DAC_2 (HIGH THRESHOLD) ********************************/"<<endl;
    cout<<"/*****************************************************************/"<<endl;}

  else{

  TString root_directory = argv[1];
  TString root_filename = argv[2];
  int dac = atoi(argv[3]);

  bool slab13_inversion = 0;
  float width_thr = 3;
  bool display = 0;
  bool print = 1;

  TString root_file = root_directory + root_filename;
  TString scurve_file = root_directory + "scurve_" + root_filename;

  cout<<endl;
  cout<<"Read file "<<root_file<<endl;
  cout<<"Wish to write "<<scurve_file<<endl;
  cout<<endl;

  int nbChannel = 0;
  float nbEvt = 0;
  int threshold = 0;
  int chipid = 0;
  int hardid = 0;
  int value = 0;
  int content = 0;
  int dac_min = 1025;
  int dac_max = 0;
  int scurve_max = 0;
  int deriv = 0;

  int N = 144;
 
  TH1I * hvalue = new TH1I("hvalue","",5,0,5);
  TH1I * hnpulse_threshold = new TH1I("hnpulse_threshold","",1024,0,1024);

  TH1F * hc_inflex = new TH1F("hc_inflex","Scurve derivative mean;channel number;mean (DAC)",9216,0,9216);
  TH1F * hc_width = new TH1F("hc_width","Scurve derivative RMS;channel number;RMS (DAC)",9216,0,9216);
  TH1F * hc_inflex_width = new TH1F("hc_inflex_width","Scurve derivative mean with RMS as error bar;channel number;mean (DAC)",9216,0,9216);

  TH1F * hinflex = new TH1F("hinflex","Scurve derivative mean;mean (DAC)",1000,0,1000);
  TH1F * hwidth  = new TH1F("hwidth","Scurve derivative RMS;RMS (DAC)",1000,0,100);

  hvalue->SetFillColor(2);
  hnpulse_threshold->SetFillColor(2);
  hc_inflex->SetLineColor(4);
  hc_width->SetLineColor(2);
  hinflex->SetFillColor(4);
  hwidth->SetFillColor(2);
  hinflex->SetLineColor(4);
  hwidth->SetLineColor(2);

  TH1I * hscurve[N][64];
  TH1I * hscurve_deriv[N][64];

  float inflex_tab[N][64];
  float inflex_err_tab[N][64];
  float width_tab[N][64];
  float width_err_tab[N][64];
  bool noisy[N][64];

  TString name,title;

  bool newchip[144];
  bool newchannel[144][64];

  for (int i=0;i<144;i++)
    {
      newchip[i] = 1;
      for (int j=0;j<64;j++)
        {
          inflex_tab[i][j] = 0;
          inflex_err_tab[i][j] = 0;

          width_tab[i][j] = 0;
          width_err_tab[i][j] = 0;

          noisy[i][j] = 0;
          newchannel[i][j] = 1;
        }
    }



  cout<<""<<endl;
  cout<<"******* READ CALIB DATA *******"<<endl;
  cout<<""<<endl;
  
  
  TFile f(root_file);
  TIter nextkey(f.GetListOfKeys());
  TKey *key;
  
  while (key = (TKey*)nextkey()) {

    TTree *tree = (TTree*)key->ReadObj();                
    MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");
                if ( run != NULL )
                {
                        const MTChip&  chip = run->GetOneChip();
                        try
                                {       const MTMicrorocChip &microChip = dynamic_cast<const MTMicrorocChip &> (chip);
                                        if (dac==0){threshold = microChip.GetThresholdDac_0();}
                                        if (dac==1){threshold = microChip.GetThresholdDac_1();}
                                        if (dac==2){threshold = microChip.GetThresholdDac_2();}
                                }

                        catch (...) {}
                }


    cout<<"DAC value = "<<threshold<<" \r "<<flush;

    /*Find DAC range*/
    if (threshold<dac_min){dac_min = threshold;}
    if (threshold>dac_max){dac_max = threshold;}

    MTEvent *evt =  new MTEvent();
    TBranch *branch= tree->GetBranch("MTEvent");
    branch->SetAddress(&evt);
    
    MTChannel* channel = NULL;
    
    nbEvt = tree->GetEntries();

    hnpulse_threshold->SetBinContent(threshold+1,nbEvt);


    /*Create and fill histogram for each ASU channel*/

    for ( int evtNum = 0; evtNum < nbEvt ; evtNum++){

      tree->GetEntry(evtNum);
      nbChannel = evt->GetNchannel();

      for(int i=0;i<nbChannel ;i++){    

        channel = (MTChannel*)evt->GetChannels()->UncheckedAt(i);
        chipid = channel->GetChipId();
        hardid = channel->GetHardId();

        /*2nd microroc m2 prototype: Invert slab1 and slab3*/

        if (slab13_inversion){
          if (chipid<=48){chipid += 96;}
          else{
            if (chipid>96){chipid -= 96;}}}

        if (newchip[chipid-1]){

          cout<<"Chip "<<chipid<<" found"<<endl;

          newchip[chipid-1] = 0;}

        if (newchannel[chipid-1][hardid]){

          name  = Form("hscurve_%i_%i",chipid,hardid);
          title = Form("SCurve - chip %i - channel %i",chipid,hardid);
          hscurve[chipid-1][hardid] = new TH1I(name,title,1024,0,1024);
          hscurve[chipid-1][hardid]->SetFillColor(2);

          name  = Form("hscurve_deriv_%i_%i",chipid,hardid);
          title = Form("Differential SCurve - chip %i - channel %i",chipid,hardid);
          hscurve_deriv[chipid-1][hardid] = new TH1I(name,title,1024,0.5,1024.5);
          hscurve_deriv[chipid-1][hardid]->SetFillColor(4);

          value = channel->GetDigitalValue();
          hvalue->Fill(value);
          if (value==(dac+1)){hscurve[chipid-1][hardid]->Fill(threshold);}

          newchannel[chipid-1][hardid] = 0;}

        else{

          value = channel->GetDigitalValue();
          hvalue->Fill(value);
          if (value==(dac+1)){hscurve[chipid-1][hardid]->Fill(threshold);}}}}


    tree->Delete();}



  /*Calculate derivative of Scurve*/

  cout<<endl;
  cout<<"Threshold DAC range = ["<<dac_min<<";"<<dac_max<<"] DAC units"<<endl;
  cout<<endl;

  for (int i=0;i<144;i++){
    for (int j=0;j<64;j++){
      if (newchannel[i][j]==0){

        for (int k=dac_min+1;k<dac_max;k++){

          deriv = hscurve[i][j]->GetBinContent(k)-hscurve[i][j]->GetBinContent(k+1);
          if (deriv>0){hscurve_deriv[i][j]->SetBinContent(k,deriv);}}}}}








  cout<<""<<endl;
  cout<<""<<endl;
  cout<<"******* PROCESS CALIB DATA *******"<<endl;
  cout<<""<<endl;


  float inflex = 0;
  float width = 0;
  float inflex_err = 0;
  float width_err = 0;


  for (int i=0;i<N;i++)
    {
      for (int j=0;j<64;j++)
        {
          /*If SCurve exist*/
          if (newchannel[i][j]==0)
            {
              inflex = hscurve_deriv[i][j]->GetMean();
              inflex_err = hscurve_deriv[i][j]->GetMeanError();
              width = hscurve_deriv[i][j]->GetRMS();
              width_err = hscurve_deriv[i][j]->GetRMSError();
  
              hc_inflex_width->SetBinContent(i*64+j+1,inflex);
              hc_inflex_width->SetBinError(i*64+j+1,width);

              hc_inflex->SetBinContent(i*64+j+1,inflex);
              hc_inflex->SetBinError(i*64+j+1,inflex_err);

              hc_width->SetBinContent(i*64+j+1,width);
              hc_width->SetBinError(i*64+j+1,width_err);

              hinflex->Fill(inflex);
              hwidth->Fill(width);

              if (width>width_thr || hscurve[i][j]->GetBinContent(dac_max)!=0)
                {
                  noisy[i][j] = 1;
                }
            }
        }
    }






  
  
 

  cout<<""<<endl;
  cout<<""<<endl;
  cout<<"******* WRITE CALIB DATA *******"<<endl;
  cout<<""<<endl;



  TString outfile = root_directory + "scurve_" + root_filename;
  for (int i=0;i<4;i++){outfile.Chop();}
  outfile += "txt";
  ofstream out(outfile.Data());

  TFile * tf = new TFile(scurve_file,"RECREATE");

  hvalue->Write();
  hnpulse_threshold->Write();

  for (int i=0;i<N;i++)
    {
    if (!newchip[i])
      {     
        for (int j=0;j<64;j++)
          {
            if (!newchannel[i][j])
              {
                hscurve[i][j]->Write();
                hscurve_deriv[i][j]->Write();

                inflex = hc_inflex->GetBinContent(i*64+j+1);
                inflex_err = hc_inflex->GetBinError(i*64+j+1);

                width = hc_width->GetBinContent(i*64+j+1);
                width_err = hc_width->GetBinError(i*64+j+1);

                cout<<i+1<<" "<<j<<" "<<inflex<<" "<<inflex_err<<" "<<width<<" "<<width_err<<" \r "<<flush;
                out<<i+1<<" "<<j<<" "<<inflex<<" "<<inflex_err<<" "<<width<<" "<<width_err<<endl;                   
              }
          }
      }
    }
  

  cout<<endl;
  cout<<endl;
  cout<<"Noisy channel with pedestal RMS > "<<width_thr<<" DAC"<<endl;

  TString outfile2 = root_directory + "noisy_" + root_filename;
  for (int i=0;i<4;i++){outfile2.Chop();}
  outfile2 += "txt";
  ofstream out2(outfile2.Data());

  for (int i=0;i<N;i++)
    {
    if (!newchip[i])
      {     
        for (int j=0;j<64;j++)
          {
            if (!newchannel[i][j])
              {
                if (noisy[i][j])
                  {
                    inflex = hc_inflex->GetBinContent(i*64+j+1);
                    inflex_err = hc_inflex->GetBinError(i*64+j+1);

                    width = hc_width->GetBinContent(i*64+j+1);
                    width_err = hc_width->GetBinError(i*64+j+1);

                    cout<<"  "<<i+1<<" "<<j<<" "<<inflex<<" "<<inflex_err<<" "<<width<<" "<<width_err<<endl;
                    out2<<i+1<<" "<<j<<" "<<inflex<<" "<<inflex_err<<" "<<width<<" "<<width_err<<endl;
                  }
              }
          }
      }
    }
  


  hc_inflex_width->Write();
  hc_inflex->Write();
  hc_width->Write();
  hinflex->Write();
  hwidth->Write();




  if (print)
    {

      cout<<""<<endl;
      cout<<""<<endl;
      cout<<"******* PRINT CALIB DATA *******"<<endl;
      cout<<""<<endl;

      int channel_min = 0;
      int channel_max = 0;

      float value_min = 0;
      float value_max = 0;

      for (int i=0;i<9216;i++)
        {
          if (hc_inflex->GetBinContent(i+1)!=0)
            {
              channel_min = i-100;
              i=9216;
            }
        }

      for (int i=9216;i>0;i--)
        {
          if (hc_inflex->GetBinContent(i)!=0)
            {
              channel_max = i-1+100;
              i=0;
            }
        }

      title = Form("Summary plots of calibration run : %s",root_filename.Data());
      TCanvas * c0 = new TCanvas("c0",title,10,10,1200,800);
      c0->Divide(2,2);

      c0->cd(1);
      value_min = hinflex->GetMean() - hinflex->GetRMS()*5;
      value_max = hinflex->GetMean() + hinflex->GetRMS()*5;
      hc_inflex->GetYaxis()->SetTitleOffset(1.5);
      hc_inflex->GetXaxis()->SetRangeUser(channel_min-10,channel_max+10);
      hc_inflex->GetYaxis()->SetRangeUser(value_min,value_max);
      hc_inflex->Draw();

      c0->cd(3);
      hinflex->GetXaxis()->SetRangeUser(value_min,value_max);
      hinflex->Draw();

      c0->cd(2);
      value_min = 0;
      value_max = hwidth->GetMean() + hwidth->GetRMS()*5;
      hc_width->GetYaxis()->SetTitleOffset(1.5);
      hc_width->GetXaxis()->SetRangeUser(channel_min-10,channel_max+10);
      hc_width->GetYaxis()->SetRangeUser(value_min,value_max);
      hc_width->Draw();
      
      c0->cd(4);
      hwidth->GetXaxis()->SetRangeUser(value_min,value_max);
      hwidth->Draw();
      
      TString print_filename = root_directory + "plots_" + root_filename;
      for (int i=0;i<4;i++){print_filename.Chop();}
      print_filename += "png";

      cout<<"Print file "<<print_filename<<endl;
      c0->Print(print_filename);
    }


  tf->Close();

  cout<<endl;
  cout<<endl;
  }




    return 0 ;}

---