/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/root/square_event_2D.cpp

Go to the documentation of this file.

//############################################
//
//--- Identify and count the square events ---
//
//############################################


#include <map>
#include <iostream>
#include "TH1I.h"
#include "TH2I.h"
#include "TCanvas.h"
#include "TFile.h"
#include <TROOT.h>
#include <TRint.h>
#include <TKey.h>
#include "TTree.h"
#include <TApplication.h>

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


using namespace std;



void add(map<unsigned int,map<unsigned int ,int> >&channels_map, unsigned int asu, unsigned int hit_time){

  if ( channels_map.find(asu) == channels_map.end()) { // if this asu doesn't exist yet
    
    map<unsigned int, int> nbHitMap;
    nbHitMap[hit_time]=1;
    channels_map.insert(make_pair(asu,nbHitMap));
  }
  else{ // if this asu already exist 

    map<unsigned int, int> &nbHitMap =  channels_map.find(asu)->second;
    if (nbHitMap.find(hit_time) == nbHitMap.end()){ //if this hit_time doesn't exist yet

      nbHitMap.insert(make_pair(hit_time,1)); 
    }
    else{ //if this abs time already exist

      nbHitMap[hit_time] = nbHitMap[hit_time]++;
    }
  }  
}



int main(void){  
  
  //-----------------
  //User's parameters
  
  //do you want to plot (0 = no plot, 1 = plots)
  bool plot = 1;
  
  //Do you want some info ? (0 = not at all, 1 = yes a little, 2 = yes a lot, 3 = yes all : for debug purpose)
  Int_t blabla = 2;
  
  //threshold (number of hits in a given asu at a given time (maximum=1536))
  Int_t threshold = 400;

  //energy threshold 
  Int_t E_threshold = 123;

  //nb plots
  Int_t nb_plots = 30;

  //timecut (time interval after trigger for tagging an event as square and remove it from efficency calculation)
  //  Int_t timecut = 50;

  //input data
  //string file = "04082011_2236_MGT.root";  
  //string file = "05082011_1055_MGT.root";  
  //string file = "05082011_2248_MGT.root";
  //
  //string file = "06082011_0002_MGT.root";  
  //string file = "06082011_0034_MGT.root";  
  //string file = "06082011_0105_MGT.root";  
  //string file = "06082011_0124_MGT.root";  
  //string file = "06082011_0152_MGT.root";  
  //string file = "06082011_0213_MGT.root";  
  //string file = "06082011_0231_MGT.root";  
  //string file = "06082011_0514_MGT.root";  
  //string file = "06082011_0536_MGT.root";  
  //string file = "06082011_0554_MGT.root";  
  //string file = "06082011_0652_MGT.root";  
  //string file = "06082011_0727_MGT.root";  
  //string file = "06082011_0829_MGT.root";  
  string file = "/lapp_data/LC/Detecteurs/MicroMegas/data/TB2012/SPS_H2_may_2012/Production/CaloHits/CaloHit_DHCAL_714448_I4_0.slcio.root";
 //-----------------
  

  //To see the plots in "live"
  TApplication *theApp;
  theApp = new TRint("App", 0, NULL);
  TCanvas *c0 = new TCanvas("c0", "Prototype m2", 100, 50, 1000, 600);
  TH2I *hxy = new TH2I("hxy","Hit position distribution (No time cut);y (cm);x (cm)",96,0,96,96,0,96);
  TH2I *hxy2 = new TH2I("hxy2","Hit position distribution (at the time where it is square);y (cm);x (cm)",96,0,96,96,0,96);
  TH2I *hxy3 = new TH2I("hxy3","Hit position distribution (at t_square - 3 bcid);y (cm);x (cm)",96,0,96,96,0,96);
  TH2I *hxy4 = new TH2I("hxy4","Hit position distribution (at t_square - 2 bcid);y (cm);x (cm)",96,0,96,96,0,96);       
  TH2I *hxy5 = new TH2I("hxy5","Hit position distribution (at t_square - 1 bcid);y (cm);x (cm)",96,0,96,96,0,96);
  TH2I *hxy6 = new TH2I("hxy6","Hit position distribution (at t_square + 1 bcid);y (cm);x (cm)",96,0,96,96,0,96);
  TH2I *hxy7 = new TH2I("hxy7","Hit position distribution (at t_square + 2 bcid);y (cm);x (cm)",96,0,96,96,0,96);
  TH2I *hxy8 = new TH2I("hxy8","Hit position distribution (at t_square + 3 bcid);y (cm);x (cm)",96,0,96,96,0,96);
        
  
  //our data file
  TFile *f = new TFile(file.c_str());
  if(blabla==1||blabla==2||blabla==3){
    cout<<endl;
    cout<<"reading File : "<<file<<endl;

  }
  
  //key iterator
  TIter nextkey(f->GetListOfKeys());
  TKey *key;
  UInt_t counter = 0;
   
  if(blabla==1||blabla==2||blabla==3) cout<<endl;
  

  //loop on all the keys
  while (key = (TKey*)nextkey()) {
    
    counter++;
    if(blabla==1||blabla==2||blabla==3) cout<<"key number = "<<counter<<endl;
    
    TTree *tree = (TTree*)key->ReadObj();                
    MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");

    Int_t nEvent = tree->GetEntries();
    if(blabla==1||blabla==2||blabla==3) cout <<"Nevent = "<<nEvent<<endl;
    
    MTEvent *evt = new MTEvent();
    
    TBranch *branch= tree->GetBranch("MTEvent");
    branch->SetAddress(&evt);
    if(blabla==2||blabla==3) f->Print();
    
    //at each new key, declare the event counters and indicators
    bool IsSquare = 0;
    Int_t list_square_events[nEvent][4];
    //    Int_t list_square_events_timecut[nEvent][4];
    UInt_t list_event[nEvent];
    //    UInt_t list_event_timecut[nEvent];
    UInt_t nb_tot = 0;
    //    UInt_t nb_timecut = 0;
    UInt_t begin_t = 0;
    UInt_t end_t = 0;

    
    //loop on all the events of the key
    for (int i=0;i<nEvent-1;i++){
      
      tree->GetEntry(i);
      
      UInt_t nchannel = evt->GetNchannel();
      //if(blabla==3) cout<<"Event "<<i+1<<" / "<<nEvent<<"; nb of channels=  "<<nchannel<<endl;

      //at each event declare the hits counter
      map<unsigned int,map<unsigned int,int> > channels_map;// key1=asu (board_id), key2=hit_time      
      
      //at each event declare the time interval
      Int_t tmin=0;
      Int_t tmax=0;
      
      //at each event reset the square indicator
      IsSquare=0;

      //at each event declare the time where the event is square
      UInt_t SquareTime=0;
      
      //loop on all the channels of the event
      for (int j=0;j<nchannel;j++){

        //this channel
        MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);

        //at each channel get board_id and hit_time
        UInt_t hit_time =  channel->GetBcIdAbs() - (channel->GetBcIdDif() - channel->GetBcIdHit() );
        UInt_t board_id = channel->GetBoardId();
        int DigiEng = channel->GetDigitalValue();
        //cout<<"Channel : "<<j<<"; hit on asu "<<board_id<<"; at time "<<hit_time<<"; with DigiEng "<<DigiEng<<endl;
        
        //get tmin and tmax
        if(j==0) {
          tmin=hit_time;
          tmax=hit_time;
        }
        if(hit_time<tmin) tmin=hit_time;
        if(hit_time>tmax) tmax=hit_time;

      //get begin_t and end_t
        if(i==0 && j==0) {
          begin_t=hit_time;
          end_t=hit_time;
        }
        if(hit_time<begin_t) begin_t=hit_time;
        if(hit_time>end_t) end_t=hit_time;

        //upgrade the hits counter of the good board at the good time
        if(E_threshold==123) if(DigiEng>=1) add(channels_map, board_id, hit_time);
        else if(DigiEng>=E_threshold) add(channels_map, board_id, hit_time);
      }    
      
      unsigned int timeold=0;
      
      //Loop on each board, at each time, to identify if this event is square
      map<unsigned int, map<unsigned int,int> >::iterator iter1;
      for(iter1=channels_map.begin();iter1!=channels_map.end();iter1++){
        
        unsigned int azu=iter1->first;
        map<unsigned int,int> hits_per_time=iter1->second;
        
        map<unsigned int,int>::iterator iter2;
        for(iter2=hits_per_time.begin();iter2!=hits_per_time.end();iter2++){
          
          unsigned int time = iter2->first;
          int nb_hits = iter2->second;

          //if this is a square event
          if (nb_hits>=threshold) {

            //set square event indicator to 1
            IsSquare=1;

            //set the time of this event
            SquareTime=time;

            if(blabla==3) cout <<" -> event: "<<i+1<<" = square in asu no "<<azu<<" at time "<<time<<" with n_hits = "<<nb_hits<<" DELTA T= "<<time-timeold<<endl;
            timeold=time;

            //update the event square lists :event, asu, time, nb_hits
            list_square_events[nb_tot][0] = i+1;
            list_square_events[nb_tot][1] = azu;
            list_square_events[nb_tot][2] = time;
            list_square_events[nb_tot][3] = nb_hits;
//          if(time>tmax-timecut){
//            list_square_events_timecut[nb_timecut][0] = i+1;
//            list_square_events_timecut[nb_timecut][1] = azu;
//            list_square_events_timecut[nb_timecut][2] = time;
//            list_square_events_timecut[nb_timecut][3] = nb_hits;
//          }
            
            //update the basic event square counters
            nb_tot++;
            //      if(time>tmax-timecut) nb_timecut++;
            
            //update the basic square event lists
            list_event[nb_tot-1] = i+1;   
            //      if(time>tmax-timecut) list_event_timecut[nb_tot-1] = i+1;   
          }
        }
      }

      //if we want to plot this event
      if (IsSquare && plot){

        //loop on all the channels of the event
        for (int j=0;j<nchannel;j++){
          
          //this channel
          MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
          
          //at each channel get board_id and hit_time
          UInt_t board_id = channel->GetBoardId();
          UInt_t hit_time =  channel->GetBcIdAbs() - (channel->GetBcIdDif() - channel->GetBcIdHit() );
          int DigiEng = channel->GetDigitalValue();
        
          //fill the histogramms
          if(E_threshold==123){
            if(DigiEng>=1){
              hxy->Fill(channel->GetY()/10000,channel->GetX()/10000);
              if(hit_time==SquareTime) hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000);
              if(hit_time==SquareTime-3) hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000);
              if(hit_time==SquareTime-2) hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000);
              if(hit_time==SquareTime-1) hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000);
              if(hit_time==SquareTime+1) hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000);
              if(hit_time==SquareTime+2) hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000);
              if(hit_time==SquareTime+3) hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000);
            }
            if(DigiEng>=2){
              hxy->Fill(channel->GetY()/10000,channel->GetX()/10000);
              hxy->Fill(channel->GetY()/10000,channel->GetX()/10000);
              if(hit_time==SquareTime) {
                //cout <<"2 ";
                hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000);
              }
              if(hit_time==SquareTime-3) {
                hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000);
              }
              if(hit_time==SquareTime-2) {
                hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000);
              }
              if(hit_time==SquareTime-1) {
                hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000);
              }
              if(hit_time==SquareTime+1) {
                hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000);
              }
              if(hit_time==SquareTime+2) {
                hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000);
              }
              if(hit_time==SquareTime+3) {
                hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000);
              }
            }
            if(DigiEng>=3){
              hxy->Fill(channel->GetY()/10000,channel->GetX()/10000);
              hxy->Fill(channel->GetY()/10000,channel->GetX()/10000);
              hxy->Fill(channel->GetY()/10000,channel->GetX()/10000);
              if(hit_time==SquareTime) {
                //cout <<"3 ";
                hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000);
              }
              if(hit_time==SquareTime-3) {
                hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000);
              }
              if(hit_time==SquareTime-2) {
                hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000);
              }
              if(hit_time==SquareTime-1) {
                hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000);
              }
              if(hit_time==SquareTime+1) {
                hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000);
              }
              if(hit_time==SquareTime+2) {
                hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000);
              }
              if(hit_time==SquareTime+3) {
                hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000);
                hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000);
              }
            }
          }
          else{
            hxy->Fill(channel->GetY()/10000,channel->GetX()/10000);
            if(hit_time==SquareTime) hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000);
            if(hit_time==SquareTime-3) hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000);
            if(hit_time==SquareTime-2) hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000);
            if(hit_time==SquareTime-1) hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000);
            if(hit_time==SquareTime+1) hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000);
            if(hit_time==SquareTime+2) hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000);
            if(hit_time==SquareTime+3) hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000);
          }
        }
      
        //plots
        c0->Divide(4,2);
        c0->cd(1);            
        hxy->GetXaxis()->SetTitle("X (cm)");
        hxy->GetYaxis()->SetTitle("Y (cm)");
        hxy->Draw("zcol");
        
        c0->cd(2);
        hxy3->GetXaxis()->SetTitle("X (cm)");
        hxy3->GetYaxis()->SetTitle("Y (cm)");
        hxy3->Draw("zcol");
        
        c0->cd(3);
        hxy4->GetXaxis()->SetTitle("X (cm)");
        hxy4->GetYaxis()->SetTitle("Y (cm)");
        hxy4->Draw("zcol");
        
        c0->cd(4);
        hxy5->GetXaxis()->SetTitle("X (cm)");
        hxy5->GetYaxis()->SetTitle("Y (cm)");
        hxy5->Draw("zcol");
        
        c0->cd(5);
        hxy2->GetXaxis()->SetTitle("X (cm)");
        hxy2->GetYaxis()->SetTitle("Y (cm)");
        hxy2->Draw("zcol");
        
        c0->cd(6);
        hxy6->GetXaxis()->SetTitle("X (cm)");
        hxy6->GetYaxis()->SetTitle("Y (cm)");
        hxy6->Draw("zcol");

        c0->cd(7);
        hxy7->GetXaxis()->SetTitle("X (cm)");
        hxy7->GetYaxis()->SetTitle("Y (cm)");
        hxy7->Draw("zcol");

        c0->cd(8);
        hxy8->GetXaxis()->SetTitle("X (cm)");
        hxy8->GetYaxis()->SetTitle("Y (cm)");
        hxy8->Draw("zcol");
        
        //save the plot
        char plotname[20] = "TH123_plot_evt";
        char buf[10];
        sprintf(buf,"%d",i+1);
        strcat(plotname, buf);
        strcat(plotname,".gif");
        c0->SaveAs(plotname);
        
        //clear the plot
        c0->Clear();
        
        //reset the histogramms
        hxy->Reset();
        hxy2->Reset();
        hxy3->Reset();
        hxy4->Reset();
        hxy5->Reset();
        hxy6->Reset();
        hxy7->Reset();
        hxy8->Reset();
      }
      if(nb_tot > nb_plots) break;
    }
    
    //Now the outputs: 
    if(blabla==2||blabla==3) cout<<endl;
    if(blabla==1||blabla==2||blabla==3) {

      //at each key print nb of evts square
      cout<<"In key "<<counter<<":"<<endl<<endl;;
      cout<<"   Nb square evt = "<<nb_tot<<endl;
      cout<<"   Nb square evt / Nevent tot = "<<1.0*nb_tot/nEvent*100.<<" %"<<endl;
      //      cout<<"   Nb square evt within the timecut (delta t < "<<timecut<<") = "<<nb_timecut<<endl;
      float rate = float(nb_tot)/(end_t-begin_t)*1000000000/200*60;
      cout<<"   Square event rate = "<<rate<<" event/min"<<endl<<endl;

      //at each key print print the liste of square events
      cout<<"List of square events: "<<endl;
      for(int it=0;it<nb_tot;it++) cout<<" "<<list_event[it];
      cout<<endl<<endl;

//       //at each key print print the liste of square events within the given timecut
//       cout<<"List of square events that are square within the timecut (delta t < "<<timecut<<")"<<endl;
//       for(int it=0;it<nb_timecut;it++) cout<<" "<<list_event_timecut[it];
//       cout<<endl<<endl;
    }

    if(blabla==2||blabla==3) {

      //if you want more details print also the position(asu) and time(memory order) of the square events
      cout<<"Hereafter I show the list of square events in each board/asu at each time and their number of hits:"<<endl;
      cout<<"Event, Board/asu, time, Nb_hits"<<endl;
      for(int it=0;it<nb_tot;it++) {
        for(int itt=0;itt<4;itt++) cout<<list_square_events[it][itt]<<";  ";
        cout<<endl;
      }
      cout<<endl;
      //      cout<<"Hereafter I show the list of square events that are square within the timecut (delta t < "<<timecut<<"):"<<endl;
      //      cout<<"Event, Board/asu, time, Nb_hits"<<endl;
      //      for(int it=0;it<nb_timecut;it++) {
      //        for(int itt=0;itt<4;itt++) cout<<list_square_events_timecut[it][itt]<<";  ";
      //        cout<<endl;
      //      }
      cout<<endl;
    }
    cout<<endl;
  }
  theApp->Run();
  delete theApp;
  return 0; 
}

---