Micromegas LAPP reconstruction

00001 /* @version $Revision: 1328 $ * @modifiedby $Author: jacquem $ * @lastmodified $Date: 2011-10-03 17:04:17 +0200 (Mon, 03 Oct 2011) $ */
00002 #include <map>
00003 #include "TKey.h"
00004 #include <TROOT.h>
00005 #include <TRint.h>
00006 #include <TKey.h>
00007 #include <TTree.h>
00008 #include <TFile.h>
00009 #include <TSystem.h>
00010 #include <iostream>
00011 #include <TF1.h>
00012 #include <TH2.h>
00013 #include <TGraph.h>
00014 #include <sstream>
00015 #include "Log.hh"
00016 #include <fstream>
00017 #include "TCanvas.h"
00018 #include "TApplication.h"
00019 #include "TColor.h"
00020 #include "TStyle.h"
00021 #include "TAxis.h"
00022 #include "TString.h"
00023 #include "TLegend.h"
00024 #include "Bytes.h"
00025 
00026 
00027 #include "root/MTRun.hh"
00028 #include "root/MTEvent.hh"
00029 #include "root/MTChannel.hh"
00030 #include "root/MTMicrorocChip.hh"
00031 
00032 #include "MicroException.hh"
00033 
00034 //#define DISPLAY
00035 #undef DISPLAY
00036 //#define R_DEBUG
00037 #undef R_DEBUG
00038 //#define PLOT
00039 #undef PLOT
00040 #define WRITE
00041 
00042 using namespace std;
00043 void usage(char *progname) ;
00044 #ifdef PLOT
00045         void set_style_myplain() ;
00046         TStyle *style_myplain=new TStyle("myplain","classic style");
00047 #endif
00048 ////////////////////////////////////////////////////////////////////////////////////////////
00049 
00050 
00051 
00052 
00053 
00054 int main(int argc, char**argv){
00055 
00056   /*****************************************************************/
00057   /* ARGUMENTS TO THE PROGRAM **************************************/
00058   /* 1 - DIRECTORY OF THE ROOT FILE PRODUCED BY THE RECONSTRUCTION */
00059   /* 2 - NAME OF THE ROOT FILE *************************************/
00060   /* 3 - THRESHOLD THAT IS MOVED DURING THE CALIBRATION ************/
00061   /*    0 -> DAC_0 (LOW THRESHOLD) *********************************/
00062   /*    1 -> DAC_1 (MEDIUM THRESHOLD) ******************************/
00063   /*    2 -> DAC_2 (HIGH THRESHOLD) ********************************/
00064   /*****************************************************************/
00065         if(argc != 4)
00066         {       usage(argv[0]) ;
00067                 return 1 ;
00068         }
00069         TString root_directory = argv[1];
00070         TString root_filename = argv[2];
00071         int dac = atoi(argv[3]);
00072 
00073         //bool display = 0;
00074         #ifdef PLOT
00075                 set_style_myplain() ;
00076         #endif
00077         TCanvas *c1 = new TCanvas("c1","Test MICROROC",0,0,1000,1000);
00078         TString root_file = root_directory + root_filename;
00079         TString scurve_file = root_directory + "scurve_" + root_filename;
00080 
00081 //  cout<<"Wish to write "<<scurve_file<<endl;
00082         
00083         
00084 //      cout << "Attention ! This version is modified for temperature testing ! Else, comment line 123 of source !" << endl ;   
00085          
00086         int nbChannel = 0;
00087         float nbEvt = 0;
00088         int threshold = 0;
00089         //int chipid = 0;
00090         int hardid = 0;
00091         int value = 0;
00092         int serial = 0;
00093         //int content = 0;
00094         Double_t content = 0;
00095         UInt_t chipId ;
00096     int dac_min = 1025;
00097         int dac_max = 0;
00098         
00099         TH1I * hvalue = new TH1I("hvalue","",5,0,5);
00100         TH1I * hnpulse_threshold = new TH1I("hnpulse_threshold","",1024,0,1024);
00101 
00102         TH1I * hscurve[64] ;
00103         TH1I * hscurve_deriv[64] ;
00104         TString name, title, name_deriv, title_deriv ;
00105  
00106         for (int j=0;j<64;j++)
00107         {       name  = Form("hscurve_%i",j);
00108                 title = Form("SCurve from channel %i",j);
00109                 hscurve[j] = new TH1I(name,title,1024,0,1024);
00110                 name_deriv  = Form("deriv_%i",j);
00111                 title_deriv = Form("Derivative from channel %i",j);
00112                 hscurve_deriv[j] = new TH1I(name_deriv,title_deriv,1024,0,1024);
00113         }
00114   
00115         TFile f(root_file);
00116         TIter nextkey(f.GetListOfKeys());
00117         TKey *key;
00118         while (key = (TKey*)nextkey())
00119         {       TTree *tree = (TTree*)key->ReadObj();                      
00120                 MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");
00121                 const MTChip&  chip = run->GetOneChip();
00122                 //    cout << "------------ Chip[" << chip.GetId() << "], softId: [" << chip.GetSoftId() <<"] configuration ---------------" << endl ;
00123                 
00124                 serial = chip.GetSerialNumber() ;
00125 //              cout << serial << " " ;
00126 //              if (serial == 14) {
00127                 try
00128                 {       const MTMicrorocChip &microChip = dynamic_cast<const MTMicrorocChip &> (chip);
00129                         //chipId = chip.GetId() ;
00130                         //cout << "Chip identifier: " << chipId << endl ;
00131                         if (dac==0){threshold = microChip.GetThresholdDac_0();}
00132                         if (dac==1){threshold = microChip.GetThresholdDac_1();}
00133                         if (dac==2){threshold = microChip.GetThresholdDac_2();}
00134                 }
00135                 catch (...) {}
00136                 /*Find DAC range*/
00137         if (threshold<dac_min){dac_min = threshold;}
00138         if (threshold>dac_max){dac_max = threshold;}
00139 
00140                 MTEvent *evt =  new MTEvent();
00141                 TBranch *branch= tree->GetBranch("MTEvent");
00142                 branch->SetAddress(&evt);
00143                 MTChannel* channel = NULL;
00144                 nbEvt = tree->GetEntries();
00145 //              cout << "Number of events in the tree: " << nbEvt << endl ;
00146 
00147                 hnpulse_threshold->SetBinContent(threshold+1,nbEvt);
00148 
00149                 for( int evtNum = 0; evtNum < nbEvt ; evtNum++)
00150                 {       //if (evtNum%1000 == 0) {cout<<evtNum/nbEvt*100.<<" %"<<"\r"<<flush;}
00151                         tree->GetEntry(evtNum);
00152                         nbChannel = evt->GetNchannel();
00153                         for(int i=0;i<nbChannel ;i++)
00154                         {       channel = (MTChannel*)evt->GetChannels()->UncheckedAt(i);
00155                                 chipId = channel->GetChipId();
00156                                 hardid = channel->GetHardId();
00157                                 value = channel->GetDigitalValue();
00158                                 hvalue->Fill(value);
00159                                 if (value==(dac+1)) // ie 1 for low, 2 for med and 3 for high ....
00160                                 {       hscurve[hardid]->Fill(threshold);
00161                                         #ifdef R_DEBUG
00162                                                 cout << "---- digital value:" <<  value <<endl;
00163                                                 cout << "---- hard id:" << hardid <<endl;
00164                                                 //cout << "---- x:"<< channel->GetX() <<endl;
00165                                                 //cout << "---- y:" << channel->GetY() <<endl;
00166                                                 cout << "---- chip id:" << chipId <<endl;
00167                                         #endif
00168           
00169                                 }
00170                         }
00171                 }
00172         }
00173         
00174 
00175         
00176         //===========================================================================*
00177         // Now analyze scurves  
00178         //===========================================================================*
00179 
00180         TString t_scurve_file = root_directory + root_filename ;//+ ".root";
00181         //cout << t_scurve_file<< endl;
00182 
00183         int tmin = 0;
00184         int tmed = 0;
00185         int tmax = 0;
00186         int delta = 0;
00187 
00188 //      TH1I * hscurve[64];
00189 
00190 /*      TH1F * hinflex = new TH1F("hinflex","",1024,0,1024);
00191         TH1F * hsigma = new TH1F("hsigma","",50,-2,2);
00192 
00193         TH1F * hinflex_err = new TH1F("hinflex_err","",100,0,10);
00194         TH1F * hsigma_err = new TH1F("hsigma_err","",100,0,2);
00195 */
00196         TH1F * hc_inflex = new TH1F("hc_inflex","",64,0,64);
00197         TH1F * hc_sigma = new TH1F("hc_sigma","",64,0,64);
00198 
00199         TH1F * hc_inflex_err = new TH1F("hc_inflex_err","",64,0,64);
00200         TH1F * hc_sigma_err = new TH1F("hc_sigma_err","",64,0,64);
00201 
00202         TH1F * hc_inflex_sigma = new TH1F("hc_inflex_sigma","",64,0,64);
00203 
00204 //      TString name;
00205 /*      for (int j=0;j<64;j++)
00206         {       name  = Form("hscurve_%i",j);
00207                 TH1I * h = (TH1I*)tf.Get(name);
00208                 hscurve[j] = h;
00209         }*/
00210 
00211         float scurve_max = 0;
00212         float scurve_inflex = 0;
00213         float inflex = 0;
00214         float sigma = 0;
00215         float inflex_err = 0;
00216         float sigma_err = 0;
00217         float chi2 = 0;
00218         float ndf = 0;
00219         Double_t deriv = 0;
00220         
00221         bool start_100 ; // test if a channel starts at 100% else warning
00222         bool has_flipped ; // test if each channel flips from 100% to 0% else, warning message
00223         Int_t noisy;    // test if each channel is noisy, ie has entry beyond flip point
00224         Int_t noise_limit = 2;  // number of entry accepted
00225         Int_t noise_window = 10 ;// number of bin defore dac_max to verify noise
00226         Int_t flip_limit = 30 ;
00227           
00228         for (int j=0;j<64;j++)
00229         {       noisy = 0 ;
00230                 for (int k=dac_min+1;k<=dac_max;k++) // ATTN : bin = dac + 1 !!
00231                 {       deriv = ( hscurve[j]->GetBinContent(k)) - (hscurve[j]->GetBinContent(k+1));
00232                         //cout << j << " " << k << " " << hscurve[j]->GetBinContent(k) << " " << hscurve[j]->GetBinContent(k+1) << " " << deriv << endl ;
00233                         // test if each channel has flipped
00234                         if  (k==dac_min+1) // first loop
00235                         {       if (hscurve[j]->GetBinContent(k)>=100)
00236                                         start_100 = 1 ;
00237                                 else start_100 = 0 ;
00238                         }
00239                         
00240                         if  (k>(dac_max - noise_window)) // test noisy channel
00241                         {       if (hscurve[j]->GetBinContent(k)>=noise_limit)
00242                                         noisy++ ;
00243                                 //cout << "Channel " << j << " dac " << k << " noise " << noisy << endl ;
00244                         }
00245                         
00246                         
00247                         
00248                         if  (k==dac_max-1) // last loop has flipped ?
00249                         {       if (hscurve[j]->GetBinContent(k)<=flip_limit && start_100 == 1 )
00250                                 {       has_flipped = 1 ;
00251                                         //if (hscurve[j]->GetBinContent(k)>0)
00252                                         //      noisy++;// = 0 ;
00253                                         //else noisy = 1 ;
00254                                 }
00255                                 else has_flipped = 0 ;
00256                         }                       
00257                         if (deriv > 0.)
00258                         {       hscurve_deriv[j]->SetBinContent(k,deriv);
00259                                 hscurve_deriv[j]->SetEntries((hscurve_deriv[j]->GetEntries())+deriv);
00260                                 //hscurve_deriv[j]->Fill(k,deriv);
00261                         }
00262                 }
00263                 if (has_flipped == 0)
00264                 {       cout << "****** " << argv[0] << ": Warning! Chip " << serial << " Channel " << j << " has not flipped" ;
00265                         if (start_100 == 1) cout << " (stuck to 100 %)" << endl ;
00266                         else cout << " (stuck to 0 %)" << endl ;
00267                 }
00268                 if (noisy==noise_window) cout << "****** " << argv[0] << ": Warning! Chip " << serial << " Channel " << j << " may be noisy, check root scurve file !" << endl ;
00269                 
00270          }
00271 
00272 // FIT
00273         TF1 * gaussfunc[64];
00274         for (int j=0;j<64;j++)
00275         {       name = Form("g_%i",j);
00276                 gaussfunc[j] = new TF1(name,"gaus",0,1024);
00277         }
00278         
00279         for (int j=0;j<64;j++)
00280         {       //gaussfunc[j]->SetParLimits(1, dac_min, dac_max) ;
00281                 gaussfunc[j]->SetParameter(0, hscurve_deriv[j]->GetMaximumBin());
00282                 gaussfunc[j]->SetParLimits(0, 1, 100);
00283                 gaussfunc[j]->SetParameter(1, hscurve_deriv[j]->GetMean() ); //GetMaximumBin()
00284                 //gaussfunc[j]->FixParameter(1, hscurve_deriv[j]->GetMean() );
00285                 gaussfunc[j]->SetParLimits(1, dac_min, dac_max) ; //hscurve_deriv[j]->GetMaximumBin()-10, hscurve_deriv[j]->GetMaximumBin() + 10) ;
00286                 gaussfunc[j]->SetParameter(2, hscurve_deriv[j]->GetRMS() );
00287                 gaussfunc[j]->SetParLimits(2, 0.01, hscurve_deriv[j]->GetRMS()) ;
00288 //              if (j==48)
00289 //              hscurve_deriv[j]->Fit(gaussfunc[j],"BV", "",hscurve_deriv[j]->GetMean()-15, hscurve_deriv[j]->GetMean() + 15 );
00290 //              else
00291                                 hscurve_deriv[j]->Fit(gaussfunc[j],"WBQ", "",hscurve_deriv[j]->GetMean()-20, hscurve_deriv[j]->GetMean() + 20 );
00292                 //hscurve_deriv[j]->Fit(gaussfunc[j],"Q");
00293                 inflex = gaussfunc[j]->GetParameter(1);
00294                 if (inflex<=dac_min || inflex>=dac_max)
00295                 {       cout << "****** " << argv[0] << ": Warning! Chip " << serial << " Channel " << j << ", fit may be not correct ! Inflexion: " << inflex << " Dac min: " << dac_min << " Dac max: " << dac_max << endl ;
00296                 }
00297                 sigma = gaussfunc[j]->GetParameter(2);
00298                 inflex_err = gaussfunc[j]->GetParError(1);
00299                 sigma_err =gaussfunc[j]->GetParError(2);
00300                 
00301                 hc_inflex->SetBinContent(j+1,inflex);
00302                 hc_sigma->SetBinContent(j+1,sigma);             
00303 
00304                 hc_inflex_err->SetBinContent(j+1,inflex_err);
00305                 hc_sigma_err->SetBinContent(j+1,sigma_err);
00306 
00307                 //cout << "Channel " << j << " init at " << hscurve_deriv[j]->GetMean() << " " << gaussfunc[j]->GetParameter(1) << endl ;
00308                 // cout << j << " " << inflex << " " << sigma << endl ;
00309         }
00310 
00311 
00312         //////////////////////////////////////////////////////////////////////////////////////////
00313         // Write Scurves
00314         //////////////////////////////////////////////////////////////////////////////////////////
00315         TFile * tf = new TFile(scurve_file,"RECREATE");
00316         hvalue->Write();
00317         hnpulse_threshold->Write();
00318         for (int j=0;j<64;j++)
00319         {       hscurve[j]->Write();
00320                 hscurve_deriv[j]->Write();
00321                 //gaussfunc[j]->Write();
00322         }
00323         tf->Close();
00324         
00325         
00326         #ifdef WRITE
00327                 root_filename.ReplaceAll(".root", "");
00328                 //TString outfile = root_directory + "../RESULTS/" +  "parameters_" + root_filename + ".txt";
00329                 TString outfile = root_directory + "../RESULTS/" +  root_filename + ".txt";
00330                 ofstream out ;
00331                 out.open(outfile) ;
00332                 for (int j=0;j<64;j++)
00333                 {       out<<j<<" "<<hc_inflex->GetBinContent(j+1)<<" "<<hc_inflex_err->GetBinContent(j+1)<<" "<<hc_sigma->GetBinContent(j+1)<<" "<<hc_sigma_err->GetBinContent(j+1)<<endl;
00334                 }
00335         #endif
00336         #ifdef PLOT
00337                 TApplication theApp("App", &argc, argv);
00338                 TString adj;
00339                 if (dac==0){adj = "Low";}
00340                 if (dac==1){adj = "Medium";}
00341                 if (dac==2){adj = "High";}
00342                 
00343                 c1->Divide(8,8);
00344                 //c1->cd();
00345                 TString htitle = adj + " threshold";
00346                 for (int j=0;j<64;j++)
00347                 {       //c0->cd(j/8,j%8);
00348                         c1->cd(j+1);
00349                         hscurve[j]->SetLineWidth(2);
00350                         hscurve[j]->Draw("");
00351                         c1->Update() ;
00352                 }
00353                 theApp.Run(kTRUE) ;
00354                 c1->Update() ;
00355         #endif
00356     return 0 ;
00357 }
00358 
00359 
00360 
00361 //**********************************************************************************************
00362 //**********************************************************************************************
00363 
00364 
00365 void usage(char *progname)
00366 {
00367         printf("Usage : %s directory rootfile [0|1|2]", progname);
00368         printf("        Where [0|1|2] is the dac number\n");
00369         printf("        scurve_rootfile.root is created\n");
00370 }
00371 
00372 //**********************************************************************************************
00373 //**********************************************************************************************
00374 #ifdef PLOT
00375 void set_style_myplain()
00376 {
00377         style_myplain->SetCanvasBorderMode(0);
00378         style_myplain->SetCanvasColor(0);
00379         style_myplain->SetDrawBorder(0);
00380         style_myplain->SetPadBorderMode(0);
00381         style_myplain->SetPadColor(10);
00382         style_myplain->SetFrameLineColor(1);
00383         style_myplain->SetFrameFillColor(5);
00384         style_myplain->SetFrameFillStyle(0);
00385         style_myplain->SetFrameBorderMode(0);
00386 
00387         //style_myplain->SetFrameLineWidth(5); //test for setframelinecolor !
00388 
00389         style_myplain->SetLegendBorderSize(1);
00390 
00391         style_myplain->SetTextColor(231);
00392 
00393 //      style_myplain->SetFillColor(0); // If enable : impossible to fill histo with color !!! see note in the end !
00394         style_myplain->SetLineColor(231);
00395 
00396         style_myplain->SetStatColor(10);
00397         style_myplain->SetStatTextColor(1);
00398         style_myplain->SetStatBorderSize(1) ;
00399         style_myplain->SetStatX(0.9);
00400         style_myplain->SetStatY(0.9);
00401         style_myplain->SetOptStat("");//"emr" by default
00402         style_myplain->SetOptFit(0111);
00403 
00404         style_myplain->SetTitleTextColor(1);
00405         style_myplain->SetTitleBorderSize(0);
00406         style_myplain->SetTitleAlign(23);
00407         style_myplain->SetTitleX(0.5);
00408 
00409         style_myplain->SetAxisColor(1,"x");     // bleu
00410         style_myplain->SetAxisColor(1,"y");
00411         style_myplain->SetAxisColor(1,"z");
00412         style_myplain->SetLabelColor(1,"x");
00413         style_myplain->SetLabelColor(1,"x");
00414         style_myplain->SetLabelColor(1,"y");
00415         style_myplain->SetLabelColor(1,"z");
00416         style_myplain->SetLabelSize(0.03,"x");
00417         style_myplain->SetLabelSize(0.03,"y");
00418         style_myplain->SetLabelSize(0.03,"z");
00419         style_myplain->SetTitleColor(1,"x");
00420         style_myplain->SetTitleColor(1,"y");
00421         style_myplain->SetTitleColor(1,"z");
00422         style_myplain->SetTitleFillColor(10);
00423 
00424         style_myplain->SetHistFillColor(10);
00425         //style_myplain->SetHistFillStyle(0);
00426         style_myplain->SetHistFillStyle(1001);
00427         style_myplain->SetHistLineColor(1);                     // 231 !!!
00428 
00429         style_myplain->SetFuncColor(kOrange+10);        // dark orange
00430         style_myplain->SetPalette(1);                           // temp�rature like
00431         
00432 //      style_myplain->SetPaperSize(100,100);  // for Hi Res plots only !
00433 
00434         style_myplain->cd();
00435 
00436         gROOT->SetStyle("myplain"); 
00437 
00438         gROOT->ForceStyle(true);
00439 
00440 
00441 
00442 // frame style not applied ! why ????????
00443 // do it manually : right click on frame then "use current style"
00444 }
00445 #endif
/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/Renaud/test_microroc_2.cpp
