#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 chip_event.cpp:
Go to the source code of this file.
Functions | |
| void | add (map< unsigned int, map< unsigned int, int > > &channels_map, unsigned int chip, unsigned int hit_time) |
| int | main (void) |
| void add | ( | map< unsigned int, map< unsigned int, int > > & | channels_map, | |
| unsigned int | chip, | |||
| unsigned int | hit_time | |||
| ) |
Definition at line 29 of file chip_event.cpp.
00029 { 00030 00031 if ( channels_map.find(chip) == channels_map.end()) { // if this chip doesn't exist yet 00032 00033 map<unsigned int, int> nbHitMap; 00034 nbHitMap[hit_time]=1; 00035 channels_map.insert(make_pair(chip,nbHitMap)); 00036 } 00037 else{ // if this chip already exist 00038 00039 map<unsigned int, int> &nbHitMap = channels_map.find(chip)->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 chip_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 chip at a given time (maximum=64)) 00065 Int_t threshold = 40; 00066 00067 //timecut (time interval after trigger for tagging an event as chip 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 //string file = "06082011_0002_MGT.root"; 00076 //string file = "06082011_0034_MGT.root"; 00077 //string file = "06082011_0105_MGT.root"; 00078 //string file = "06082011_0124_MGT.root"; 00079 //string file = "06082011_0152_MGT.root"; 00080 //string file = "06082011_0213_MGT.root"; 00081 //string file = "06082011_0231_MGT.root"; 00082 //string file = "06082011_0514_MGT.root"; 00083 //string file = "06082011_0536_MGT.root"; 00084 //string file = "06082011_0554_MGT.root"; 00085 //string file = "06082011_0652_MGT.root"; 00086 //string file = "06082011_0727_MGT.root"; 00087 string file = "06082011_0829_MGT.root"; 00088 //----------------- 00089 00090 00091 //To see the plots in "live" 00092 TApplication *theApp; 00093 theApp = new TRint("App", 0, NULL); 00094 TCanvas *c0 = new TCanvas("c0", "Prototype m2", 100, 50, 1000, 600); 00095 TH2I *hxy = new TH2I("hxy","Hit position distribution (No time cut);y (cm);x (cm)",96,0,96,96,0,96); 00096 TH2I *hxy2 = new TH2I("hxy2","Hit position distribution (within time cut);y (cm);x (cm)",96,0,96,96,0,96); 00097 TH2I *hxy3 = new TH2I("hxy3","Hit position distribution (at the time where it is chip);y (cm);x (cm)",96,0,96,96,0,96); 00098 //TH1I *hdt = new TH1I("hdt","Hit time distribution (No time cut);N_hits;Time",10003,-1,10001); 00099 //TH1I *hdt2 = new TH1I("hdt2","Hit time distribution (within time cut);N_hits;Time",103,-1,101); 00100 //TH1I *hdt3 = new TH1I("hdt3","Hit time distribution (around the time where it is chip);N_hits;Time",23,-1,21); 00101 00102 00103 //our data file 00104 TFile *f = new TFile(file.c_str()); 00105 if(blabla==1||blabla==2||blabla==3){ 00106 cout<<endl; 00107 cout<<"reading File : "<<file<<endl; 00108 } 00109 00110 //key iterator 00111 TIter nextkey(f->GetListOfKeys()); 00112 TKey *key; 00113 UInt_t counter = 0; 00114 00115 if(blabla==1||blabla==2||blabla==3) cout<<endl; 00116 00117 00118 //loop on all the keys 00119 while (key = (TKey*)nextkey()) { 00120 00121 counter++; 00122 if(blabla==1||blabla==2||blabla==3) cout<<"key number = "<<counter<<endl; 00123 00124 TTree *tree = (TTree*)key->ReadObj(); 00125 MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun"); 00126 00127 Int_t nEvent = tree->GetEntries(); 00128 if(blabla==1||blabla==2||blabla==3) cout <<"Nevent = "<<nEvent<<endl; 00129 00130 MTEvent *evt = new MTEvent(); 00131 00132 TBranch *branch= tree->GetBranch("MTEvent"); 00133 branch->SetAddress(&evt); 00134 00135 //at each new key, declare the event counters and indicators 00136 bool IsChip = 0; 00137 Int_t list_chip_events[nEvent][4]; 00138 Int_t list_chip_events_timecut[nEvent][4]; 00139 UInt_t list_event[nEvent]; 00140 UInt_t list_event_timecut[nEvent]; 00141 UInt_t nb_tot = 0; 00142 UInt_t nb_timecut = 0; 00143 UInt_t begin_t = 0; 00144 UInt_t end_t = 0; 00145 00146 00147 //loop on all the events of the key 00148 for (int i=0;i<nEvent-1;i++){ 00149 00150 tree->GetEntry(i); 00151 00152 UInt_t nchannel = evt->GetNchannel(); 00153 //if(blabla==3) cout<<"Event "<<i+1<<" / "<<nEvent<<"; nb of channels= "<<nchannel<<endl; 00154 00155 //at each event declare the hits counter 00156 map<unsigned int,map<unsigned int,int> > channels_map;// key1=chip (chip_id), key2=hit_time 00157 00158 //at each event declare the time interval 00159 Int_t tmin=0; 00160 Int_t tmax=0; 00161 00162 //at each event reset the chip indicator 00163 IsChip=0; 00164 00165 //at each event declare the time where the event is chip 00166 UInt_t ChipTime=0; 00167 00168 //loop on all the channels of the event 00169 for (int j=0;j<nchannel;j++){ 00170 00171 //this channel 00172 MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j); 00173 00174 //at each channel get chip_id and hit_time 00175 UInt_t hit_time = channel->GetBcIdAbs() - (channel->GetBcIdDif() - channel->GetBcIdHit() ); 00176 UInt_t chip_id = channel->GetChipId(); 00177 //if(blabla==3) cout<<"Channel : "<<j<<"; hit on chip "<<chip_id<<"; at time "<<hit_time<<endl; 00178 00179 //get tmin and tmax 00180 if(j==0) { 00181 tmin=hit_time; 00182 tmax=hit_time; 00183 } 00184 if(hit_time<tmin) tmin=hit_time; 00185 if(hit_time>tmax) tmax=hit_time; 00186 00187 //get begin_t and end_t 00188 if(i==0 && j==0) { 00189 begin_t=hit_time; 00190 end_t=hit_time; 00191 } 00192 if(hit_time<begin_t) begin_t=hit_time; 00193 if(hit_time>end_t) end_t=hit_time; 00194 00195 //upgrade the hits counter of the good chip at the good time 00196 add(channels_map, chip_id, hit_time); 00197 } 00198 00199 unsigned int timeold=0; 00200 00201 //Loop on each chip, at each time, to identify if this event is chip 00202 map<unsigned int, map<unsigned int,int> >::iterator iter1; 00203 for(iter1=channels_map.begin();iter1!=channels_map.end();iter1++){ 00204 00205 unsigned int chiip=iter1->first; 00206 map<unsigned int,int> hits_per_time=iter1->second; 00207 00208 map<unsigned int,int>::iterator iter2; 00209 for(iter2=hits_per_time.begin();iter2!=hits_per_time.end();iter2++){ 00210 00211 unsigned int time = iter2->first; 00212 int nb_hits = iter2->second; 00213 00214 //if this is a chip event 00215 if (nb_hits>=threshold) { 00216 00217 //set chip event indicator to 1 00218 //this.set->IsChip=true 00219 IsChip=1; 00220 00221 //set the time of this event 00222 ChipTime=time; 00223 00224 if(blabla==3) cout <<" -> event: "<<i+1<<" = chip in chip no "<<chiip<<" at time "<<time<<" with n_hits = "<<nb_hits<<" DELTA T= "<<time-timeold<<endl; 00225 timeold=time; 00226 00227 //update the event chip lists :event, chip, time, nb_hits 00228 list_chip_events[nb_tot][0] = i+1; 00229 list_chip_events[nb_tot][1] = chiip; 00230 list_chip_events[nb_tot][2] = time; 00231 list_chip_events[nb_tot][3] = nb_hits; 00232 if(time>tmax-timecut){ 00233 list_chip_events_timecut[nb_timecut][0] = i+1; 00234 list_chip_events_timecut[nb_timecut][1] = chiip; 00235 list_chip_events_timecut[nb_timecut][2] = time; 00236 list_chip_events_timecut[nb_timecut][3] = nb_hits; 00237 } 00238 00239 //update the basic event chip counters 00240 nb_tot++; 00241 if(time>tmax-timecut) nb_timecut++; 00242 00243 //update the basic chip event lists 00244 list_event[nb_tot-1] = i+1; 00245 if(time>tmax-timecut) list_event_timecut[nb_tot-1] = i+1; 00246 } 00247 } 00248 } 00249 00250 //if we want to plot this event 00251 if (IsChip && plot){ 00252 00253 TH1I *hdt = new TH1I("hdt","Hit time distribution (No time cut);N_hits;Time",tmax-tmin+3,-1,tmax-tmin+1); 00254 TH1I *hdt2 = new TH1I("hdt2","Hit time distribution (within time cut);N_hits;Time",timecut+3,tmax-timecut-1,tmax+1); 00255 TH1I *hdt3 = new TH1I("hdt3","Hit time distribution (around the time where it is chip);N_hits;Time",23,ChipTime-tmin-11,ChipTime-tmin+11); 00256 00257 //loop on all the channels of the event 00258 for (int j=0;j<nchannel;j++){ 00259 00260 //this channel 00261 MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j); 00262 00263 //at each channel get chip_id and hit_time 00264 UInt_t chip_id = channel->GetChipId(); 00265 UInt_t hit_time = channel->GetBcIdAbs() - (channel->GetBcIdDif() - channel->GetBcIdHit() ); 00266 00267 //fill the histogramms 00268 hxy->Fill(channel->GetY()/10000,channel->GetX()/10000); 00269 if(hit_time>tmax-timecut) hxy2->Fill(channel->GetY()/10000,channel->GetX()/10000); 00270 if(hit_time==ChipTime) hxy3->Fill(channel->GetY()/10000,channel->GetX()/10000); 00271 hdt->Fill(hit_time-tmin); 00272 if(hit_time>tmax-timecut) hdt2->Fill(hit_time-tmin); 00273 if(hit_time>ChipTime-100 && hit_time<ChipTime+500) hdt3->Fill(hit_time-tmin); 00274 } 00275 00276 //plots 00277 c0->Divide(3,2); 00278 c0->cd(1); 00279 hxy->GetXaxis()->SetTitle("X (cm)"); 00280 hxy->GetYaxis()->SetTitle("Y (cm)"); 00281 hxy->Draw("zcol"); 00282 00283 c0->cd(2); 00284 hxy2->GetXaxis()->SetTitle("X (cm)"); 00285 hxy2->GetYaxis()->SetTitle("Y (cm)"); 00286 hxy2->Draw("zcol"); 00287 00288 c0->cd(3); 00289 hxy3->GetXaxis()->SetTitle("X (cm)"); 00290 hxy3->GetYaxis()->SetTitle("Y (cm)"); 00291 hxy3->Draw("zcol"); 00292 00293 c0->cd(4); 00294 //hdt->SetAxisRange(-1,tmax-tmin+1); 00295 //hdt->SetMinimum(-1); 00296 //hdt->SetMaximum(tmax-tmin+1); 00297 hdt->GetXaxis()->SetTitle("time after trigger (BcId)"); 00298 hdt->GetYaxis()->SetTitle("Nb hits"); 00299 hdt->Draw(); 00300 00301 c0->cd(5); 00302 //hdt2->SetAxisRange(-1,timecut+1); 00303 //hdt2->SetMinimum(-1); 00304 //hdt2->SetMaximum(timecut+1); 00305 hdt2->GetXaxis()->SetTitle("time after trigger (BcId)"); 00306 hdt2->GetYaxis()->SetTitle("Nb hits"); 00307 hdt2->Draw(); 00308 00309 c0->cd(6); 00310 //hdt3->SetAxisRange(ChipTime-tmin-11,ChipTime-tmin+11); 00311 //hdt3->SetMinimum(ChipTime-tmin-11); 00312 //hdt3->SetMaximum(ChipTime-tmin+11); 00313 hdt3->GetXaxis()->SetTitle("time after trigger (BcId)"); 00314 hdt3->GetYaxis()->SetTitle("Nb hits"); 00315 hdt3->Draw(); 00316 00317 //save the plot 00318 char plotname[20] = "plot_evt"; 00319 char buf[10]; 00320 sprintf(buf,"%d",i+1); 00321 strcat(plotname, buf); 00322 strcat(plotname,".gif"); 00323 c0->SaveAs(plotname); 00324 00325 //clear the plot 00326 c0->Clear(); 00327 00328 //reset the histogramms 00329 hxy->Reset(); 00330 hxy2->Reset(); 00331 hxy3->Reset(); 00332 //hdt->Reset(); 00333 //hdt2->Reset(); 00334 //hdt3->Reset(); 00335 //hxy->Delete(); 00336 //hxy2->Delete(); 00337 //hxy3->Delete(); 00338 hdt->Delete(); 00339 hdt2->Delete(); 00340 hdt3->Delete(); 00341 } 00342 } 00343 00344 //Now the outputs: 00345 if(blabla==2||blabla==3) cout<<endl; 00346 if(blabla==1||blabla==2||blabla==3) { 00347 00348 //at each key print nb of evts chip 00349 cout<<"In key "<<counter<<":"<<endl<<endl;; 00350 cout<<" Nb chip evt = "<<nb_tot<<endl; 00351 cout<<" Nb chip evt / Nevent tot = "<<1.0*nb_tot/nEvent*100.<<" %"<<endl; 00352 cout<<" Nb chip evt within the timecut (delta t < "<<timecut<<") = "<<nb_timecut<<endl; 00353 float rate = float(nb_tot)/(end_t-begin_t)*1000000000/200*60; 00354 cout<<" Chip event rate = "<<rate<<" event/min"<<endl<<endl; 00355 00356 //at each key print print the liste of chip events 00357 cout<<"List of chip events: "<<endl; 00358 for(int it=0;it<nb_tot;it++) cout<<" "<<list_event[it]; 00359 cout<<endl<<endl; 00360 00361 //at each key print print the liste of chip events within the given timecut 00362 cout<<"List of chip events that are chip within the timecut (delta t < "<<timecut<<")"<<endl; 00363 for(int it=0;it<nb_timecut;it++) cout<<" "<<list_event_timecut[it]; 00364 cout<<endl<<endl; 00365 } 00366 00367 if(blabla==2||blabla==3) { 00368 00369 //if you want more details print also the position(chip) and time(memory order) of the chip events 00370 cout<<"Hereafter I show the list of chip events in each chip/chip at each time and their number of hits:"<<endl; 00371 cout<<"Event, chip, time, Nb_hits"<<endl; 00372 for(int it=0;it<nb_tot;it++) { 00373 for(int itt=0;itt<4;itt++) cout<<list_chip_events[it][itt]<<"; "; 00374 cout<<endl; 00375 } 00376 cout<<endl; 00377 cout<<"Hereafter I show the list of chip events that are chip within the timecut (delta t < "<<timecut<<"):"<<endl; 00378 cout<<"Event, chip, time, Nb_hits"<<endl; 00379 for(int it=0;it<nb_timecut;it++) { 00380 for(int itt=0;itt<4;itt++) cout<<list_chip_events_timecut[it][itt]<<"; "; 00381 cout<<endl; 00382 } 00383 cout<<endl; 00384 } 00385 cout<<endl; 00386 } 00387 theApp->Run(); 00388 delete theApp; 00389 return 0; 00390 }
1.4.7