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

Go to the documentation of this file.

//##################################
//
//--- Identify the border 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 = 1;
  
  //threshold (number of hits in a given asu at a given time (maximum=160))
  Int_t threshold = 80;

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

  //Micromegas chamber
  UInt_t chb = 1;
  //-----------------
  

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

    
    //loop on all the events of the key
    for (int i=0;i<nEvent-1;i++){
      
      tree->GetEntry(i);
      
      //an event can be border only if it is not square
      UInt_t A;
      UInt_t T;
      if(evt->IsSquare(450,A,T,chb)) {
        if(blabla==2||blabla==3) cout<<"evt "<<i+1<<" square"<<endl;
        continue;
      }
      
      UInt_t nchannel = evt->GetNchannel();
      if(blabla==2||blabla==3) cout<<"Event "<<i+1<<" / "<<nEvent<<"; nb of channels=  "<<nchannel<<endl;
      
      //at each event declare the hits counter
      map<unsigned int,map<unsigned long,unsigned int> > channels_map;// key1=asu (board_id), key2=hit_time      
      
      //at each event reset the square indicator
      IsBorder=0;
      Asu = -1;
      AbsTime = -1;

      //loop on all the channels of the event
      for (unsigned int j=0;j<nchannel;j++){

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

        if (channel->GetChamberId()==chb)
          { 
            //at each channel get board_id and hit_time
            ULong_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 the position of this hit/channel
            //Int_t xpos = channel->GetX()/1000;
            //Int_t ypos = channel->GetY()/1000;
            if(blabla==3) cout <<"X= "<<channel->GetX()<<"; Y= "<<channel->GetY()<<endl;
            
            //if the hit is not on the borders of an asu go to next hit
            if (channel->GetX() != 0000 && channel->GetX() != 310000 && channel->GetX() != 320000 && channel->GetX() != 640000 && channel->GetX() != 650000 && channel->GetX() != 960000 && channel->GetY() != 000000 && channel->GetY() != 470000 && channel->GetY() != 480000 && channel->GetY() != 960000) continue;
            
            //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 long, unsigned int> nbHitMap;
              nbHitMap[hit_time]=1;
              channels_map.insert(make_pair(board_id,nbHitMap));
            }
            else{ // if this asu already exist 
              
              map<unsigned long, unsigned 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 border
        map<unsigned int, map<unsigned long,unsigned int> >::iterator iter1;
        for(iter1=channels_map.begin();iter1!=channels_map.end();iter1++){
          
          unsigned int azu=iter1->first;
          map<unsigned long,unsigned int> hits_per_time=iter1->second;
          
          map<unsigned long,unsigned int>::iterator iter2;
          for(iter2=hits_per_time.begin();iter2!=hits_per_time.end();iter2++){
            
            unsigned long time = iter2->first;
            unsigned int nb_hits = iter2->second;
            
            //if this is a border event
            if (nb_hits>=threshold) {
              
              //set square event indicator to 1
              IsBorder = 1;
              Asu = azu;
              AbsTime = time;
              
              if(blabla==1||blabla==2||blabla==3) cout <<" -> event: "<<i+1<<" = border 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; 
}

---