#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_2D.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_2D.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_2D.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 = 400; 00066 00067 //energy threshold 00068 Int_t E_threshold = 123; 00069 00070 //nb plots 00071 Int_t nb_plots = 30; 00072 00073 //timecut (time interval after trigger for tagging an event as square and remove it from efficency calculation) 00074 // Int_t timecut = 50; 00075 00076 //input data 00077 //string file = "04082011_2236_MGT.root"; 00078 //string file = "05082011_1055_MGT.root"; 00079 //string file = "05082011_2248_MGT.root"; 00080 // 00081 //string file = "06082011_0002_MGT.root"; 00082 //string file = "06082011_0034_MGT.root"; 00083 //string file = "06082011_0105_MGT.root"; 00084 //string file = "06082011_0124_MGT.root"; 00085 //string file = "06082011_0152_MGT.root"; 00086 //string file = "06082011_0213_MGT.root"; 00087 //string file = "06082011_0231_MGT.root"; 00088 //string file = "06082011_0514_MGT.root"; 00089 //string file = "06082011_0536_MGT.root"; 00090 //string file = "06082011_0554_MGT.root"; 00091 //string file = "06082011_0652_MGT.root"; 00092 //string file = "06082011_0727_MGT.root"; 00093 //string file = "06082011_0829_MGT.root"; 00094 string file = "/lapp_data/LC/Detecteurs/MicroMegas/data/TB2012/SPS_H2_may_2012/Production/CaloHits/CaloHit_DHCAL_714448_I4_0.slcio.root"; 00095 //----------------- 00096 00097 00098 //To see the plots in "live" 00099 TApplication *theApp; 00100 theApp = new TRint("App", 0, NULL); 00101 TCanvas *c0 = new TCanvas("c0", "Prototype m2", 100, 50, 1000, 600); 00102 TH2I *hxy = new TH2I("hxy","Hit position distribution (No time cut);y (cm);x (cm)",96,0,96,96,0,96); 00103 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); 00104 TH2I *hxy3 = new TH2I("hxy3","Hit position distribution (at t_square - 3 bcid);y (cm);x (cm)",96,0,96,96,0,96); 00105 TH2I *hxy4 = new TH2I("hxy4","Hit position distribution (at t_square - 2 bcid);y (cm);x (cm)",96,0,96,96,0,96); 00106 TH2I *hxy5 = new TH2I("hxy5","Hit position distribution (at t_square - 1 bcid);y (cm);x (cm)",96,0,96,96,0,96); 00107 TH2I *hxy6 = new TH2I("hxy6","Hit position distribution (at t_square + 1 bcid);y (cm);x (cm)",96,0,96,96,0,96); 00108 TH2I *hxy7 = new TH2I("hxy7","Hit position distribution (at t_square + 2 bcid);y (cm);x (cm)",96,0,96,96,0,96); 00109 TH2I *hxy8 = new TH2I("hxy8","Hit position distribution (at t_square + 3 bcid);y (cm);x (cm)",96,0,96,96,0,96); 00110 00111 00112 //our data file 00113 TFile *f = new TFile(file.c_str()); 00114 if(blabla==1||blabla==2||blabla==3){ 00115 cout<<endl; 00116 cout<<"reading File : "<<file<<endl; 00117 00118 } 00119 00120 //key iterator 00121 TIter nextkey(f->GetListOfKeys()); 00122 TKey *key; 00123 UInt_t counter = 0; 00124 00125 if(blabla==1||blabla==2||blabla==3) cout<<endl; 00126 00127 00128 //loop on all the keys 00129 while (key = (TKey*)nextkey()) { 00130 00131 counter++; 00132 if(blabla==1||blabla==2||blabla==3) cout<<"key number = "<<counter<<endl; 00133 00134 TTree *tree = (TTree*)key->ReadObj(); 00135 MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun"); 00136 00137 Int_t nEvent = tree->GetEntries(); 00138 if(blabla==1||blabla==2||blabla==3) cout <<"Nevent = "<<nEvent<<endl; 00139 00140 MTEvent *evt = new MTEvent(); 00141 00142 TBranch *branch= tree->GetBranch("MTEvent"); 00143 branch->SetAddress(&evt); 00144 if(blabla==2||blabla==3) f->Print(); 00145 00146 //at each new key, declare the event counters and indicators 00147 bool IsSquare = 0; 00148 Int_t list_square_events[nEvent][4]; 00149 // Int_t list_square_events_timecut[nEvent][4]; 00150 UInt_t list_event[nEvent]; 00151 // UInt_t list_event_timecut[nEvent]; 00152 UInt_t nb_tot = 0; 00153 // UInt_t nb_timecut = 0; 00154 UInt_t begin_t = 0; 00155 UInt_t end_t = 0; 00156 00157 00158 //loop on all the events of the key 00159 for (int i=0;i<nEvent-1;i++){ 00160 00161 tree->GetEntry(i); 00162 00163 UInt_t nchannel = evt->GetNchannel(); 00164 //if(blabla==3) cout<<"Event "<<i+1<<" / "<<nEvent<<"; nb of channels= "<<nchannel<<endl; 00165 00166 //at each event declare the hits counter 00167 map<unsigned int,map<unsigned int,int> > channels_map;// key1=asu (board_id), key2=hit_time 00168 00169 //at each event declare the time interval 00170 Int_t tmin=0; 00171 Int_t tmax=0; 00172 00173 //at each event reset the square indicator 00174 IsSquare=0; 00175 00176 //at each event declare the time where the event is square 00177 UInt_t SquareTime=0; 00178 00179 //loop on all the channels of the event 00180 for (int j=0;j<nchannel;j++){ 00181 00182 //this channel 00183 MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j); 00184 00185 //at each channel get board_id and hit_time 00186 UInt_t hit_time = channel->GetBcIdAbs() - (channel->GetBcIdDif() - channel->GetBcIdHit() ); 00187 UInt_t board_id = channel->GetBoardId(); 00188 int DigiEng = channel->GetDigitalValue(); 00189 //cout<<"Channel : "<<j<<"; hit on asu "<<board_id<<"; at time "<<hit_time<<"; with DigiEng "<<DigiEng<<endl; 00190 00191 //get tmin and tmax 00192 if(j==0) { 00193 tmin=hit_time; 00194 tmax=hit_time; 00195 } 00196 if(hit_time<tmin) tmin=hit_time; 00197 if(hit_time>tmax) tmax=hit_time; 00198 00199 //get begin_t and end_t 00200 if(i==0 && j==0) { 00201 begin_t=hit_time; 00202 end_t=hit_time; 00203 } 00204 if(hit_time<begin_t) begin_t=hit_time; 00205 if(hit_time>end_t) end_t=hit_time; 00206 00207 //upgrade the hits counter of the good board at the good time 00208 if(E_threshold==123) if(DigiEng>=1) add(channels_map, board_id, hit_time); 00209 else if(DigiEng>=E_threshold) add(channels_map, board_id, hit_time); 00210 } 00211 00212 unsigned int timeold=0; 00213 00214 //Loop on each board, at each time, to identify if this event is square 00215 map<unsigned int, map<unsigned int,int> >::iterator iter1; 00216 for(iter1=channels_map.begin();iter1!=channels_map.end();iter1++){ 00217 00218 unsigned int azu=iter1->first; 00219 map<unsigned int,int> hits_per_time=iter1->second; 00220 00221 map<unsigned int,int>::iterator iter2; 00222 for(iter2=hits_per_time.begin();iter2!=hits_per_time.end();iter2++){ 00223 00224 unsigned int time = iter2->first; 00225 int nb_hits = iter2->second; 00226 00227 //if this is a square event 00228 if (nb_hits>=threshold) { 00229 00230 //set square event indicator to 1 00231 IsSquare=1; 00232 00233 //set the time of this event 00234 SquareTime=time; 00235 00236 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; 00237 timeold=time; 00238 00239 //update the event square lists :event, asu, time, nb_hits 00240 list_square_events[nb_tot][0] = i+1; 00241 list_square_events[nb_tot][1] = azu; 00242 list_square_events[nb_tot][2] = time; 00243 list_square_events[nb_tot][3] = nb_hits; 00244 // if(time>tmax-timecut){ 00245 // list_square_events_timecut[nb_timecut][0] = i+1; 00246 // list_square_events_timecut[nb_timecut][1] = azu; 00247 // list_square_events_timecut[nb_timecut][2] = time; 00248 // list_square_events_timecut[nb_timecut][3] = nb_hits; 00249 // } 00250 00251 //update the basic event square counters 00252 nb_tot++; 00253 // if(time>tmax-timecut) nb_timecut++; 00254 00255 //update the basic square event lists 00256 list_event[nb_tot-1] = i+1; 00257 // if(time>tmax-timecut) list_event_timecut[nb_tot-1] = i+1; 00258 } 00259 } 00260 } 00261 00262 //if we want to plot this event 00263 if (IsSquare && plot){ 00264 00265 //loop on all the channels of the event 00266 for (int j=0;j<nchannel;j++){ 00267 00268 //this channel 00269 MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j); 00270 00271 //at each channel get board_id and hit_time 00272 UInt_t board_id = channel->GetBoardId(); 00273 UInt_t hit_time = channel->GetBcIdAbs() - (channel->GetBcIdDif() - channel->GetBcIdHit() ); 00274 int DigiEng = channel->GetDigitalValue(); 00275 00276 //fill the histogramms 00277 if(E_threshold==123){ 00278 if(DigiEng>=1){ 00279 hxy->Fill(channel->GetY()/10000,channel->GetX()/10000); 00280 if(hit_time==SquareTime) hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000); 00281 if(hit_time==SquareTime-3) hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000); 00282 if(hit_time==SquareTime-2) hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000); 00283 if(hit_time==SquareTime-1) hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000); 00284 if(hit_time==SquareTime+1) hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000); 00285 if(hit_time==SquareTime+2) hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000); 00286 if(hit_time==SquareTime+3) hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000); 00287 } 00288 if(DigiEng>=2){ 00289 hxy->Fill(channel->GetY()/10000,channel->GetX()/10000); 00290 hxy->Fill(channel->GetY()/10000,channel->GetX()/10000); 00291 if(hit_time==SquareTime) { 00292 //cout <<"2 "; 00293 hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000); 00294 hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000); 00295 } 00296 if(hit_time==SquareTime-3) { 00297 hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000); 00298 hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000); 00299 } 00300 if(hit_time==SquareTime-2) { 00301 hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000); 00302 hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000); 00303 } 00304 if(hit_time==SquareTime-1) { 00305 hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000); 00306 hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000); 00307 } 00308 if(hit_time==SquareTime+1) { 00309 hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000); 00310 hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000); 00311 } 00312 if(hit_time==SquareTime+2) { 00313 hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000); 00314 hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000); 00315 } 00316 if(hit_time==SquareTime+3) { 00317 hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000); 00318 hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000); 00319 } 00320 } 00321 if(DigiEng>=3){ 00322 hxy->Fill(channel->GetY()/10000,channel->GetX()/10000); 00323 hxy->Fill(channel->GetY()/10000,channel->GetX()/10000); 00324 hxy->Fill(channel->GetY()/10000,channel->GetX()/10000); 00325 if(hit_time==SquareTime) { 00326 //cout <<"3 "; 00327 hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000); 00328 hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000); 00329 hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000); 00330 } 00331 if(hit_time==SquareTime-3) { 00332 hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000); 00333 hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000); 00334 hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000); 00335 } 00336 if(hit_time==SquareTime-2) { 00337 hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000); 00338 hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000); 00339 hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000); 00340 } 00341 if(hit_time==SquareTime-1) { 00342 hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000); 00343 hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000); 00344 hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000); 00345 } 00346 if(hit_time==SquareTime+1) { 00347 hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000); 00348 hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000); 00349 hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000); 00350 } 00351 if(hit_time==SquareTime+2) { 00352 hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000); 00353 hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000); 00354 hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000); 00355 } 00356 if(hit_time==SquareTime+3) { 00357 hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000); 00358 hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000); 00359 hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000); 00360 } 00361 } 00362 } 00363 else{ 00364 hxy->Fill(channel->GetY()/10000,channel->GetX()/10000); 00365 if(hit_time==SquareTime) hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000); 00366 if(hit_time==SquareTime-3) hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000); 00367 if(hit_time==SquareTime-2) hxy4->Fill(channel->GetY()/10000,channel->GetX()/10000); 00368 if(hit_time==SquareTime-1) hxy5->Fill(channel->GetY()/10000,channel->GetX()/10000); 00369 if(hit_time==SquareTime+1) hxy6->Fill(channel->GetY()/10000,channel->GetX()/10000); 00370 if(hit_time==SquareTime+2) hxy7->Fill(channel->GetY()/10000,channel->GetX()/10000); 00371 if(hit_time==SquareTime+3) hxy8->Fill(channel->GetY()/10000,channel->GetX()/10000); 00372 } 00373 } 00374 00375 //plots 00376 c0->Divide(4,2); 00377 c0->cd(1); 00378 hxy->GetXaxis()->SetTitle("X (cm)"); 00379 hxy->GetYaxis()->SetTitle("Y (cm)"); 00380 hxy->Draw("zcol"); 00381 00382 c0->cd(2); 00383 hxy3->GetXaxis()->SetTitle("X (cm)"); 00384 hxy3->GetYaxis()->SetTitle("Y (cm)"); 00385 hxy3->Draw("zcol"); 00386 00387 c0->cd(3); 00388 hxy4->GetXaxis()->SetTitle("X (cm)"); 00389 hxy4->GetYaxis()->SetTitle("Y (cm)"); 00390 hxy4->Draw("zcol"); 00391 00392 c0->cd(4); 00393 hxy5->GetXaxis()->SetTitle("X (cm)"); 00394 hxy5->GetYaxis()->SetTitle("Y (cm)"); 00395 hxy5->Draw("zcol"); 00396 00397 c0->cd(5); 00398 hxy2->GetXaxis()->SetTitle("X (cm)"); 00399 hxy2->GetYaxis()->SetTitle("Y (cm)"); 00400 hxy2->Draw("zcol"); 00401 00402 c0->cd(6); 00403 hxy6->GetXaxis()->SetTitle("X (cm)"); 00404 hxy6->GetYaxis()->SetTitle("Y (cm)"); 00405 hxy6->Draw("zcol"); 00406 00407 c0->cd(7); 00408 hxy7->GetXaxis()->SetTitle("X (cm)"); 00409 hxy7->GetYaxis()->SetTitle("Y (cm)"); 00410 hxy7->Draw("zcol"); 00411 00412 c0->cd(8); 00413 hxy8->GetXaxis()->SetTitle("X (cm)"); 00414 hxy8->GetYaxis()->SetTitle("Y (cm)"); 00415 hxy8->Draw("zcol"); 00416 00417 //save the plot 00418 char plotname[20] = "TH123_plot_evt"; 00419 char buf[10]; 00420 sprintf(buf,"%d",i+1); 00421 strcat(plotname, buf); 00422 strcat(plotname,".gif"); 00423 c0->SaveAs(plotname); 00424 00425 //clear the plot 00426 c0->Clear(); 00427 00428 //reset the histogramms 00429 hxy->Reset(); 00430 hxy2->Reset(); 00431 hxy3->Reset(); 00432 hxy4->Reset(); 00433 hxy5->Reset(); 00434 hxy6->Reset(); 00435 hxy7->Reset(); 00436 hxy8->Reset(); 00437 } 00438 if(nb_tot > nb_plots) break; 00439 } 00440 00441 //Now the outputs: 00442 if(blabla==2||blabla==3) cout<<endl; 00443 if(blabla==1||blabla==2||blabla==3) { 00444 00445 //at each key print nb of evts square 00446 cout<<"In key "<<counter<<":"<<endl<<endl;; 00447 cout<<" Nb square evt = "<<nb_tot<<endl; 00448 cout<<" Nb square evt / Nevent tot = "<<1.0*nb_tot/nEvent*100.<<" %"<<endl; 00449 // cout<<" Nb square evt within the timecut (delta t < "<<timecut<<") = "<<nb_timecut<<endl; 00450 float rate = float(nb_tot)/(end_t-begin_t)*1000000000/200*60; 00451 cout<<" Square event rate = "<<rate<<" event/min"<<endl<<endl; 00452 00453 //at each key print print the liste of square events 00454 cout<<"List of square events: "<<endl; 00455 for(int it=0;it<nb_tot;it++) cout<<" "<<list_event[it]; 00456 cout<<endl<<endl; 00457 00458 // //at each key print print the liste of square events within the given timecut 00459 // cout<<"List of square events that are square within the timecut (delta t < "<<timecut<<")"<<endl; 00460 // for(int it=0;it<nb_timecut;it++) cout<<" "<<list_event_timecut[it]; 00461 // cout<<endl<<endl; 00462 } 00463 00464 if(blabla==2||blabla==3) { 00465 00466 //if you want more details print also the position(asu) and time(memory order) of the square events 00467 cout<<"Hereafter I show the list of square events in each board/asu at each time and their number of hits:"<<endl; 00468 cout<<"Event, Board/asu, time, Nb_hits"<<endl; 00469 for(int it=0;it<nb_tot;it++) { 00470 for(int itt=0;itt<4;itt++) cout<<list_square_events[it][itt]<<"; "; 00471 cout<<endl; 00472 } 00473 cout<<endl; 00474 // cout<<"Hereafter I show the list of square events that are square within the timecut (delta t < "<<timecut<<"):"<<endl; 00475 // cout<<"Event, Board/asu, time, Nb_hits"<<endl; 00476 // for(int it=0;it<nb_timecut;it++) { 00477 // for(int itt=0;itt<4;itt++) cout<<list_square_events_timecut[it][itt]<<"; "; 00478 // cout<<endl; 00479 // } 00480 cout<<endl; 00481 } 00482 cout<<endl; 00483 } 00484 theApp->Run(); 00485 delete theApp; 00486 return 0; 00487 }
1.4.7