Micromegas LAPP reconstruction

00001 //##################################
00002 //
00003 //--- Identify the border events ---
00004 //
00005 //##################################
00006 
00007 
00008 #include <map>
00009 #include <iostream>
00010 #include "TH1I.h"
00011 #include "TH2I.h"
00012 #include "TCanvas.h"
00013 #include "TFile.h"
00014 #include <TROOT.h>
00015 #include <TRint.h>
00016 #include <TKey.h>
00017 #include "TTree.h"
00018 #include <TApplication.h>
00019 
00020 #include "root/MTRun.hh"
00021 #include "root/MTEvent.hh"
00022 #include "root/MTChannel.hh"
00023 
00024 using namespace std;
00025 
00026 
00027 void add(map<unsigned int,map<unsigned long ,unsigned int> >&channels_map, unsigned int asu, unsigned long hit_time){
00028 
00029   if ( channels_map.find(asu) == channels_map.end()) { // if this asu doesn't exist yet
00030     
00031     map<unsigned long, unsigned int> nbHitMap;
00032     nbHitMap[hit_time]=1;
00033     channels_map.insert(make_pair(asu,nbHitMap));
00034   }
00035   else{ // if this asu already exist 
00036 
00037     map<unsigned long, unsigned int> &nbHitMap =  channels_map.find(asu)->second;
00038     if (nbHitMap.find(hit_time) == nbHitMap.end()){ //if this hit_time doesn't exist yet
00039 
00040       nbHitMap.insert(make_pair(hit_time,1)); 
00041     }
00042     else{ //if this abs time already exist
00043 
00044       nbHitMap[hit_time] = nbHitMap[hit_time]++;
00045     }
00046   }  
00047 }
00048 
00049 
00050 
00051 int main(void){  
00052   
00053   //-----------------
00054   //User's parameters
00055   
00056   //Do you want some info ? (0 = not at all, 1 = yes a little, 2 = yes a lot, 3 = yes all : for debug purpose)
00057   Int_t blabla = 1;
00058   
00059   //Do you want some plots ? (0 = not, 1 = yes)
00060   Int_t plot = 1;
00061   
00062   //threshold (number of hits in a given asu at a given time (maximum=160))
00063   Int_t threshold = 50;
00064 
00065   //timecut (time interval after trigger for tagging an event as border and remove it from efficency calculation)
00066   Int_t timecut = 50;
00067 
00068   //Micromegas chamber
00069   UInt_t chb = 1;
00070 
00071   //input data
00072   string file = "05082011_1055_MGT.root";  
00073   //string file = "06082011_0829_MGT.root";  
00074   //-----------------
00075   
00076 
00077   //To see the plots in "live"
00078   TApplication *theApp;
00079   theApp = new TRint("App", 0, NULL);
00080   TCanvas *c0 = new TCanvas("c0", "Prototype m2", 100, 50, 1000, 600);
00081   TH2I *hxy = new TH2I("hxy","Hit position distribution (No time cut);y (cm);x (cm)",96,0,96,96,0,96);
00082   TH2I *hxy2 = new TH2I("hxy2","Hit position distribution (within time cut);y (cm);x (cm)",96,0,96,96,0,96);
00083   TH2I *hxy3 = new TH2I("hxy3","Hit position distribution (at the time where it is border);y (cm);x (cm)",96,0,96,96,0,96);
00084         
00085   
00086   //our data file
00087   TFile *f = new TFile(file.c_str());
00088   if(blabla==1||blabla==2||blabla==3){
00089     cout<<endl;
00090     cout<<"reading File : "<<file<<endl;
00091   }
00092   
00093   //key iterator
00094   TIter nextkey(f->GetListOfKeys());
00095   TKey *key;
00096   UInt_t counter = 0;
00097    
00098   if(blabla==1||blabla==2||blabla==3) cout<<endl;
00099   
00100 
00101   //loop on all the keys
00102   while (key = (TKey*)nextkey()) {
00103     
00104     counter++;
00105     if(blabla==1||blabla==2||blabla==3) cout<<"key number = "<<counter<<endl;
00106     
00107     TTree *tree = (TTree*)key->ReadObj();                
00108     MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");
00109 
00110     Int_t nEvent = tree->GetEntries();
00111     if(blabla==1||blabla==2||blabla==3) cout <<"Nevent = "<<nEvent<<endl<<endl;
00112     
00113     MTEvent *evt = new MTEvent();
00114     
00115     TBranch *branch= tree->GetBranch("MTEvent");
00116     branch->SetAddress(&evt);
00117     
00118     //at each new key, declare the event counters and indicators
00119     bool IsBorder;
00120     Int_t Asu;
00121     ULong_t BorderTime;
00122     UInt_t nb_tot = 0;
00123     UInt_t nb_timecut = 0;
00124     ULong_t begin_t = 0;
00125     ULong_t end_t = 0;
00126 
00127     
00128     //loop on all the events of the key
00129     for (int i=0;i<nEvent-1;i++){
00130       
00131       tree->GetEntry(i);
00132       
00133       //an event can be border only if it is not square
00134       UInt_t A;
00135       UInt_t T;
00136       if(evt->IsSquare(300,A,T,chb)) {
00137         if(blabla==2||blabla==3) cout<<"evt "<<i+1<<" square"<<endl;
00138         continue;
00139       }
00140 
00141       UInt_t nchannel = evt->GetNchannel();
00142       if(blabla==2||blabla==3) cout<<"Event "<<i+1<<" / "<<nEvent<<"; nb of channels=  "<<nchannel<<endl;
00143 
00144       //at each event declare the hits counter
00145       map<unsigned int,map<unsigned long,unsigned int> > channels_map;// key1=asu (board_id), key2=hit_time      
00146       
00147       //at each event declare the time interval
00148       ULong_t tmin=2;
00149       ULong_t tmax=2;
00150       
00151       //at each event reset the border indicator
00152       IsBorder=0;
00153       Asu = -1;
00154       BorderTime = -1;
00155 
00156       //loop on all the channels of the event
00157       for (unsigned int j=0;j<nchannel;j++){
00158 
00159         //this channel
00160         MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00161 
00162         //at each channel get hit_time
00163         ULong_t hit_time =  channel->GetBcIdAbs() - (channel->GetBcIdDif() - channel->GetBcIdHit() );
00164         
00165         //get tmin and tmax of this event
00166         if(j==0) {
00167           tmin=hit_time;
00168           tmax=hit_time;
00169         }
00170         if(hit_time<tmin) tmin=hit_time;
00171         if(hit_time>tmax) tmax=hit_time;
00172 
00173         //get begin_t and end_t
00174         if(i==0 && j==0) {
00175           begin_t=hit_time;
00176           end_t=hit_time;
00177         }
00178         if(hit_time<begin_t) begin_t=hit_time;
00179         if(hit_time>end_t) end_t=hit_time;
00180 
00181         //if the hit is not on the borders of an asu go to next hit
00182         if(blabla==3) cout <<"X= "<<channel->GetX()<<"; Y= "<<channel->GetY()<<endl;
00183         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;
00184 
00185         //at each channel get board_id and hit_time
00186         UInt_t board_id = channel->GetBoardId();
00187         if(blabla==3) cout<<"Channel : "<<j<<"; hit on asu "<<board_id<<"; at time "<<hit_time<<endl;
00188         
00189         //upgrade the hits counter of the good board at the good time
00190         add(channels_map, board_id, hit_time);
00191       }      
00192       
00193       unsigned int timeold=0;
00194       
00195       //Loop on each board, at each time, to identify if this event is border
00196       map<unsigned int, map<unsigned long,unsigned int> >::iterator iter1;
00197       for(iter1=channels_map.begin();iter1!=channels_map.end();iter1++){
00198         
00199         unsigned int azu=iter1->first;
00200         map<unsigned long,unsigned int> hits_per_time=iter1->second;
00201         
00202         map<unsigned long,unsigned int>::iterator iter2;
00203         for(iter2=hits_per_time.begin();iter2!=hits_per_time.end();iter2++){
00204           
00205           unsigned long time = iter2->first;
00206           unsigned int nb_hits = iter2->second;
00207 
00208           //if this is a border event
00209           if (nb_hits>=threshold) {
00210 
00211             //set border event indicator to 1
00212             IsBorder = 1;
00213             Asu = azu;
00214             BorderTime = time;
00215 
00216             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;
00217           }
00218         }
00219       }
00220       //if we want to plot this event
00221       if (IsBorder && plot){
00222         
00223         TH1I *hdt = new TH1I("hdt","Hit time distribution (No time cut);N_hits;Time",tmax-tmin+3,-1,tmax-tmin+1);
00224         TH1I *hdt2 = new TH1I("hdt2","Hit time distribution (within time cut);N_hits;Time",timecut+3,tmax-timecut-1,tmax+1);
00225         TH1I *hdt3 = new TH1I("hdt3","Hit time distribution (around the time where it is border);N_hits;Time",23,BorderTime-tmin-11,BorderTime-tmin+11);
00226 
00227         //loop on all the channels of the event
00228         for (unsigned int j=0;j<nchannel;j++){
00229           
00230           //this channel
00231           MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
00232           
00233           //at each channel get board_id and hit_time
00234           UInt_t board_id = channel->GetBoardId();
00235           ULong_t hit_time =  channel->GetBcIdAbs() - (channel->GetBcIdDif() - channel->GetBcIdHit() );
00236           
00237           //fill the histogramms
00238           hxy->Fill(channel->GetY()/10000,channel->GetX()/10000);
00239           if(hit_time>tmax-timecut) hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000);
00240           if(hit_time==BorderTime) hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000);
00241           hdt->Fill(hit_time-tmin);
00242           if(hit_time>tmax-timecut) hdt2->Fill(hit_time-tmin);
00243           if(hit_time>BorderTime-100 && hit_time<BorderTime+500) hdt3->Fill(hit_time-tmin);
00244         }
00245       
00246         //plots
00247         c0->Divide(3,2);
00248         c0->cd(1);            
00249         hxy->GetXaxis()->SetTitle("X (cm)");
00250         hxy->GetYaxis()->SetTitle("Y (cm)");
00251         hxy->Draw("zcol");
00252         
00253         c0->cd(2);
00254         hxy2->GetXaxis()->SetTitle("X (cm)");
00255         hxy2->GetYaxis()->SetTitle("Y (cm)");
00256         hxy2->Draw("zcol");
00257         
00258         c0->cd(3);
00259         hxy3->GetXaxis()->SetTitle("X (cm)");
00260         hxy3->GetYaxis()->SetTitle("Y (cm)");
00261         hxy3->Draw("zcol");
00262         
00263         c0->cd(4);
00264         //hdt->SetAxisRange(-1,tmax-tmin+1);
00265         //hdt->SetMinimum(-1);
00266         //hdt->SetMaximum(tmax-tmin+1);
00267         hdt->GetXaxis()->SetTitle("time after trigger (BcId)");
00268         hdt->GetYaxis()->SetTitle("Nb hits");
00269         hdt->Draw();
00270         
00271         c0->cd(5);
00272         //hdt2->SetAxisRange(-1,timecut+1);
00273         //hdt2->SetMinimum(-1);
00274         //hdt2->SetMaximum(timecut+1);
00275         hdt2->GetXaxis()->SetTitle("time after trigger (BcId)");
00276         hdt2->GetYaxis()->SetTitle("Nb hits");
00277         hdt2->Draw();
00278         
00279         c0->cd(6);
00280         //hdt3->SetAxisRange(BorderTime-tmin-11,BorderTime-tmin+11);
00281         //hdt3->SetMinimum(BorderTime-tmin-11);
00282         //hdt3->SetMaximum(BorderTime-tmin+11);
00283         hdt3->GetXaxis()->SetTitle("time after trigger (BcId)");
00284         hdt3->GetYaxis()->SetTitle("Nb hits");
00285         hdt3->Draw();
00286         
00287         //save the plot
00288         //char plotname[100] = "/home1/girardd/Micromegas/Analyse/Damien/plots/plot_evt";
00289         char plotname[100] = "plot_evt";
00290         char buf[10];
00291         sprintf(buf,"%d",i+1);
00292         strcat(plotname, buf);
00293         strcat(plotname,".gif");
00294       
00295         //if(IsBorder==1) cout<<"DEBUG1 name:"<<plotname<<" max:"<<tmax<<" min:"<<tmin<<" cut:"<<timecut<<" bordertime:"<<BorderTime<<endl;
00296         
00297         c0->SaveAs(plotname);
00298                 
00299         //clear the plot
00300         c0->Clear();
00301         
00302         //reset the histogramms
00303         hxy->Reset();
00304         hxy2->Reset();
00305         hxy3->Reset();
00306         hdt->Delete();
00307         hdt2->Delete();
00308         hdt3->Delete();
00309       }
00310     }
00311     if(blabla==2||blabla==3) cout<<endl;      
00312   }
00313   if(blabla==1||blabla==2||blabla==3) cout<<endl;
00314   theApp->Run();
00315   delete theApp;
00316   return 0; 
00317 }
00318
/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/root/border_event.cpp
