Micromegas LAPP reconstruction

int main	(	int	argc,
		char **	argv
	)
Definition at line 33 of file Calib_asu_microroc.cpp.
00033                               {
00034 
00035   if (argc!=4){
00036     
00037     cout<<"/*****************************************************************/"<<endl;
00038     cout<<"/* ARGUMENTS TO THE PROGRAM **************************************/"<<endl;
00039     cout<<"/* 1 - DIRECTORY OF THE ROOT FILE PRODUCED BY THE RECONSTRUCTION */"<<endl;
00040     cout<<"/* 2 - NAME OF THE ROOT FILE *************************************/"<<endl;
00041     cout<<"/* 3 - THRESHOLD THAT IS MOVED DURING THE CALIBRATION ************/"<<endl;
00042     cout<<"/*    0 -> DAC_0 (LOW THRESHOLD) *********************************/"<<endl;
00043     cout<<"/*    1 -> DAC_1 (MEDIUM THRESHOLD) ******************************/"<<endl;
00044     cout<<"/*    2 -> DAC_2 (HIGH THRESHOLD) ********************************/"<<endl;
00045     cout<<"/*****************************************************************/"<<endl;}
00046 
00047   else{
00048 
00049   TString root_directory = argv[1];
00050   TString root_filename = argv[2];
00051   int dac = atoi(argv[3]);
00052 
00053   bool slab13_inversion = 0;
00054   float width_thr = 3;
00055   bool display = 0;
00056   bool print = 1;
00057 
00058   TString root_file = root_directory + root_filename;
00059   TString scurve_file = root_directory + "scurve_" + root_filename;
00060 
00061   cout<<endl;
00062   cout<<"Read file "<<root_file<<endl;
00063   cout<<"Wish to write "<<scurve_file<<endl;
00064   cout<<endl;
00065 
00066   int nbChannel = 0;
00067   float nbEvt = 0;
00068   int threshold = 0;
00069   int chipid = 0;
00070   int hardid = 0;
00071   int value = 0;
00072   int content = 0;
00073   int dac_min = 1025;
00074   int dac_max = 0;
00075   int scurve_max = 0;
00076   int deriv = 0;
00077 
00078   int N = 144;
00079  
00080   TH1I * hvalue = new TH1I("hvalue","",5,0,5);
00081   TH1I * hnpulse_threshold = new TH1I("hnpulse_threshold","",1024,0,1024);
00082 
00083   TH1F * hc_inflex = new TH1F("hc_inflex","Scurve derivative mean;channel number;mean (DAC)",9216,0,9216);
00084   TH1F * hc_width = new TH1F("hc_width","Scurve derivative RMS;channel number;RMS (DAC)",9216,0,9216);
00085   TH1F * hc_inflex_width = new TH1F("hc_inflex_width","Scurve derivative mean with RMS as error bar;channel number;mean (DAC)",9216,0,9216);
00086 
00087   TH1F * hinflex = new TH1F("hinflex","Scurve derivative mean;mean (DAC)",1000,0,1000);
00088   TH1F * hwidth  = new TH1F("hwidth","Scurve derivative RMS;RMS (DAC)",1000,0,100);
00089 
00090   hvalue->SetFillColor(2);
00091   hnpulse_threshold->SetFillColor(2);
00092   hc_inflex->SetLineColor(4);
00093   hc_width->SetLineColor(2);
00094   hinflex->SetFillColor(4);
00095   hwidth->SetFillColor(2);
00096   hinflex->SetLineColor(4);
00097   hwidth->SetLineColor(2);
00098 
00099   TH1I * hscurve[N][64];
00100   TH1I * hscurve_deriv[N][64];
00101 
00102   float inflex_tab[N][64];
00103   float inflex_err_tab[N][64];
00104   float width_tab[N][64];
00105   float width_err_tab[N][64];
00106   bool noisy[N][64];
00107 
00108   TString name,title;
00109 
00110   bool newchip[144];
00111   bool newchannel[144][64];
00112 
00113   for (int i=0;i<144;i++)
00114     {
00115       newchip[i] = 1;
00116       for (int j=0;j<64;j++)
00117         {
00118           inflex_tab[i][j] = 0;
00119           inflex_err_tab[i][j] = 0;
00120 
00121           width_tab[i][j] = 0;
00122           width_err_tab[i][j] = 0;
00123 
00124           noisy[i][j] = 0;
00125           newchannel[i][j] = 1;
00126         }
00127     }
00128 
00129 
00130 
00131   cout<<""<<endl;
00132   cout<<"******* READ CALIB DATA *******"<<endl;
00133   cout<<""<<endl;
00134   
00135   
00136   TFile f(root_file);
00137   TIter nextkey(f.GetListOfKeys());
00138   TKey *key;
00139   
00140   while (key = (TKey*)nextkey()) {
00141 
00142     TTree *tree = (TTree*)key->ReadObj();                
00143     MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");
00144                 if ( run != NULL )
00145                 {
00146                         const MTChip&  chip = run->GetOneChip();
00147                         try
00148                                 {       const MTMicrorocChip &microChip = dynamic_cast<const MTMicrorocChip &> (chip);
00149                                         if (dac==0){threshold = microChip.GetThresholdDac_0();}
00150                                         if (dac==1){threshold = microChip.GetThresholdDac_1();}
00151                                         if (dac==2){threshold = microChip.GetThresholdDac_2();}
00152                                 }
00153 
00154                         catch (...) {}
00155                 }
00156 
00157 
00158     cout<<"DAC value = "<<threshold<<" \r "<<flush;
00159 
00160     /*Find DAC range*/
00161     if (threshold<dac_min){dac_min = threshold;}
00162     if (threshold>dac_max){dac_max = threshold;}
00163 
00164     MTEvent *evt =  new MTEvent();
00165     TBranch *branch= tree->GetBranch("MTEvent");
00166     branch->SetAddress(&evt);
00167     
00168     MTChannel* channel = NULL;
00169     
00170     nbEvt = tree->GetEntries();
00171 
00172     hnpulse_threshold->SetBinContent(threshold+1,nbEvt);
00173 
00174 
00175     /*Create and fill histogram for each ASU channel*/
00176 
00177     for ( int evtNum = 0; evtNum < nbEvt ; evtNum++){
00178 
00179       tree->GetEntry(evtNum);
00180       nbChannel = evt->GetNchannel();
00181 
00182       for(int i=0;i<nbChannel ;i++){    
00183 
00184         channel = (MTChannel*)evt->GetChannels()->UncheckedAt(i);
00185         chipid = channel->GetChipId();
00186         hardid = channel->GetHardId();
00187 
00188         /*2nd microroc m2 prototype: Invert slab1 and slab3*/
00189 
00190         if (slab13_inversion){
00191           if (chipid<=48){chipid += 96;}
00192           else{
00193             if (chipid>96){chipid -= 96;}}}
00194 
00195         if (newchip[chipid-1]){
00196 
00197           cout<<"Chip "<<chipid<<" found"<<endl;
00198 
00199           newchip[chipid-1] = 0;}
00200 
00201         if (newchannel[chipid-1][hardid]){
00202 
00203           name  = Form("hscurve_%i_%i",chipid,hardid);
00204           title = Form("SCurve - chip %i - channel %i",chipid,hardid);
00205           hscurve[chipid-1][hardid] = new TH1I(name,title,1024,0,1024);
00206           hscurve[chipid-1][hardid]->SetFillColor(2);
00207 
00208           name  = Form("hscurve_deriv_%i_%i",chipid,hardid);
00209           title = Form("Differential SCurve - chip %i - channel %i",chipid,hardid);
00210           hscurve_deriv[chipid-1][hardid] = new TH1I(name,title,1024,0.5,1024.5);
00211           hscurve_deriv[chipid-1][hardid]->SetFillColor(4);
00212 
00213           value = channel->GetDigitalValue();
00214           hvalue->Fill(value);
00215           if (value==(dac+1)){hscurve[chipid-1][hardid]->Fill(threshold);}
00216 
00217           newchannel[chipid-1][hardid] = 0;}
00218 
00219         else{
00220 
00221           value = channel->GetDigitalValue();
00222           hvalue->Fill(value);
00223           if (value==(dac+1)){hscurve[chipid-1][hardid]->Fill(threshold);}}}}
00224 
00225 
00226     tree->Delete();}
00227 
00228 
00229 
00230   /*Calculate derivative of Scurve*/
00231 
00232   cout<<endl;
00233   cout<<"Threshold DAC range = ["<<dac_min<<";"<<dac_max<<"] DAC units"<<endl;
00234   cout<<endl;
00235 
00236   for (int i=0;i<144;i++){
00237     for (int j=0;j<64;j++){
00238       if (newchannel[i][j]==0){
00239 
00240         for (int k=dac_min+1;k<dac_max;k++){
00241 
00242           deriv = hscurve[i][j]->GetBinContent(k)-hscurve[i][j]->GetBinContent(k+1);
00243           if (deriv>0){hscurve_deriv[i][j]->SetBinContent(k,deriv);}}}}}
00244 
00245 
00246 
00247 
00248 
00249 
00250 
00251 
00252   cout<<""<<endl;
00253   cout<<""<<endl;
00254   cout<<"******* PROCESS CALIB DATA *******"<<endl;
00255   cout<<""<<endl;
00256 
00257 
00258   float inflex = 0;
00259   float width = 0;
00260   float inflex_err = 0;
00261   float width_err = 0;
00262 
00263 
00264   for (int i=0;i<N;i++)
00265     {
00266       for (int j=0;j<64;j++)
00267         {
00268           /*If SCurve exist*/
00269           if (newchannel[i][j]==0)
00270             {
00271               inflex = hscurve_deriv[i][j]->GetMean();
00272               inflex_err = hscurve_deriv[i][j]->GetMeanError();
00273               width = hscurve_deriv[i][j]->GetRMS();
00274               width_err = hscurve_deriv[i][j]->GetRMSError();
00275   
00276               hc_inflex_width->SetBinContent(i*64+j+1,inflex);
00277               hc_inflex_width->SetBinError(i*64+j+1,width);
00278 
00279               hc_inflex->SetBinContent(i*64+j+1,inflex);
00280               hc_inflex->SetBinError(i*64+j+1,inflex_err);
00281 
00282               hc_width->SetBinContent(i*64+j+1,width);
00283               hc_width->SetBinError(i*64+j+1,width_err);
00284 
00285               hinflex->Fill(inflex);
00286               hwidth->Fill(width);
00287 
00288               if (width>width_thr || hscurve[i][j]->GetBinContent(dac_max)!=0)
00289                 {
00290                   noisy[i][j] = 1;
00291                 }
00292             }
00293         }
00294     }
00295 
00296 
00297 
00298 
00299 
00300 
00301   
00302   
00303  
00304 
00305   cout<<""<<endl;
00306   cout<<""<<endl;
00307   cout<<"******* WRITE CALIB DATA *******"<<endl;
00308   cout<<""<<endl;
00309 
00310 
00311 
00312   TString outfile = root_directory + "scurve_" + root_filename;
00313   for (int i=0;i<4;i++){outfile.Chop();}
00314   outfile += "txt";
00315   ofstream out(outfile.Data());
00316 
00317   TFile * tf = new TFile(scurve_file,"RECREATE");
00318 
00319   hvalue->Write();
00320   hnpulse_threshold->Write();
00321 
00322   for (int i=0;i<N;i++)
00323     {
00324     if (!newchip[i])
00325       {     
00326         for (int j=0;j<64;j++)
00327           {
00328             if (!newchannel[i][j])
00329               {
00330                 hscurve[i][j]->Write();
00331                 hscurve_deriv[i][j]->Write();
00332 
00333                 inflex = hc_inflex->GetBinContent(i*64+j+1);
00334                 inflex_err = hc_inflex->GetBinError(i*64+j+1);
00335 
00336                 width = hc_width->GetBinContent(i*64+j+1);
00337                 width_err = hc_width->GetBinError(i*64+j+1);
00338 
00339                 cout<<i+1<<" "<<j<<" "<<inflex<<" "<<inflex_err<<" "<<width<<" "<<width_err<<" \r "<<flush;
00340                 out<<i+1<<" "<<j<<" "<<inflex<<" "<<inflex_err<<" "<<width<<" "<<width_err<<endl;                   
00341               }
00342           }
00343       }
00344     }
00345   
00346 
00347   cout<<endl;
00348   cout<<endl;
00349   cout<<"Noisy channel with pedestal RMS > "<<width_thr<<" DAC"<<endl;
00350 
00351   TString outfile2 = root_directory + "noisy_" + root_filename;
00352   for (int i=0;i<4;i++){outfile2.Chop();}
00353   outfile2 += "txt";
00354   ofstream out2(outfile2.Data());
00355 
00356   for (int i=0;i<N;i++)
00357     {
00358     if (!newchip[i])
00359       {     
00360         for (int j=0;j<64;j++)
00361           {
00362             if (!newchannel[i][j])
00363               {
00364                 if (noisy[i][j])
00365                   {
00366                     inflex = hc_inflex->GetBinContent(i*64+j+1);
00367                     inflex_err = hc_inflex->GetBinError(i*64+j+1);
00368 
00369                     width = hc_width->GetBinContent(i*64+j+1);
00370                     width_err = hc_width->GetBinError(i*64+j+1);
00371 
00372                     cout<<"  "<<i+1<<" "<<j<<" "<<inflex<<" "<<inflex_err<<" "<<width<<" "<<width_err<<endl;
00373                     out2<<i+1<<" "<<j<<" "<<inflex<<" "<<inflex_err<<" "<<width<<" "<<width_err<<endl;
00374                   }
00375               }
00376           }
00377       }
00378     }
00379   
00380 
00381 
00382   hc_inflex_width->Write();
00383   hc_inflex->Write();
00384   hc_width->Write();
00385   hinflex->Write();
00386   hwidth->Write();
00387 
00388 
00389 
00390 
00391   if (print)
00392     {
00393 
00394       cout<<""<<endl;
00395       cout<<""<<endl;
00396       cout<<"******* PRINT CALIB DATA *******"<<endl;
00397       cout<<""<<endl;
00398 
00399       int channel_min = 0;
00400       int channel_max = 0;
00401 
00402       float value_min = 0;
00403       float value_max = 0;
00404 
00405       for (int i=0;i<9216;i++)
00406         {
00407           if (hc_inflex->GetBinContent(i+1)!=0)
00408             {
00409               channel_min = i-100;
00410               i=9216;
00411             }
00412         }
00413 
00414       for (int i=9216;i>0;i--)
00415         {
00416           if (hc_inflex->GetBinContent(i)!=0)
00417             {
00418               channel_max = i-1+100;
00419               i=0;
00420             }
00421         }
00422 
00423       title = Form("Summary plots of calibration run : %s",root_filename.Data());
00424       TCanvas * c0 = new TCanvas("c0",title,10,10,1200,800);
00425       c0->Divide(2,2);
00426 
00427       c0->cd(1);
00428       value_min = hinflex->GetMean() - hinflex->GetRMS()*5;
00429       value_max = hinflex->GetMean() + hinflex->GetRMS()*5;
00430       hc_inflex->GetYaxis()->SetTitleOffset(1.5);
00431       hc_inflex->GetXaxis()->SetRangeUser(channel_min-10,channel_max+10);
00432       hc_inflex->GetYaxis()->SetRangeUser(value_min,value_max);
00433       hc_inflex->Draw();
00434 
00435       c0->cd(3);
00436       hinflex->GetXaxis()->SetRangeUser(value_min,value_max);
00437       hinflex->Draw();
00438 
00439       c0->cd(2);
00440       value_min = 0;
00441       value_max = hwidth->GetMean() + hwidth->GetRMS()*5;
00442       hc_width->GetYaxis()->SetTitleOffset(1.5);
00443       hc_width->GetXaxis()->SetRangeUser(channel_min-10,channel_max+10);
00444       hc_width->GetYaxis()->SetRangeUser(value_min,value_max);
00445       hc_width->Draw();
00446       
00447       c0->cd(4);
00448       hwidth->GetXaxis()->SetRangeUser(value_min,value_max);
00449       hwidth->Draw();
00450       
00451       TString print_filename = root_directory + "plots_" + root_filename;
00452       for (int i=0;i<4;i++){print_filename.Chop();}
00453       print_filename += "png";
00454 
00455       cout<<"Print file "<<print_filename<<endl;
00456       c0->Print(print_filename);
00457     }
00458 
00459 
00460   tf->Close();
00461 
00462   cout<<endl;
00463   cout<<endl;
00464   }
00465 
00466 
00467 
00468 
00469     return 0 ;}
/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/Microroc/Calibration/Calib_asu_microroc.cpp File Reference

Functions

Function Documentation
