Micromegas LAPP reconstruction

00001 //############################################
00002 //
00003 //--- Identify and count the square events ---
00004 //
00005 //############################################
00006 
00007 
00008 #include <map>
00009 #include "TKey.h"
00010 #include <TProfile.h>
00011 #include <TROOT.h>
00012 #include <TRint.h>
00013 #include <TKey.h>
00014 #include <TList.h>
00015 #include <TTree.h>
00016 #include <TFile.h>
00017 #include <TChain.h>
00018 #include "TChainElement.h"
00019 #include <TMath.h>
00020 #include <stdlib.h>
00021 #include <TSystem.h>
00022 #include <iostream>
00023 #include <sstream>
00024 #include "Log.hh"
00025 #include "MicroException.hh"
00026 #include <string>
00027 #include "mysql/Mysql.hh"
00028 #include "root/CaloRun.hh"
00029 #include "root/CaloEvent.hh"
00030 #include "root/CaloHit.hh"
00031 #include "root/MTRun.hh"
00032 #include "root/MTEvent.hh"
00033 #include "root/MTChannel.hh"
00034 #include <TCanvas.h>
00035 #include "TGraphErrors.h"
00036 #include "TGraphAsymmErrors.h"
00037 #include "TH1I.h"
00038 #include "TH2I.h"
00039 #include <TApplication.h>
00040 
00041 using namespace std;
00042 
00043 const int maxTrees = 200;  // How many trees to read
00044 
00045 
00046 void add(map<unsigned int,map<unsigned int ,int> >&channels_map, unsigned int asu, unsigned int hit_time){
00047 
00048   if ( channels_map.find(asu) == channels_map.end()) { // if this asu doesn't exist yet
00049     
00050     map<unsigned int, int> nbHitMap;
00051     nbHitMap[hit_time]=1;
00052     channels_map.insert(make_pair(asu,nbHitMap));
00053   }
00054   else{ // if this asu already exist 
00055 
00056     map<unsigned int, int> &nbHitMap =  channels_map.find(asu)->second;
00057     if (nbHitMap.find(hit_time) == nbHitMap.end()){ //if this hit_time doesn't exist yet
00058 
00059       nbHitMap.insert(make_pair(hit_time,1)); 
00060     }
00061     else{ //if this abs time already exist
00062 
00063       nbHitMap[hit_time] = nbHitMap[hit_time]++;
00064     }
00065   }  
00066 }
00067 
00068 
00069 
00070 int main(int argc, char**argv){
00071 
00072   if ( argc !=2  ) {
00073     FILE_LOG(logERROR)  << "usage: example rootFile " << endl;
00074     exit(1);
00075   }
00076 
00077   //-----------------
00078   //User's parameters
00079   
00080   //do you want to plot (0 = no plot, 1 = plots)
00081   bool plot = 1;
00082   
00083   //Do you want some info ? (0 = not at all, 1 = yes a little, 2 = yes a lot, 3 = yes all : for debug purpose)
00084   Int_t blabla = 2;
00085   
00086   //threshold (number of hits in a given asu at a given time (maximum=1536))
00087   Int_t threshold = 400;
00088 
00089   //energy threshold 
00090   Int_t E_threshold = 123;
00091 
00092   //nb plots
00093   Int_t nb_plots = 30;
00094 
00095   //timecut (time interval after trigger for tagging an event as square and remove it from efficency calculation)
00096   //  Int_t timecut = 50;
00097 
00098   //input data
00099   //string file = "/lapp_data/LC/Detecteurs/MicroMegas/data/TB2012/SPS_H2_may_2012/Production/CaloHits/CaloHit_DHCAL_714448_I4_0.slcio.root";
00100  //-----------------
00101   
00102 
00103   //To see the plots in "live"
00104   TApplication *theApp;
00105   theApp = new TRint("App", 0, NULL);
00106   TCanvas *c0 = new TCanvas("c0", "Prototype m2", 50, 50, 1200, 400);
00107   TH2I *hxy = new TH2I("hxy","Hit position distribution (No time cut);y (cm);x (cm)",96,0,96,96,0,96);
00108   TH2I *hxy2 = new TH2I("hxy2","Hit position distribution (at the time where it is square);y (cm);x (cm)",96,0,96,96,0,96);
00109   TH2I *hxy3 = new TH2I("hxy3","Hit position distribution (at t_square - 3 bcid);y (cm);x (cm)",96,0,96,96,0,96);
00110   TH2I *hxy4 = new TH2I("hxy4","Hit position distribution (at t_square - 2 bcid);y (cm);x (cm)",96,0,96,96,0,96);       
00111   TH2I *hxy5 = new TH2I("hxy5","Hit position distribution (at t_square - 1 bcid);y (cm);x (cm)",96,0,96,96,0,96);
00112   TH2I *hxy6 = new TH2I("hxy6","Hit position distribution (at t_square + 1 bcid);y (cm);x (cm)",96,0,96,96,0,96);
00113   TH2I *hxy7 = new TH2I("hxy7","Hit position distribution (at t_square + 2 bcid);y (cm);x (cm)",96,0,96,96,0,96);
00114   TH2I *hxy8 = new TH2I("hxy8","Hit position distribution (at t_square + 3 bcid);y (cm);x (cm)",96,0,96,96,0,96);
00115   TH2I *hxy9 = new TH2I("hxy9","Hit position distribution (at t_square - 6 bcid);y (cm);x (cm)",96,0,96,96,0,96);
00116   TH2I *hxy10 = new TH2I("hxy10","Hit position distribution (at t_square - 5 bcid);y (cm);x (cm)",96,0,96,96,0,96);
00117   TH2I *hxy11 = new TH2I("hxy11","Hit position distribution (at t_square - 4 bcid);y (cm);x (cm)",96,0,96,96,0,96);
00118   TH2I *hxy12 = new TH2I("hxy12","Hit position distribution (at t_square + 4 bcid);y (cm);x (cm)",96,0,96,96,0,96);
00119   TH2I *hxy13 = new TH2I("hxy13","Hit position distribution (at t_square + 5 bcid);y (cm);x (cm)",96,0,96,96,0,96);
00120   TH2I *hxy14 = new TH2I("hxy14","Hit position distribution (at t_square + 6 bcid);y (cm);x (cm)",96,0,96,96,0,96);
00121         
00122   
00123   int NbTrees = 0;
00124 
00125   //data reading structure
00126   string rootName;
00127   rootName.assign(argv[1]);
00128   TChain chain("T");
00129   for (int i =0; i<5 ; i++){
00130     TString SubFile = Form("_I%i",i);
00131     TString fname = Form("/lapp_data/LC/Detecteurs/MicroMegas/data/TB2012/SPS_H2_may_2012/Production/CaloHits/CaloHit_DHCAL_"+rootName+SubFile+"_0.slcio.root");
00132   chain.Add(fname);
00133   cout << fname << endl;
00134   }
00135 
00136   TObjArray *fileElements=chain.GetListOfFiles();
00137   TIter next(fileElements);
00138   TChainElement *chEl=0;
00139   int NchainFiles=0;
00140   while (( chEl=(TChainElement*)next() )) 
00141     {
00142       TFile f(chEl->GetTitle());
00143       TIter nextkey(f.GetListOfKeys());
00144       TKey *key;  
00145       UInt_t counter = 0;
00146       NchainFiles++;
00147       cout << " Reading a new file in the chain " << NchainFiles << endl;
00148       while (key = (TKey*)nextkey()) 
00149         {
00150           TString keyname = key->GetName();
00151           if ((keyname == "hResidualX") || (keyname == "hResidualY")){}
00152           else
00153             {
00154               TTree *tree = (TTree*)key->ReadObj();                
00155               cout << "---------------------------- NEW TTREE -----------------"<< endl;
00156               if(blabla==3) f.Print();
00157               NbTrees++;
00158               if( NbTrees > maxTrees ) continue;
00159               CaloRun *run=NULL;
00160 
00161               try {
00162                 CaloRun *run=dynamic_cast<CaloRun *>(tree->GetUserInfo()->First());
00163                 if( run != NULL ){cout << "Run[" << run->GetRunId() <<  "], temperature[" << run->GetTemperature() <<  "] pression[" << run->GetPressure() << "] " << endl ; }
00164 
00165               }
00166               catch ( ... ) {}
00167               CaloEvent *evt =  new CaloEvent();
00168               TBranch *branch= tree->GetBranch("CaloEvent");
00169               branch->SetAddress(&evt);
00170 
00171               CaloHit* channel =NULL;
00172               int nEvent = tree->GetEntries();
00173               cout << " Nb of events to read " << nEvent << endl;
00174 
00175               //at each new key, declare the event counters and indicators
00176               bool IsSquare = 0;
00177               Int_t list_square_events[nEvent][4];
00178               UInt_t list_event[nEvent];
00179               UInt_t nb_tot = 0;
00180               UInt_t begin_t = 0;
00181               UInt_t end_t = 0;
00182 
00183               for (int i=0; i<nEvent-1; i++){
00184               //for(int i=0 ; i<10 ; i++){
00185 
00186                 tree->GetEntry(i);
00187                 int nchannel = evt->GetNHits();
00188                 
00189                 //if(blabla==3) cout<<"Event "<<i+1<<" / "<<nEvent<<"; nb of channels=  "<<nchannel<<endl;
00190                 
00191                 //at each event declare the hits counter
00192                 map<unsigned int,map<unsigned int,int> > channels_map;// key1=asu (board_id), key2=hit_time      
00193                 
00194                 //at each event declare the time interval
00195                 Int_t tmin=0;
00196                 Int_t tmax=0;
00197                 
00198                 //at each event reset the square indicator
00199                 IsSquare=0;
00200                 
00201                 //at each event declare the time where the event is square
00202                 UInt_t SquareTime=0;
00203 
00204                 for(int j=0 ; j<nchannel ; j++){
00205                   //for(int j=0 ; j<10 ; j++){
00206                       channel = (CaloHit*)evt->GetHits()->UncheckedAt(j);
00207                       //if(channel->GetLayer()<48) continue;
00208                       if(channel->GetLayer()==48 || channel->GetLayer()==49){
00209 
00210                         //cout<<channel->GetLayer()<<" ";
00211                         
00212                         //at each channel get board_id and hit_time
00213                         //UInt_t hit_time =  channel->GetDeltaT();
00214                         UInt_t hit_time =  evt->GetDeltaTmax() - channel->GetDeltaT() ;
00215                         UInt_t board_id = channel->GetSolftId().GetBoardId();
00216                         int DigiEng = channel->GetThreshold();
00217                         //cout<<"Channel : "<<j<<"; hit on asu "<<board_id<<"; at time "<<hit_time<<"; with DigiEng "<<DigiEng<<endl;
00218                         
00219                         //get tmin and tmax
00220                         if(j==0) {
00221                           tmin=hit_time;
00222                           tmax=hit_time;
00223                         }
00224                         if(hit_time<tmin) tmin=hit_time;
00225                         if(hit_time>tmax) tmax=hit_time;
00226                         
00227                         //get begin_t and end_t
00228                         if(i==0 && j==0) {
00229                           begin_t=hit_time;
00230                           end_t=hit_time;
00231                         }
00232                         if(hit_time<begin_t) begin_t=hit_time;
00233                         if(hit_time>end_t) end_t=hit_time;
00234                         
00235                         //upgrade the hits counter of the good board at the good time
00236                         if(E_threshold==123) if(DigiEng>=1) add(channels_map, board_id, hit_time);
00237                         else if(DigiEng>=E_threshold) add(channels_map, board_id, hit_time);
00238                       } 
00239                   }
00240                 unsigned int timeold=0;
00241                 
00242                 //Loop on each board, at each time, to identify if this event is square
00243                 map<unsigned int, map<unsigned int,int> >::iterator iter1;
00244                 for(iter1=channels_map.begin();iter1!=channels_map.end();iter1++){
00245                   
00246                   unsigned int azu=iter1->first;
00247                   map<unsigned int,int> hits_per_time=iter1->second;
00248                   
00249                   map<unsigned int,int>::iterator iter2;
00250                   for(iter2=hits_per_time.begin();iter2!=hits_per_time.end();iter2++){
00251                     
00252                     unsigned int time = iter2->first;
00253                     int nb_hits = iter2->second;
00254                     
00255                     //if this is a square event
00256                     if (nb_hits>=threshold) {
00257                       
00258                       //set square event indicator to 1
00259                       IsSquare=1;
00260                       
00261                       //set the time of this event
00262                       SquareTime=time;
00263                       
00264                       if(blabla==3) cout <<" -> event: "<<i+1<<" = square in asu no "<<azu<<" at time "<<time<<" with n_hits = "<<nb_hits<<" DELTA T= "<<time-timeold<<endl;
00265                       timeold=time;
00266                       
00267                       //update the event square lists :event, asu, time, nb_hits
00268                       list_square_events[nb_tot][0] = i+1;
00269                       list_square_events[nb_tot][1] = azu;
00270                       list_square_events[nb_tot][2] = time;
00271                       list_square_events[nb_tot][3] = nb_hits;
00272                       //            if(time>tmax-timecut){
00273                       //              list_square_events_timecut[nb_timecut][0] = i+1;
00274                       //              list_square_events_timecut[nb_timecut][1] = azu;
00275                       //              list_square_events_timecut[nb_timecut][2] = time;
00276                       //              list_square_events_timecut[nb_timecut][3] = nb_hits;
00277                       //            }
00278                       
00279                       //update the basic event square counters
00280                       nb_tot++;
00281                       //            if(time>tmax-timecut) nb_timecut++;
00282                       
00283                       //update the basic square event lists
00284                       list_event[nb_tot-1] = i+1;   
00285                       //            if(time>tmax-timecut) list_event_timecut[nb_tot-1] = i+1;   
00286                     }
00287                   }
00288                 }
00289                 
00290                 //if we want to plot this event
00291                 if (IsSquare && plot){
00292                   
00293                   //loop on all the channels of the event
00294                   for (int j=0;j<nchannel;j++){
00295                     
00296                     //this channel
00297                     channel = (CaloHit*)evt->GetHits()->UncheckedAt(j);
00298                     
00299                     //at each channel get board_id and hit_time
00300                     //UInt_t hit_time =  channel->GetDeltaT();
00301                     UInt_t hit_time =  evt->GetDeltaTmax() - channel->GetDeltaT() ;
00302                     UInt_t board_id = channel->GetSolftId().GetBoardId();
00303                     int DigiEng = channel->GetThreshold();
00304                     
00305                     //fill the histogramms
00306                     if(E_threshold==123){
00307                       if(DigiEng>=1){
00308                         hxy->Fill(channel->GetY()/10000,channel->GetX()/10000);
00309                         if(hit_time==SquareTime) hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000);
00310                         if(hit_time==SquareTime-3) hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000);
00311                         if(hit_time==SquareTime-2) hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000);
00312                         if(hit_time==SquareTime-1) hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000);
00313                         if(hit_time==SquareTime+1) hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000);
00314                         if(hit_time==SquareTime+2) hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000);
00315                         if(hit_time==SquareTime+3) hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000);
00316                         if(hit_time==SquareTime-6) hxy9->Fill(channel->GetY()/10000,channel->GetX()/10000);
00317                         if(hit_time==SquareTime-5) hxy10->Fill(channel->GetY()/10000,channel->GetX()/10000);
00318                         if(hit_time==SquareTime-4) hxy11->Fill(channel->GetY()/10000,channel->GetX()/10000);
00319                         if(hit_time==SquareTime+4) hxy12->Fill(channel->GetY()/10000,channel->GetX()/10000);
00320                         if(hit_time==SquareTime+5) hxy13->Fill(channel->GetY()/10000,channel->GetX()/10000);
00321                         if(hit_time==SquareTime+6) hxy14->Fill(channel->GetY()/10000,channel->GetX()/10000);
00322                       }
00323                       if(DigiEng>=2){
00324                         hxy->Fill(channel->GetY()/10000,channel->GetX()/10000);
00325                         hxy->Fill(channel->GetY()/10000,channel->GetX()/10000);
00326                         if(hit_time==SquareTime) {
00327                           hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000);
00328                           hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000);
00329                         }
00330                         if(hit_time==SquareTime-3) {
00331                           hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000);
00332                           hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000);
00333                         }
00334                         if(hit_time==SquareTime-2) {
00335                           hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000);
00336                           hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000);
00337                         }
00338                         if(hit_time==SquareTime-1) {
00339                           hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000);
00340                           hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000);
00341                         }
00342                         if(hit_time==SquareTime+1) {
00343                           hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000);
00344                           hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000);
00345                         }
00346                         if(hit_time==SquareTime+2) {
00347                           hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000);
00348                           hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000);
00349                         }
00350                         if(hit_time==SquareTime+3) {
00351                           hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000);
00352                           hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000);
00353                         }
00354                         if(hit_time==SquareTime-6) {
00355                           hxy9->Fill(channel->GetY()/10000,channel->GetX()/10000);
00356                           hxy9->Fill(channel->GetY()/10000,channel->GetX()/10000);
00357                         }
00358                         if(hit_time==SquareTime-5) {
00359                           hxy10->Fill(channel->GetY()/10000,channel->GetX()/10000);
00360                           hxy10->Fill(channel->GetY()/10000,channel->GetX()/10000);
00361                         }
00362                         if(hit_time==SquareTime-4) {
00363                           hxy11->Fill(channel->GetY()/10000,channel->GetX()/10000);
00364                           hxy11->Fill(channel->GetY()/10000,channel->GetX()/10000);
00365                         }
00366                         if(hit_time==SquareTime+4) {
00367                           hxy12->Fill(channel->GetY()/10000,channel->GetX()/10000);
00368                           hxy12->Fill(channel->GetY()/10000,channel->GetX()/10000);
00369                         }
00370                         if(hit_time==SquareTime+5) {
00371                           hxy13->Fill(channel->GetY()/10000,channel->GetX()/10000);
00372                           hxy13->Fill(channel->GetY()/10000,channel->GetX()/10000);
00373                         }
00374                         if(hit_time==SquareTime+6) {
00375                           hxy14->Fill(channel->GetY()/10000,channel->GetX()/10000);
00376                           hxy14->Fill(channel->GetY()/10000,channel->GetX()/10000);
00377                         }
00378                       }
00379                       if(DigiEng>=3){
00380                         hxy->Fill(channel->GetY()/10000,channel->GetX()/10000);
00381                         hxy->Fill(channel->GetY()/10000,channel->GetX()/10000);
00382                         hxy->Fill(channel->GetY()/10000,channel->GetX()/10000);
00383                         if(hit_time==SquareTime) {
00384                           hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000);
00385                           hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000);
00386                           hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000);
00387                         }
00388                         if(hit_time==SquareTime-3) {
00389                           hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000);
00390                           hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000);
00391                           hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000);
00392                         }
00393                         if(hit_time==SquareTime-2) {
00394                           hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000);
00395                           hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000);
00396                           hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000);
00397                         }
00398                         if(hit_time==SquareTime-1) {
00399                           hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000);
00400                           hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000);
00401                           hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000);
00402                         }
00403                         if(hit_time==SquareTime+1) {
00404                           hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000);
00405                           hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000);
00406                           hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000);
00407                         }
00408                         if(hit_time==SquareTime+2) {
00409                           hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000);
00410                           hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000);
00411                           hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000);
00412                         }
00413                         if(hit_time==SquareTime+3) {
00414                           hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000);
00415                           hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000);
00416                           hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000);
00417                         }
00418                         if(hit_time==SquareTime-6) {
00419                           hxy9->Fill(channel->GetY()/10000,channel->GetX()/10000);
00420                           hxy9->Fill(channel->GetY()/10000,channel->GetX()/10000);
00421                           hxy9->Fill(channel->GetY()/10000,channel->GetX()/10000);
00422                         }
00423                         if(hit_time==SquareTime-5) {
00424                           hxy10->Fill(channel->GetY()/10000,channel->GetX()/10000);
00425                           hxy10->Fill(channel->GetY()/10000,channel->GetX()/10000);
00426                           hxy10->Fill(channel->GetY()/10000,channel->GetX()/10000);
00427                         }
00428                         if(hit_time==SquareTime-4) {
00429                           hxy11->Fill(channel->GetY()/10000,channel->GetX()/10000);
00430                           hxy11->Fill(channel->GetY()/10000,channel->GetX()/10000);
00431                           hxy11->Fill(channel->GetY()/10000,channel->GetX()/10000);
00432                         }
00433                         if(hit_time==SquareTime+4) {
00434                           hxy12->Fill(channel->GetY()/10000,channel->GetX()/10000);
00435                           hxy12->Fill(channel->GetY()/10000,channel->GetX()/10000);
00436                           hxy12->Fill(channel->GetY()/10000,channel->GetX()/10000);
00437                         }
00438                         if(hit_time==SquareTime+5) {
00439                           hxy13->Fill(channel->GetY()/10000,channel->GetX()/10000);
00440                           hxy13->Fill(channel->GetY()/10000,channel->GetX()/10000);
00441                           hxy13->Fill(channel->GetY()/10000,channel->GetX()/10000);
00442                         }
00443                         if(hit_time==SquareTime+6) {
00444                           hxy14->Fill(channel->GetY()/10000,channel->GetX()/10000);
00445                           hxy14->Fill(channel->GetY()/10000,channel->GetX()/10000);
00446                           hxy14->Fill(channel->GetY()/10000,channel->GetX()/10000);
00447                         }
00448                       }
00449                     }
00450                     else{
00451                       hxy->Fill(channel->GetY()/10000,channel->GetX()/10000);
00452                       if(hit_time==SquareTime) hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000);
00453                       if(hit_time==SquareTime-3) hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000);
00454                       if(hit_time==SquareTime-2) hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000);
00455                       if(hit_time==SquareTime-1) hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000);
00456                       if(hit_time==SquareTime+1) hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000);
00457                       if(hit_time==SquareTime+2) hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000);
00458                       if(hit_time==SquareTime+3) hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000);
00459                       if(hit_time==SquareTime-6) hxy9->Fill(channel->GetY()/10000,channel->GetX()/10000);
00460                       if(hit_time==SquareTime-5) hxy10->Fill(channel->GetY()/10000,channel->GetX()/10000);
00461                       if(hit_time==SquareTime-4) hxy11->Fill(channel->GetY()/10000,channel->GetX()/10000);
00462                       if(hit_time==SquareTime+4) hxy12->Fill(channel->GetY()/10000,channel->GetX()/10000);
00463                       if(hit_time==SquareTime+5) hxy13->Fill(channel->GetY()/10000,channel->GetX()/10000);
00464                       if(hit_time==SquareTime+6) hxy14->Fill(channel->GetY()/10000,channel->GetX()/10000);
00465                     }
00466                   }
00467                   
00468                   //plots
00469                   c0->Divide(7,2);
00470                   c0->cd(1);          
00471                   hxy->GetXaxis()->SetTitle("X (cm)");
00472                   hxy->GetYaxis()->SetTitle("Y (cm)");
00473                   hxy->Draw("zcol");
00474                   
00475                   c0->cd(2);
00476                   hxy9->GetXaxis()->SetTitle("X (cm)");
00477                   hxy9->GetYaxis()->SetTitle("Y (cm)");
00478                   hxy9->Draw("zcol");
00479                   
00480                   c0->cd(3);
00481                   hxy10->GetXaxis()->SetTitle("X (cm)");
00482                   hxy10->GetYaxis()->SetTitle("Y (cm)");
00483                   hxy10->Draw("zcol");
00484                   
00485                   c0->cd(4);
00486                   hxy11->GetXaxis()->SetTitle("X (cm)");
00487                   hxy11->GetYaxis()->SetTitle("Y (cm)");
00488                   hxy11->Draw("zcol");
00489                   
00490                   c0->cd(5);
00491                   hxy3->GetXaxis()->SetTitle("X (cm)");
00492                   hxy3->GetYaxis()->SetTitle("Y (cm)");
00493                   hxy3->Draw("zcol");
00494                   
00495                   c0->cd(6);
00496                   hxy4->GetXaxis()->SetTitle("X (cm)");
00497                   hxy4->GetYaxis()->SetTitle("Y (cm)");
00498                   hxy4->Draw("zcol");
00499                   
00500                   c0->cd(7);
00501                   hxy5->GetXaxis()->SetTitle("X (cm)");
00502                   hxy5->GetYaxis()->SetTitle("Y (cm)");
00503                   hxy5->Draw("zcol");
00504                   
00505                   c0->cd(8);
00506                   hxy2->GetXaxis()->SetTitle("X (cm)");
00507                   hxy2->GetYaxis()->SetTitle("Y (cm)");
00508                   hxy2->Draw("zcol");
00509                   
00510                   c0->cd(9);
00511                   hxy6->GetXaxis()->SetTitle("X (cm)");
00512                   hxy6->GetYaxis()->SetTitle("Y (cm)");
00513                   hxy6->Draw("zcol");
00514                   
00515                   c0->cd(10);
00516                   hxy7->GetXaxis()->SetTitle("X (cm)");
00517                   hxy7->GetYaxis()->SetTitle("Y (cm)");
00518                   hxy7->Draw("zcol");
00519                   
00520                   c0->cd(11);
00521                   hxy8->GetXaxis()->SetTitle("X (cm)");
00522                   hxy8->GetYaxis()->SetTitle("Y (cm)");
00523                   hxy8->Draw("zcol");
00524                   
00525                   c0->cd(12);
00526                   hxy12->GetXaxis()->SetTitle("X (cm)");
00527                   hxy12->GetYaxis()->SetTitle("Y (cm)");
00528                   hxy12->Draw("zcol");
00529                   
00530                   c0->cd(13);
00531                   hxy13->GetXaxis()->SetTitle("X (cm)");
00532                   hxy13->GetYaxis()->SetTitle("Y (cm)");
00533                   hxy13->Draw("zcol");
00534                   
00535                   c0->cd(14);
00536                   hxy14->GetXaxis()->SetTitle("X (cm)");
00537                   hxy14->GetYaxis()->SetTitle("Y (cm)");
00538                   hxy14->Draw("zcol");
00539                   
00540                   //save the plot
00541                   char plotname[20] = "TH123_plot_evt";
00542                   char buf[10];
00543                   sprintf(buf,"%d",i+1);
00544                   strcat(plotname, buf);
00545                   strcat(plotname,".gif");
00546                   c0->SaveAs(plotname);
00547                   
00548                   //clear the plot
00549                   c0->Clear();
00550                   
00551                   //reset the histogramms
00552                   hxy->Reset();
00553                   hxy2->Reset();
00554                   hxy3->Reset();
00555                   hxy4->Reset();
00556                   hxy5->Reset();
00557                   hxy6->Reset();
00558                   hxy7->Reset();
00559                   hxy8->Reset();
00560                   hxy9->Reset();
00561                   hxy10->Reset();
00562                   hxy11->Reset();
00563                   hxy12->Reset();
00564                   hxy13->Reset();
00565                   hxy14->Reset();
00566                 }
00567                 if(nb_tot > nb_plots) break;
00568                 
00569               }
00570               
00571               //Now the outputs: 
00572               if(blabla==2||blabla==3) cout<<endl;
00573               if(blabla==1||blabla==2||blabla==3) {
00574                 
00575                 //at each key print nb of evts square
00576                 cout<<"In key "<<counter<<":"<<endl<<endl;;
00577                 cout<<"   Nb square evt = "<<nb_tot<<endl;
00578                 cout<<"   Nb square evt / Nevent tot = "<<1.0*nb_tot/nEvent*100.<<" %"<<endl;
00579                 //      cout<<"   Nb square evt within the timecut (delta t < "<<timecut<<") = "<<nb_timecut<<endl;
00580                 float rate = float(nb_tot)/(end_t-begin_t)*1000000000/200*60;
00581                 cout<<"   Square event rate = "<<rate<<" event/min"<<endl<<endl;
00582                 
00583                 //at each key print print the liste of square events
00584                 cout<<"List of square events: "<<endl;
00585                 for(int it=0;it<nb_tot;it++) cout<<" "<<list_event[it];
00586                 cout<<endl<<endl;
00587                 
00588                 //       //at each key print print the liste of square events within the given timecut
00589                 //       cout<<"List of square events that are square within the timecut (delta t < "<<timecut<<")"<<endl;
00590                 //       for(int it=0;it<nb_timecut;it++) cout<<" "<<list_event_timecut[it];
00591                 //       cout<<endl<<endl;
00592               }
00593               
00594               if(blabla==2||blabla==3) {
00595                 
00596                 //if you want more details print also the position(asu) and time(memory order) of the square events
00597                 cout<<"Hereafter I show the list of square events in each board/asu at each time and their number of hits:"<<endl;
00598                 cout<<"Event, Board/asu, time, Nb_hits"<<endl;
00599                 for(int it=0;it<nb_tot;it++) {
00600                   for(int itt=0;itt<4;itt++) cout<<list_square_events[it][itt]<<";  ";
00601                   cout<<endl;
00602                 }
00603                 cout<<endl;
00604                 //      cout<<"Hereafter I show the list of square events that are square within the timecut (delta t < "<<timecut<<"):"<<endl;
00605                 //      cout<<"Event, Board/asu, time, Nb_hits"<<endl;
00606                 //      for(int it=0;it<nb_timecut;it++) {
00607                 //      for(int itt=0;itt<4;itt++) cout<<list_square_events_timecut[it][itt]<<";  ";
00608                 //      cout<<endl;
00609                 //      }
00610                 cout<<endl;
00611               }
00612               cout<<endl;
00613             }
00614         }
00615     }
00616   theApp->Run();
00617   delete theApp;
00618   return 0; 
00619 }
00620
/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/root/square_event_2D_2012.cpp
