/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/root/square_event.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 = 1000;

  //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";  
  //-----------------
  

  //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 (within time cut);y (cm);x (cm)",96,0,96,96,0,96);
  TH2I *hxy3 = new TH2I("hxy3","Hit position distribution (at the time where it is square);y (cm);x (cm)",96,0,96,96,0,96);
  //TH1I *hdt = new TH1I("hdt","Hit time distribution (No time cut);N_hits;Time",10003,-1,10001);
  //TH1I *hdt2 = new TH1I("hdt2","Hit time distribution (within time cut);N_hits;Time",103,-1,101);
  //TH1I *hdt3 = new TH1I("hdt3","Hit time distribution (around the time where it is square);N_hits;Time",23,-1,21);
        
  
  //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);
    
    //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;

    
    //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();
        //if(blabla==3) cout<<"Channel : "<<j<<"; hit on asu "<<board_id<<"; at time "<<hit_time<<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;

        //upgrade the hits counter of the good board at the good time
        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
            //this.set->IsSquare=true
            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){

        TH1I *hdt = new TH1I("hdt","Hit time distribution (No time cut);N_hits;Time",tmax-tmin+3,-1,tmax-tmin+1);
        TH1I *hdt2 = new TH1I("hdt2","Hit time distribution (within time cut);N_hits;Time",timecut+3,tmax-timecut-1,tmax+1);
        TH1I *hdt3 = new TH1I("hdt3","Hit time distribution (around the time where it is square);N_hits;Time",23,SquareTime-tmin-11,SquareTime-tmin+11);

        //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() );
          
          //fill the histogramms
          hxy->Fill(channel->GetY()/10000,channel->GetX()/10000);
          if(hit_time>tmax-timecut) hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000);
          if(hit_time==SquareTime) hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000);
          hdt->Fill(hit_time-tmin);
          if(hit_time>tmax-timecut) hdt2->Fill(hit_time-tmin);
          if(hit_time>SquareTime-100 && hit_time<SquareTime+500) hdt3->Fill(hit_time-tmin);
        }
      
        //plots
        c0->Divide(3,2);
        c0->cd(1);            
        hxy->GetXaxis()->SetTitle("X (cm)");
        hxy->GetYaxis()->SetTitle("Y (cm)");
        hxy->Draw("zcol");
        
        c0->cd(2);
        hxy2->GetXaxis()->SetTitle("X (cm)");
        hxy2->GetYaxis()->SetTitle("Y (cm)");
        hxy2->Draw("zcol");
        
        c0->cd(3);
        hxy3->GetXaxis()->SetTitle("X (cm)");
        hxy3->GetYaxis()->SetTitle("Y (cm)");
        hxy3->Draw("zcol");
        
        c0->cd(4);
        //hdt->SetAxisRange(-1,tmax-tmin+1);
        //hdt->SetMinimum(-1);
        //hdt->SetMaximum(tmax-tmin+1);
        hdt->GetXaxis()->SetTitle("time after trigger (BcId)");
        hdt->GetYaxis()->SetTitle("Nb hits");
        hdt->Draw();
        
        c0->cd(5);
        //hdt2->SetAxisRange(-1,timecut+1);
        //hdt2->SetMinimum(-1);
        //hdt2->SetMaximum(timecut+1);
        hdt2->GetXaxis()->SetTitle("time after trigger (BcId)");
        hdt2->GetYaxis()->SetTitle("Nb hits");
        hdt2->Draw();
        
        c0->cd(6);
        //hdt3->SetAxisRange(SquareTime-tmin-11,SquareTime-tmin+11);
        //hdt3->SetMinimum(SquareTime-tmin-11);
        //hdt3->SetMaximum(SquareTime-tmin+11);
        hdt3->GetXaxis()->SetTitle("time after trigger (BcId)");
        hdt3->GetYaxis()->SetTitle("Nb hits");
        hdt3->Draw();
        
        //save the plot
        char plotname[20] = "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();
        //hdt->Reset();
        //hdt2->Reset();
        //hdt3->Reset();
        //hxy->Delete();
        //hxy2->Delete();
        //hxy3->Delete();
        hdt->Delete();
        hdt2->Delete();
        hdt3->Delete();
      }
    }
    
    //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 within the timecut (delta t < "<<timecut<<") = "<<nb_timecut<<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; 
}

---