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

Go to the documentation of this file.

//############################################
//
//--- Identify and count the chip 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 chip, unsigned int hit_time){

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

    map<unsigned int, int> &nbHitMap =  channels_map.find(chip)->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 chip at a given time (maximum=64))
  Int_t threshold = 40;

  //timecut (time interval after trigger for tagging an event as chip 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";  
  //-----------------
  

  //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 chip);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 chip);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 IsChip = 0;
    Int_t list_chip_events[nEvent][4];
    Int_t list_chip_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=chip (chip_id), key2=hit_time      
      
      //at each event declare the time interval
      Int_t tmin=0;
      Int_t tmax=0;
      
      //at each event reset the chip indicator
      IsChip=0;

      //at each event declare the time where the event is chip
      UInt_t ChipTime=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 chip_id and hit_time
        UInt_t hit_time =  channel->GetBcIdAbs() - (channel->GetBcIdDif() - channel->GetBcIdHit() );
        UInt_t chip_id = channel->GetChipId();
        //if(blabla==3) cout<<"Channel : "<<j<<"; hit on chip "<<chip_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;

      //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 chip at the good time
        add(channels_map, chip_id, hit_time);
      }    
      
      unsigned int timeold=0;
      
      //Loop on each chip, at each time, to identify if this event is chip
      map<unsigned int, map<unsigned int,int> >::iterator iter1;
      for(iter1=channels_map.begin();iter1!=channels_map.end();iter1++){
        
        unsigned int chiip=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 chip event
          if (nb_hits>=threshold) {

            //set chip event indicator to 1
            //this.set->IsChip=true
            IsChip=1;

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

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

            //update the event chip lists :event, chip, time, nb_hits
            list_chip_events[nb_tot][0] = i+1;
            list_chip_events[nb_tot][1] = chiip;
            list_chip_events[nb_tot][2] = time;
            list_chip_events[nb_tot][3] = nb_hits;
            if(time>tmax-timecut){
              list_chip_events_timecut[nb_timecut][0] = i+1;
              list_chip_events_timecut[nb_timecut][1] = chiip;
              list_chip_events_timecut[nb_timecut][2] = time;
              list_chip_events_timecut[nb_timecut][3] = nb_hits;
            }
            
            //update the basic event chip counters
            nb_tot++;
            if(time>tmax-timecut) nb_timecut++;
            
            //update the basic chip 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 (IsChip && 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 chip);N_hits;Time",23,ChipTime-tmin-11,ChipTime-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 chip_id and hit_time
          UInt_t chip_id = channel->GetChipId();
          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==ChipTime) 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>ChipTime-100 && hit_time<ChipTime+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(ChipTime-tmin-11,ChipTime-tmin+11);
        //hdt3->SetMinimum(ChipTime-tmin-11);
        //hdt3->SetMaximum(ChipTime-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 chip
      cout<<"In key "<<counter<<":"<<endl<<endl;;
      cout<<"   Nb chip evt = "<<nb_tot<<endl;
      cout<<"   Nb chip evt / Nevent tot = "<<1.0*nb_tot/nEvent*100.<<" %"<<endl;
      cout<<"   Nb chip evt within the timecut (delta t < "<<timecut<<") = "<<nb_timecut<<endl;
      float rate = float(nb_tot)/(end_t-begin_t)*1000000000/200*60;
      cout<<"   Chip event rate = "<<rate<<" event/min"<<endl<<endl;

      //at each key print print the liste of chip events
      cout<<"List of chip 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 chip events within the given timecut
      cout<<"List of chip events that are chip 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(chip) and time(memory order) of the chip events
      cout<<"Hereafter I show the list of chip events in each chip/chip at each time and their number of hits:"<<endl;
      cout<<"Event, chip, time, Nb_hits"<<endl;
      for(int it=0;it<nb_tot;it++) {
        for(int itt=0;itt<4;itt++) cout<<list_chip_events[it][itt]<<";  ";
        cout<<endl;
      }
      cout<<endl;
      cout<<"Hereafter I show the list of chip events that are chip within the timecut (delta t < "<<timecut<<"):"<<endl;
      cout<<"Event, chip, time, Nb_hits"<<endl;
      for(int it=0;it<nb_timecut;it++) {
        for(int itt=0;itt<4;itt++) cout<<list_chip_events_timecut[it][itt]<<";  ";
        cout<<endl;
      }
      cout<<endl;
    }
    cout<<endl;
  }
  theApp->Run();
  delete theApp;
  return 0; 
}

---