/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/Microroc/Calibration/Align_noise_rate_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 <TreeClass.hh>
#include <TGraphErrors.h>
#include <TF1.h>
#include <TH1F.h>
#include <TH1I.h>
#include <TH2I.h>
#include <TCanvas.h>
#include <sstream>
#include "Log.hh"

#include "root/CaloRun.hh"
#include "root/CaloEvent.hh"
#include "root/MTRun.hh"
#include "root/MTEvent.hh"
#include "root/MTChannel.hh"


#include "MicroException.hh"


using namespace std;


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

  if (argc!=7)
    {
      cout<<"/*****************************************************************************************/"<<endl;
      cout<<"/* ARGUMENTS TO THE PROGRAM **************************************************************/"<<endl;
      cout<<"/* 1 - DIRECTORY CONTAINING THE ROOTFILE *************************************************/"<<endl;
      cout<<"/* 2 - NAME OF MERGED ROOTFILE ***********************************************************/"<<endl;
      cout<<"/* 3 - M2 PROTOTYPE NUMBER (1 OR 2)*******************************************************/"<<endl;
      cout<<"/* 4 - TRIGGER RATE (Hz)******************************************************************/"<<endl;
      cout<<"/* 5 - TARGET NOISE RATE (Hz)*************************************************************/"<<endl;
      cout<<"/* 6 - ACTIONS: readonly - correctoffset - increaseoffset ********************************/"<<endl;
      cout<<"/*****************************************************************************************/"<<endl;
    }

  else
    {

      TString dir  = argv[1];
      TString rootfilename  = argv[2];
      TString m2  = atoi(argv[3]);
      float trigger_rate = atof(argv[4]);
      float target_rate = atof(argv[5]);
      TString action = argv[6];

      if (!(m2==1 || m2==2))
        {
          cout<<endl;
          cout<<"M2 NUMBER SHOULD BE 1 OR 2"<<endl;
          cout<<endl;
        }

      else
        {

          TString filename = dir + rootfilename;
          int nevent = 0;

          TH1F * hnhit = new TH1F("hnhit",";nhit",1000,0,10000);
          TH1F * hc_nhit = new TH1F("hc_nhit",";channel;number of hits",9216,0,9216);
          TH1F * hc_rate = new TH1F("hc_rate",";channel;rate (Hz)",9216,0,9216);
          TH1F * hrate = new TH1F("hrate",";rate (Hz)",100,0,1);
          TH1I * hoffset = new TH1I("hoffset","offset distribution;offset (DAC)",16,0,16);
          TH2I * hxy = new TH2I("hxy","hit distribution;y (pad);x (pad)",96,0,96,96,0,96);
          
          //READ ROOTFILE
      
          TFile *f = new TFile(filename);
          cout<<"reading File : "<<filename<<endl;

          TIter nextkey(f->GetListOfKeys());
          TKey *key;

          while (key = (TKey*)nextkey())
            {
              TTree *tree = (TTree*)key->ReadObj();                
              MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");

              nevent = tree->GetEntries();
              cout <<"nEvent="<<nevent<<endl;

              MTEvent *evt = new MTEvent();

              TBranch *branch= tree->GetBranch("MTEvent");
              branch->SetAddress(&evt);

              for (int i=0;i<nevent;i++)
                { 
                  tree->GetEntry(i);
                  int nchannel = evt->GetNchannel();
                  hnhit->Fill(nchannel);

                  MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(0);    

                  for (int j=0;j<nchannel;j++)
                    {
                      MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);

                      hxy->Fill(channel->GetY()*1e-4+48,channel->GetX()*1e-4+48,1);

                      int hardid = channel->GetHardId();
                      int chipid = channel->GetChipId();
                      int boardid = channel->GetBoardId();
                
                      if (0)
                        {
                          if (m2==1)
                            {
                              if (!(boardid==13 || boardid==14))
                                {
                                  if (boardid==15 || boardid==16){chipid += 48;}
                                  else{chipid += 96;}
                                }
                            }
                          
                          if (m2==2)
                            {
                              if (!(boardid==24 || boardid==25))
                                {
                                  if (boardid==22 || boardid==23){chipid += 48;}
                                  else{chipid += 96;}
                                }
                            }
                        }

                      hc_nhit->Fill((chipid-1)*64+hardid);
                    }
                }
            }


          //CALCULATE RATE FROM NUMBER OF HITS

          for (int i=0;i<9216;i++)
            {
              int entry = hc_nhit->GetBinContent(i+1);
              
              if (entry!=0)
                {
                  float rate = entry / 1. / nevent * trigger_rate ;
                  hc_rate->SetBinContent(i+1,rate);
                  hrate->Fill(rate);
                  //hrate->Fill(rate);
                }
            }


          
          //DECLARE FILES TO WRITE TO AND READ FROM

          TString slab1;
          TString slab2;
          TString slab3;

          if (m2 == 1)
            {
              slab1 = "ASU1314";
              slab2 = "ASU1516";
              slab3 = "ASU1718";
            }
          
          if (m2 == 2)
            {
              slab1 = "ASU2425";
              slab2 = "ASU2223";
              slab3 = "ASU2021";
            }


          TString infile11 = dir + slab1 + "_gains_noise_align.txt";
          TString infile12 = dir + slab1 + "_offset_set.txt";

          TString infile21 = dir + slab2 + "_gains_noise_align.txt";
          TString infile22 = dir + slab2 + "_offset_set.txt";

          TString infile31 = dir + slab3 + "_gains_noise_align.txt";
          TString infile32 = dir + slab3 + "_offset_set.txt";

          ifstream in11(infile11.Data());
          ifstream in12(infile12.Data());

          ifstream in21(infile21.Data());
          ifstream in22(infile22.Data());

          ifstream in31(infile31.Data());
          ifstream in32(infile32.Data());

          cout<<endl;
          cout<<"Read offset status (set/not-set) from files:"<<endl;
          cout<<"  ->  "<<infile12<<endl;
          cout<<"  ->  "<<infile22<<endl;
          cout<<"  ->  "<<infile32<<endl;
          cout<<endl;
          cout<<"Write offset map to files"<<endl;
          cout<<"  ->  "<<infile11<<endl;
          cout<<"  ->  "<<infile21<<endl;
          cout<<"  ->  "<<infile31<<endl;
          cout<<endl;


          //READ WHICH CHANNELS HAVE THEIR OFFSET ALREADY SET
          
          bool set[144][64];
          int offset[144][64];

          for (int i=0;i<144;i++)
            {
              for (int j=0;j<64;j++)
                {
                  set[i][j] = 0;
                  offset[i][j] = -1;
                }
            }

          int i,j,off;

          //Read offset maps
          while(in11>>i>>j>>off){offset[i-1][j] = off;}
          while(in21>>i>>j>>off){offset[i-1][j] = off;}   
          while(in31>>i>>j>>off){offset[i-1][j] = off;}

          //Read which channel have their offset already set
          while(in12>>i>>j){set[i-1][j] = 1;}
          while(in22>>i>>j){set[i-1][j] = 1;}
          while(in32>>i>>j){set[i-1][j] = 1;}


          //For non-set offsets
          //Determine if offset should be raised or lowered
          //by comparing measured rate to target rate

          int nthr = 0;

          //loop over m2 channels
          for (int i=0;i<144;i++)
            {
              for (int j=0;j<64;j++)
                {
 
                  off = offset[i][j];
                  float rate = hc_rate->GetBinContent(i*64+j+1);

                  if (rate > target_rate)
                    {
                      nthr++;
                      //cout<<i+1<<" "<<j<<" rate = "<<rate<<" Hz - offset:"<<offset[i][j]<<" -> "<<offset[i][j]-1<<endl;
                    }

                  if (action != "readonly")
                    {
                      //if rate > target, lower offset and define it as set
                      if (rate > target_rate)
                        {
                          offset[i][j]--;
                          set[i][j] = 1;
                          if (action == "correctoffset"){cout<<i+1<<" "<<j<<" "<<offset[i][j]+1<<" -> "<<offset[i][j]<<endl;}
                        }

                      if (action == "increaseoffset")
                        {
                          if (!set[i][j])
                            {
                              offset[i][j]++;
                              cout<<i+1<<" "<<j<<" "<<offset[i][j]-1<<" -> "<<offset[i][j]<<endl;
                            }
                        }
                    }
                }
            }



          //Re-write offset map and channel set files

          if (action == "correctoffset" || action == "increaseoffset")
            {

              TString outfile11 = infile11;
              TString outfile12 = infile12;

              TString outfile21 = infile21;
              TString outfile22 = infile22;

              TString outfile31 = infile31;
              TString outfile32 = infile32;

              ofstream out11(outfile11.Data());
              ofstream out12(outfile12.Data());

              ofstream out21(outfile21.Data());
              ofstream out22(outfile22.Data());

              ofstream out31(outfile31.Data());
              ofstream out32(outfile32.Data());


              for (int i=0;i<144;i++)
                {
                  for (int j=0;j<64;j++)
                    {
                      //slab1      
                      if (i<48)
                        {
                          out11<<i+1<<" "<<j<<" "<<offset[i][j]<<endl;
                          hoffset->Fill(offset[i][j]);
                          if (set[i][j])
                            {
                              out12<<i+1<<" "<<j<<endl;
                            }
                        }
                      else
                        {
                          //slab2
                          if (i<96)
                            {
                              out21<<i+1<<" "<<j<<" "<<offset[i][j]<<endl;
                              hoffset->Fill(offset[i][j]);
                              if (set[i][j])
                                {
                                  out22<<i+1<<" "<<j<<endl;
                                }
                            }
                          //slab3
                          else
                            {
                              out31<<i+1<<" "<<j<<" "<<offset[i][j]<<endl;
                              hoffset->Fill(offset[i][j]);
                              if (set[i][j])
                                {
                                  out32<<i+1<<" "<<j<<endl;
                                }
                            }
                        }
                    }
                }
            }

          cout<<endl;
          cout<<"Number of channel with noise rate large than "<<target_rate<<" Hz  = "<<nthr<<endl;
          cout<<"Average number of hits per readout = "<<hnhit->GetMean()<<endl;
          cout<<endl;


          
          TCanvas * c0 = new TCanvas("c0","",10,10,1000,500);

          cout<<endl;
          cout<<"Write histos to files: "<<endl;

          hnhit->Draw();
          TString picfilename = dir + "hnhit.png";
          cout<<picfilename<<endl;
          c0->Print(picfilename);

          hoffset->Draw();
          picfilename = dir + "hoffset.png";
          cout<<picfilename<<endl;
          c0->Print(picfilename);

          hxy->Draw("zcol");
          picfilename = dir + "hxy.png";
          cout<<picfilename<<endl;
          c0->Print(picfilename);

          c0->SetLogy();
          hc_rate->Draw();
          picfilename = dir + "hc_rate.png";
          cout<<picfilename<<endl;
          c0->Print(picfilename);

          hrate->Draw();
          picfilename = dir + "hrate.png";
          cout<<picfilename<<endl;
          c0->Print(picfilename);

        }


    }
}

---