/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/Damien/evt_carre.C

Go to the documentation of this file.

//####################################################
//
//--- Identifier et compter les évennements carrés ---
//
//####################################################
{  
  
gROOT->Reset();
gROOT->SetStyle("Bold");
gStyle->SetPalette(1);
gSystem->Load("~/Micromegas/Analyse/trunk/lib/libMicro.so");

gSystem->Load("~/Micromegas/Analyse/trunk/src/analyse/Damien/load.C");

//using namespace std;

  //-----------------
  //User's parameters
  
  //which event do you want to plot (0 = no plot, a number = the event you want)
  bool plot = 1;
  
  //Do you want some info ? (0 = no, 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 = 500;

  //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;
    
//Là y'a un truc bizarre:
//si la déclaration de hxy est dans le if -> ca rame, si elle est hors du if -> ca tourne bien
//je comprends pas du tout pourquoi ??? 
  TH2I *hxy = new TH2I("hxy","Hit position distribution;y (cm);x (cm)",96,0,96,96,0,96);
  TH1I *hdt = new TH1I("hdt","Hit time distribution;N_hits;delta_t (/200ns)",128,0,128);
  if(plot){
    //TH2I *hxy = new TH2I("hxy","Hit position distribution;y (cm);x (cm)",96,0,96,96,0,96);
    TCanvas *c0 = new TCanvas("c0", "Prototype m2", 100, 50, 1000, 500);
    c0->Divide(2,1);
  }  
  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);
    
    UInt_t list_event[nEvent];
    
    //at each new key, declare the event counters
    std::map<unsigned int,std::map<unsigned int,int> > nb_square;// key1=asu (board_id), key2=abs_time      
    UInt_t nb_tot;

    
    //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
      std::map<unsigned int,std::map<unsigned int,int> > channels_map;// key1=asu (board_id), key2=abs_time      
      

      //loop on all the channels of the event
      for (int j=0;j<nchannel;j++){
        
        //create this channel
        MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);

        //at each channel declare board_id and abs_time
        UInt_t board_id = channel->GetBoardId();
        UInt_t abs_time =  channel->GetBcId_Abs() - (channel->GetBcId_Dif() - channel->GetBcId_Hit() );
        
        //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
          std::map<unsigned int, int> nbHitMap;
          nbHitMap[abs_time]=1;
          channels_map.insert(make_pair(board_id,nbHitMap));
        }
        else{ // if this asu already exist 
          std::map<unsigned int, int> &nbHitMap =  channels_map.find(board_id)->second;
          if (nbHitMap.find(abs_time) == nbHitMap.end()){ //if this abs_time doesn't exist yet
            nbHitMap.insert(make_pair(abs_time,1)); 
          }
          else{ //if this abs time already exist
            nbHitMap[abs_time] = nbHitMap[abs_time]++;
          }
        }  
        
        if(blabla==3) cout<<"Channel : "<<channel<<"; hit on asu "<<board_id<<"; at memory order "<<mem_order<<endl;
        
        //if we want to plot this event, fill the histogramm
        if (plot){
          hxy->Fill(channel->GetY()/10000,channel->GetX()/10000);
          hdt->Fill(abs_time);
        }
      }


      //On each board, at each time, identify the events which are square
      std::map<unsigned int, std::map<unsigned int,int> >::iterator iter1;
      for(iter1=channels_map.begin();iter1!=channels_map.end();iter1++){
        
        unsigned int azu=iter1->first;
        std::map<unsigned int,int> hits_per_abstime=iter1->second;
        
        std::map<unsigned int,int>::iterator iter2;
        for(iter2=hits_per_abstime.begin();iter2!=hits_per_abstime.end();iter2++){
          
          unsigned int abstime = iter2->first;
          int nb_hits = iter2->second;
          
          if (nb_hits>= threshold) {
            if(blabla==2||blabla==3) cout <<" -> event: "<<i+1<<" = square in asu no "<<azu<<" at abs_time "<<abstime<<" with n_hits = "<<nb_hits<<endl;

            //update the event square counter of the good board at the good time
              if ( nb_square.find(board_id) == nb_square.end()) { // if this asu doesn't exist yet
                std::map<unsigned int, int> nbHitMap;
                nbHitMap[abs_time]=1;
                nb_square.insert(make_pair(board_id,nbHitMap));
              }
              else{ // if this asu already exist 
                std::map<unsigned int, int> &nbHitMap =  nb_square.find(board_id)->second;
                if (nbHitMap.find(abs_time) == nbHitMap.end()){ //if this abs_time doesn't exist yet
                  nbHitMap.insert(make_pair(abs_time,1)); 
                }
                else{ //if this abs time already exist
                  nbHitMap[abs_time] = nbHitMap[abs_time]++;
                }
              }  
              
              //update the total event square counter of this key
              nb_tot++;
              
              //update the event list
              list_event[nb_tot-1] = i+1;
          }   
        }
      }
      
      //if we want to plot this event, let's do it now
      if (plot){
        c0->cd(1);
        hxy->SetTitle("Hits position");
        hxy.GetXaxis()->SetTitle("X (cm)");
        hxy.GetYaxis()->SetTitle("Y (cm)");
        hxy.GetYaxis()->SetLogy();
        hxy.GetYaxis()->SetGridy();
        hxy->Draw("zcol");

        c0->cd(2);
        hdt->SetTitle("Hits time");
        hdt.GetXaxis()->SetTitle("delta t (/ 200 ns)");
        hdt.GetYaxis()->SetTitle("Nb hits");
        hdt.GetYaxis()->SetLogy();
        hdt.GetYaxis()->SetGridy();
        hdt->Draw();
      }
    }
    

    //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<<" Nb tot evt square = "<<nb_tot<<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;
    }

    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 nb of event square in each board/asu at each time:"<<endl<<endl<<"Board/asu i"<<endl<<"memory order:";
      for (int ti=0;ti<128;ti++) cout<<" "<<ti;
      cout<<endl;

      std::map<unsigned int, std::map<unsigned int,int> >::iterator iter1;
      for(iter1=nb_square.begin();iter1!=nb_square.end();iter1++){
        
        unsigned int azu=iter1->first;
        std::map<unsigned int,int> nb_square_per_abstime=iter1->second;
        
        cout<<endl;
        cout<<"asu "<<azu<<endl;

        std::map<unsigned int,int>::iterator iter2;
        for(iter2=nb_square_per_abstime.begin();iter2!=nb_square_per_abstime.end();iter2++){
          
          unsigned int abstime = iter2->first;
          int nb_hits = iter2->second;
          
          cout<<" "<<nb_hits;
        }
        cout<<endl;
      }
      cout<<endl;
    }
    cout<<endl;
  }
  return 0; 
}

---