#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"Include dependency graph for square_event.cpp:
Go to the source code of this file.
Functions | |
| void | add (map< unsigned int, map< unsigned int, int > > &channels_map, unsigned int asu, unsigned int hit_time) |
| int | main (void) |
| void add | ( | map< unsigned int, map< unsigned int, int > > & | channels_map, | |
| unsigned int | asu, | |||
| unsigned int | hit_time | |||
| ) |
Definition at line 29 of file square_event.cpp.
00029 { 00030 00031 if ( channels_map.find(asu) == channels_map.end()) { // if this asu doesn't exist yet 00032 00033 map<unsigned int, int> nbHitMap; 00034 nbHitMap[hit_time]=1; 00035 channels_map.insert(make_pair(asu,nbHitMap)); 00036 } 00037 else{ // if this asu already exist 00038 00039 map<unsigned int, int> &nbHitMap = channels_map.find(asu)->second; 00040 if (nbHitMap.find(hit_time) == nbHitMap.end()){ //if this hit_time doesn't exist yet 00041 00042 nbHitMap.insert(make_pair(hit_time,1)); 00043 } 00044 else{ //if this abs time already exist 00045 00046 nbHitMap[hit_time] = nbHitMap[hit_time]++; 00047 } 00048 } 00049 }
| int main | ( | void | ) |
Definition at line 53 of file square_event.cpp.
00053 { 00054 00055 //----------------- 00056 //User's parameters 00057 00058 //do you want to plot (0 = no plot, 1 = plots) 00059 bool plot = 1; 00060 00061 //Do you want some info ? (0 = not at all, 1 = yes a little, 2 = yes a lot, 3 = yes all : for debug purpose) 00062 Int_t blabla = 2; 00063 00064 //threshold (number of hits in a given asu at a given time (maximum=1536)) 00065 Int_t threshold = 1000; 00066 00067 //timecut (time interval after trigger for tagging an event as square and remove it from efficency calculation) 00068 Int_t timecut = 50; 00069 00070 //input data 00071 //string file = "04082011_2236_MGT.root"; 00072 //string file = "05082011_1055_MGT.root"; 00073 string file = "05082011_2248_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 square);y (cm);x (cm)",96,0,96,96,0,96); 00084 //TH1I *hdt = new TH1I("hdt","Hit time distribution (No time cut);N_hits;Time",10003,-1,10001); 00085 //TH1I *hdt2 = new TH1I("hdt2","Hit time distribution (within time cut);N_hits;Time",103,-1,101); 00086 //TH1I *hdt3 = new TH1I("hdt3","Hit time distribution (around the time where it is square);N_hits;Time",23,-1,21); 00087 00088 00089 //our data file 00090 TFile *f = new TFile(file.c_str()); 00091 if(blabla==1||blabla==2||blabla==3){ 00092 cout<<endl; 00093 cout<<"reading File : "<<file<<endl; 00094 } 00095 00096 //key iterator 00097 TIter nextkey(f->GetListOfKeys()); 00098 TKey *key; 00099 UInt_t counter = 0; 00100 00101 if(blabla==1||blabla==2||blabla==3) cout<<endl; 00102 00103 00104 //loop on all the keys 00105 while (key = (TKey*)nextkey()) { 00106 00107 counter++; 00108 if(blabla==1||blabla==2||blabla==3) cout<<"key number = "<<counter<<endl; 00109 00110 TTree *tree = (TTree*)key->ReadObj(); 00111 MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun"); 00112 00113 Int_t nEvent = tree->GetEntries(); 00114 if(blabla==1||blabla==2||blabla==3) cout <<"Nevent = "<<nEvent<<endl; 00115 00116 MTEvent *evt = new MTEvent(); 00117 00118 TBranch *branch= tree->GetBranch("MTEvent"); 00119 branch->SetAddress(&evt); 00120 00121 //at each new key, declare the event counters and indicators 00122 bool IsSquare = 0; 00123 Int_t list_square_events[nEvent][4]; 00124 Int_t list_square_events_timecut[nEvent][4]; 00125 UInt_t list_event[nEvent]; 00126 UInt_t list_event_timecut[nEvent]; 00127 UInt_t nb_tot = 0; 00128 UInt_t nb_timecut = 0; 00129 00130 00131 //loop on all the events of the key 00132 for (int i=0;i<nEvent-1;i++){ 00133 00134 tree->GetEntry(i); 00135 00136 UInt_t nchannel = evt->GetNchannel(); 00137 //if(blabla==3) cout<<"Event "<<i+1<<" / "<<nEvent<<"; nb of channels= "<<nchannel<<endl; 00138 00139 //at each event declare the hits counter 00140 map<unsigned int,map<unsigned int,int> > channels_map;// key1=asu (board_id), key2=hit_time 00141 00142 //at each event declare the time interval 00143 Int_t tmin=0; 00144 Int_t tmax=0; 00145 00146 //at each event reset the square indicator 00147 IsSquare=0; 00148 00149 //at each event declare the time where the event is square 00150 UInt_t SquareTime=0; 00151 00152 //loop on all the channels of the event 00153 for (int j=0;j<nchannel;j++){ 00154 00155 //this channel 00156 MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j); 00157 00158 //at each channel get board_id and hit_time 00159 UInt_t hit_time = channel->GetBcIdAbs() - (channel->GetBcIdDif() - channel->GetBcIdHit() ); 00160 UInt_t board_id = channel->GetBoardId(); 00161 //if(blabla==3) cout<<"Channel : "<<j<<"; hit on asu "<<board_id<<"; at time "<<hit_time<<endl; 00162 00163 //get tmin and tmax 00164 if(j==0) { 00165 tmin=hit_time; 00166 tmax=hit_time; 00167 } 00168 if(hit_time<tmin) tmin=hit_time; 00169 if(hit_time>tmax) tmax=hit_time; 00170 00171 //upgrade the hits counter of the good board at the good time 00172 add(channels_map, board_id, hit_time); 00173 } 00174 00175 unsigned int timeold=0; 00176 00177 //Loop on each board, at each time, to identify if this event is square 00178 map<unsigned int, map<unsigned int,int> >::iterator iter1; 00179 for(iter1=channels_map.begin();iter1!=channels_map.end();iter1++){ 00180 00181 unsigned int azu=iter1->first; 00182 map<unsigned int,int> hits_per_time=iter1->second; 00183 00184 map<unsigned int,int>::iterator iter2; 00185 for(iter2=hits_per_time.begin();iter2!=hits_per_time.end();iter2++){ 00186 00187 unsigned int time = iter2->first; 00188 int nb_hits = iter2->second; 00189 00190 //if this is a square event 00191 if (nb_hits>=threshold) { 00192 00193 //set square event indicator to 1 00194 //this.set->IsSquare=true 00195 IsSquare=1; 00196 00197 //set the time of this event 00198 SquareTime=time; 00199 00200 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; 00201 timeold=time; 00202 00203 //update the event square lists :event, asu, time, nb_hits 00204 list_square_events[nb_tot][0] = i+1; 00205 list_square_events[nb_tot][1] = azu; 00206 list_square_events[nb_tot][2] = time; 00207 list_square_events[nb_tot][3] = nb_hits; 00208 if(time>tmax-timecut){ 00209 list_square_events_timecut[nb_timecut][0] = i+1; 00210 list_square_events_timecut[nb_timecut][1] = azu; 00211 list_square_events_timecut[nb_timecut][2] = time; 00212 list_square_events_timecut[nb_timecut][3] = nb_hits; 00213 } 00214 00215 //update the basic event square counters 00216 nb_tot++; 00217 if(time>tmax-timecut) nb_timecut++; 00218 00219 //update the basic square event lists 00220 list_event[nb_tot-1] = i+1; 00221 if(time>tmax-timecut) list_event_timecut[nb_tot-1] = i+1; 00222 } 00223 } 00224 } 00225 00226 //if we want to plot this event 00227 if (IsSquare && plot){ 00228 00229 TH1I *hdt = new TH1I("hdt","Hit time distribution (No time cut);N_hits;Time",tmax-tmin+3,-1,tmax-tmin+1); 00230 TH1I *hdt2 = new TH1I("hdt2","Hit time distribution (within time cut);N_hits;Time",timecut+3,tmax-timecut-1,tmax+1); 00231 TH1I *hdt3 = new TH1I("hdt3","Hit time distribution (around the time where it is square);N_hits;Time",23,SquareTime-tmin-11,SquareTime-tmin+11); 00232 00233 //loop on all the channels of the event 00234 for (int j=0;j<nchannel;j++){ 00235 00236 //this channel 00237 MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j); 00238 00239 //at each channel get board_id and hit_time 00240 UInt_t board_id = channel->GetBoardId(); 00241 UInt_t hit_time = channel->GetBcIdAbs() - (channel->GetBcIdDif() - channel->GetBcIdHit() ); 00242 00243 //fill the histogramms 00244 hxy->Fill(channel->GetY()/10000,channel->GetX()/10000); 00245 if(hit_time>tmax-timecut) hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000); 00246 if(hit_time==SquareTime) hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000); 00247 hdt->Fill(hit_time-tmin); 00248 if(hit_time>tmax-timecut) hdt2->Fill(hit_time-tmin); 00249 if(hit_time>SquareTime-100 && hit_time<SquareTime+500) hdt3->Fill(hit_time-tmin); 00250 } 00251 00252 //plots 00253 c0->Divide(3,2); 00254 c0->cd(1); 00255 hxy->GetXaxis()->SetTitle("X (cm)"); 00256 hxy->GetYaxis()->SetTitle("Y (cm)"); 00257 hxy->Draw("zcol"); 00258 00259 c0->cd(2); 00260 hxy2->GetXaxis()->SetTitle("X (cm)"); 00261 hxy2->GetYaxis()->SetTitle("Y (cm)"); 00262 hxy2->Draw("zcol"); 00263 00264 c0->cd(3); 00265 hxy3->GetXaxis()->SetTitle("X (cm)"); 00266 hxy3->GetYaxis()->SetTitle("Y (cm)"); 00267 hxy3->Draw("zcol"); 00268 00269 c0->cd(4); 00270 //hdt->SetAxisRange(-1,tmax-tmin+1); 00271 //hdt->SetMinimum(-1); 00272 //hdt->SetMaximum(tmax-tmin+1); 00273 hdt->GetXaxis()->SetTitle("time after trigger (BcId)"); 00274 hdt->GetYaxis()->SetTitle("Nb hits"); 00275 hdt->Draw(); 00276 00277 c0->cd(5); 00278 //hdt2->SetAxisRange(-1,timecut+1); 00279 //hdt2->SetMinimum(-1); 00280 //hdt2->SetMaximum(timecut+1); 00281 hdt2->GetXaxis()->SetTitle("time after trigger (BcId)"); 00282 hdt2->GetYaxis()->SetTitle("Nb hits"); 00283 hdt2->Draw(); 00284 00285 c0->cd(6); 00286 //hdt3->SetAxisRange(SquareTime-tmin-11,SquareTime-tmin+11); 00287 //hdt3->SetMinimum(SquareTime-tmin-11); 00288 //hdt3->SetMaximum(SquareTime-tmin+11); 00289 hdt3->GetXaxis()->SetTitle("time after trigger (BcId)"); 00290 hdt3->GetYaxis()->SetTitle("Nb hits"); 00291 hdt3->Draw(); 00292 00293 //save the plot 00294 char plotname[20] = "plot_evt"; 00295 char buf[10]; 00296 sprintf(buf,"%d",i+1); 00297 strcat(plotname, buf); 00298 strcat(plotname,".gif"); 00299 c0->SaveAs(plotname); 00300 00301 //clear the plot 00302 c0->Clear(); 00303 00304 //reset the histogramms 00305 hxy->Reset(); 00306 hxy2->Reset(); 00307 hxy3->Reset(); 00308 //hdt->Reset(); 00309 //hdt2->Reset(); 00310 //hdt3->Reset(); 00311 //hxy->Delete(); 00312 //hxy2->Delete(); 00313 //hxy3->Delete(); 00314 hdt->Delete(); 00315 hdt2->Delete(); 00316 hdt3->Delete(); 00317 } 00318 } 00319 00320 //Now the outputs: 00321 if(blabla==2||blabla==3) cout<<endl; 00322 if(blabla==1||blabla==2||blabla==3) { 00323 00324 //at each key print nb of evts square 00325 cout<<"In key "<<counter<<":"<<endl<<endl;; 00326 cout<<" Nb square evt= "<<nb_tot<<endl; 00327 cout<<" Nb square evt within the timecut (delta t < "<<timecut<<") = "<<nb_timecut<<endl<<endl; 00328 00329 //at each key print print the liste of square events 00330 cout<<"List of square events: "<<endl; 00331 for(int it=0;it<nb_tot;it++) cout<<" "<<list_event[it]; 00332 cout<<endl<<endl; 00333 00334 //at each key print print the liste of square events within the given timecut 00335 cout<<"List of square events that are square within the timecut (delta t < "<<timecut<<")"<<endl; 00336 for(int it=0;it<nb_timecut;it++) cout<<" "<<list_event_timecut[it]; 00337 cout<<endl<<endl; 00338 } 00339 00340 if(blabla==2||blabla==3) { 00341 00342 //if you want more details print also the position(asu) and time(memory order) of the square events 00343 cout<<"Hereafter I show the list of square events in each board/asu at each time and their number of hits:"<<endl; 00344 cout<<"Event, Board/asu, time, Nb_hits"<<endl; 00345 for(int it=0;it<nb_tot;it++) { 00346 for(int itt=0;itt<4;itt++) cout<<list_square_events[it][itt]<<"; "; 00347 cout<<endl; 00348 } 00349 cout<<endl; 00350 cout<<"Hereafter I show the list of square events that are square within the timecut (delta t < "<<timecut<<"):"<<endl; 00351 cout<<"Event, Board/asu, time, Nb_hits"<<endl; 00352 for(int it=0;it<nb_timecut;it++) { 00353 for(int itt=0;itt<4;itt++) cout<<list_square_events_timecut[it][itt]<<"; "; 00354 cout<<endl; 00355 } 00356 cout<<endl; 00357 } 00358 cout<<endl; 00359 } 00360 theApp->Run(); 00361 delete theApp; 00362 return 0; 00363 }
1.4.7