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

Go to the documentation of this file.

//##################################
//
//--- Identify the square events ---
//
//##################################


#include <map>
#include <iostream>
#include "TFile.h"
#include <TROOT.h>
#include <TRint.h>
#include <TKey.h>
#include "TTree.h"

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

using namespace std;

int main(void){  
  
  //-----------------
  //User's parameters
  
  //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 = 3;
  
  //threshold (number of hits in a given asu at a given time (maximum=1536))
  Int_t threshold = 1000;

  //input data
  string file = "05082011_1055_MGT.root";  
  //-----------------
  

  //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<<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;
    Int_t Asu;
    Int_t AbsTime;

    
    //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 reset the square indicator
      IsSquare=0;
      Asu = -1;
      AbsTime = -1;

      //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;
        
        //upgrade the hits counter of the good board at the good time
        if ( channels_map.find(board_id) == channels_map.end()) { // if this asu doesn't exist yet
          
          map<unsigned int, int> nbHitMap;
          nbHitMap[hit_time]=1;
          channels_map.insert(make_pair(board_id,nbHitMap));
        }
        else{ // if this asu already exist 
          
          map<unsigned int, int> &nbHitMap =  channels_map.find(board_id)->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]++;
          }
        }  
      }      
      
      //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;
            Asu = azu;
            AbsTime = time;

            if(blabla==2) cout<<"*";
            if(blabla==3) cout <<" -> event: "<<i+1<<" = square in asu no "<<azu<<" at abstime "<<time<<" with n_hits = "<<nb_hits<<endl;
          }
        }
      }
    }
    if(blabla==2||blabla==3) cout<<endl;
  }
    if(blabla==1||blabla==2||blabla==3) cout<<endl;
  return 0; 
}
  

---