/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 <TGraphErrors.h>
#include <TF1.h>
#include <TH1F.h>
#include <TCanvas.h>
#include <sstream>
#include "Log.hh"
#include "MTEvent.hh"
#include "MTRun.hh"
#include "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 - ASU NUMBER (13-14...)**************************************************************/"<<endl;
      cout<<"/* 4 - TRIGGER RATE (Hz)******************************************************************/"<<endl;
      cout<<"/* 5 - TARGET NOISE RATE (Hz)*************************************************************/"<<endl;
      cout<<"/* 6 - ACTIONS: R (readonly) - C (correctoffset) - I (increaseoffset) ********************/"<<endl;
      cout<<"/*****************************************************************************************/"<<endl;
    }

  else
    {

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




      if (!(asu>="13" && asu<="37"))
        {
          cout<<endl;
          cout<<"ASU NUMBER SHOULD BE BETWEEN 13 AND 37"<<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);


          
          //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);
                      
                      int hardid = channel->GetHardId();
                      int chipid = channel->GetChipId();
                
                      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);
                }
            }


          
          //DECLARE FILES TO WRITE TO AND READ FROM

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

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

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


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

          for (int i=0;i<144;i++)
            {
              for (int j=0;j<64;j++)
                {
                  chan[i][j] = 0;
                  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;chan[i-1][j] = 1;}

          //Read which channel have their offset already set
          while(in12>>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++)
                { 
                  if (chan[i][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 != "R")
                        {
                          //if rate > target, lower offset and define it as set
                          if (rate > target_rate)
                            {
                              offset[i][j]--;
                              set[i][j] = 1;
                              if (action == "C"){cout<<i+1<<" "<<j<<" "<<offset[i][j]+1<<" -> "<<offset[i][j]<<endl;}
                            }
                          
                          /*if rate <= target, increase it by 1 if not already set*/
                          else
                            {
                              if (!set[i][j])
                                {
                                  if (action == "I")
                                    {
                                      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 != "R")
            {

              TString outfile11 = infile11;
              TString outfile12 = infile12;

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

              for (int i=0;i<144;i++)
                {
                  for (int j=0;j<64;j++)
                    {
                      if (chan[i][j])
                        {
                          {
                            out11<<i+1<<" "<<j<<" "<<offset[i][j]<<endl;
                            if (set[i][j])
                              {
                                out12<<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);

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

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

        }


    }
}

---