Micromegas LAPP reconstruction

00001 /* @version $Revision: 1585 $ * @modifiedby $Author: jacquem $ * @lastmodified $Date: 2012-03-13 11:03:35 +0100 (Tue, 13 Mar 2012) $ */
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 #include "TPaveStats.h"
00026 
00027 #include "TLine.h"
00028 #include "TPaveLabel.h"
00029 
00030 
00031 #include "root/MTRun.hh"
00032 #include "root/MTEvent.hh"
00033 #include "root/MTChannel.hh"
00034 #include "root/MTMicrorocChip.hh"
00035 
00036 #include "MicroException.hh"
00037 
00038 //#define DISPLAY
00039 #undef DISPLAY
00040 //#define R_DEBUG
00041 #undef R_DEBUG
00042 #define PLOT
00043 //#undef PLOT
00044 //define WRITE
00045 #undef WRITE
00046 #define INFORENO
00047 //#undef INFORENO
00048 
00049 using namespace std;
00050 void usage(char *progname) ;
00051 #ifdef PLOT
00052         void set_style_myplain() ;
00053         TStyle *style_myplain=new TStyle("myplain","classic style");
00054 #endif
00055 ////////////////////////////////////////////////////////////////////////////////////////////
00056 
00057 int main(int argc, char**argv){
00058 
00059   /*****************************************************************/
00060   /* ARGUMENTS TO THE PROGRAM **************************************/
00061   /* 1 - DIRECTORY OF THE ROOT FILE PRODUCED BY THE RECONSTRUCTION */
00062   /* 2 - NAME OF THE ROOT FILE *************************************/
00063   /* 3 - THRESHOLD THAT IS MOVED DURING THE CALIBRATION ************/
00064   /*    0 -> DAC_0 (LOW THRESHOLD) *********************************/
00065   /*    1 -> DAC_1 (MEDIUM THRESHOLD) ******************************/
00066   /*    2 -> DAC_2 (HIGH THRESHOLD) ********************************/
00067   /*****************************************************************/
00068         if(argc != 5)
00069         {       usage(argv[0]) ;
00070                 return 1 ;
00071         }
00072         TString root_directory = argv[1];
00073         TString root_filename = argv[2];
00074         int dac = atoi(argv[3]);        
00075         int chan = atoi(argv[4]);       // chip channel
00076         
00077         //bool display = 0;
00078 /*      #ifdef PLOT
00079                 set_style_myplain() ;
00080         #endif
00081 */
00082         TString root_file = root_directory + root_filename;
00083         TString scurve_file = root_directory + "scurve_" + root_filename;
00084 
00085 //  cout<<"Wish to write "<<scurve_file<<endl;
00086 
00087         int nbChannel = 0;
00088         float nbEvt = 0;
00089         int threshold = 0;
00090         //int chipid = 0;
00091         int hardid = 0;
00092         int value = 0;
00093         //int content = 0;
00094         Double_t content = 0;
00095         UInt_t chipId ;
00096  
00097 //      TH1I * hvalue = new TH1I("hvalue","",5,0,5);
00098 //      TH1I * hnpulse_threshold = new TH1I("hnpulse_threshold","",1024,0,1024);
00099 
00100 //      TH1I * hscurve[64];
00101         TString name,title;
00102  
00103 //      for (int j=0;j<64;j++)
00104 //      {       name  = Form("hscurve_%i",j);
00105 //              title = Form("SCurve from channel %i",j);
00106 //              hscurve[j] = new TH1I(name,title,1024,0,1024);
00107 //      }
00108 
00109         name  = Form("hscurve_%i", chan);
00110         title = Form("SCurve from channel %i", chan);
00111         TH1I * hscurve = new TH1I(name,"",1025,0,1025);
00112 
00113         TFile f(root_file);
00114         TIter nextkey(f.GetListOfKeys());
00115         TKey *key;
00116         #ifdef INFORENO
00117                 cout << "Reading tree" << endl ;
00118         #endif
00119         while (key = (TKey*)nextkey())
00120         {       
00121                 TTree *tree = (TTree*)key->ReadObj();                      
00122                 MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");
00123                 const MTChip&  chip = run->GetOneChip();
00124 
00125                 try
00126                 {       const MTMicrorocChip &microChip = dynamic_cast<const MTMicrorocChip &> (chip);
00127                         //const MTHardroc2Chip &microChip = dynamic_cast<const MTHardroc2Chip &> (chip);
00128                         //cout << "threshold dac 0:[" << hr2Chip.GetThresholdDac_0()  << endl;
00129                         //cout << "threshold dac 1:[" << hr2Chip.GetThresholdDac_1()  << endl;
00130                         //cout << "threshold dac 2:[" << hr2Chip.GetThresholdDac_2()  << endl;
00131                         //chipId = chip.GetId() ;
00132                         //cout << "Chip identifier: " << chipId << endl ;
00133                         if (dac==0){threshold = microChip.GetThresholdDac_0();}
00134                         if (dac==1){threshold = microChip.GetThresholdDac_1();}
00135                         if (dac==2){threshold = microChip.GetThresholdDac_2();}
00136                 }
00137                 catch (...) {}
00138 
00139                 MTEvent *evt =  new MTEvent();
00140                 TBranch *branch= tree->GetBranch("MTEvent");
00141                 branch->SetAddress(&evt);
00142                 MTChannel* channel = NULL;
00143                 nbEvt = tree->GetEntries();
00144 //              cout << "Number of events in the tree: " << nbEvt << endl ;
00145 
00146 //              hnpulse_threshold->SetBinContent(threshold+1,nbEvt);
00147 
00148                 for( int evtNum = 0; evtNum < nbEvt ; evtNum++)
00149                 {       //if (evtNum%1000 == 0) {cout<<evtNum/nbEvt*100.<<" %"<<"\r"<<flush;}
00150                         tree->GetEntry(evtNum);
00151                         nbChannel = evt->GetNchannel();
00152                         for(int i=0;i<nbChannel ;i++)
00153                         {       channel = (MTChannel*)evt->GetChannels()->UncheckedAt(i);
00154                                 chipId = channel->GetChipId();
00155                                 hardid = channel->GetHardId();
00156                                 value = channel->GetDigitalValue();
00157 //                              hvalue->Fill(value);
00158                                 if (value==(dac+1))
00159                                 {       hscurve->Fill(threshold);
00160                                         #ifdef R_DEBUG
00161                                                 cout << "---- digital value:" <<  value <<endl;
00162                                                 cout << "---- hard id:" << hardid <<endl;
00163                                                 cout << "---- x:"<< channel->GetX() <<endl;
00164                                                 cout << "---- y:" << channel->GetY() <<endl;
00165                                                 cout << "---- chip id:" << chipId <<endl;
00166                                         #endif
00167           
00168                                 }
00169                         }
00170                 }
00171         }
00172         #ifdef INFORENO
00173                 cout << "Tree has been read !" << endl ;
00174         #endif
00175         
00176         //===========================================================================*
00177         // Now analyze scurves  
00178         #ifdef INFORENO
00179                 cout << "Fitting" << endl ;
00180         #endif
00181         TString t_scurve_file = root_directory + root_filename ;//+ ".root";
00182         //cout << t_scurve_file<< endl;
00183 
00184         int tmin = 0;
00185         int tmed = 0;
00186         int tmax = 0;
00187         int delta = 0;
00188 
00189 //      TH1I * hscurve[64];
00190 
00191         TH1F * hinflex = new TH1F("hinflex","",1024,0,1024);
00192         TH1F * hsigma = new TH1F("hsigma","",50,-2,2);
00193 
00194         TH1F * hinflex_err = new TH1F("hinflex_err","",100,0,10);
00195         TH1F * hsigma_err = new TH1F("hsigma_err","",100,0,2);
00196 
00197         TH1F * hc_inflex = new TH1F("hc_inflex","",64,0,64);
00198         TH1F * hc_sigma = new TH1F("hc_sigma","",64,0,64);
00199 
00200         TH1F * hc_inflex_err = new TH1F("hc_inflex_err","",64,0,64);
00201         TH1F * hc_sigma_err = new TH1F("hc_sigma_err","",64,0,64);
00202 
00203         TH1F * hc_inflex_sigma = new TH1F("hc_inflex_sigma","",64,0,64);
00204 
00205 //      TString name;
00206 
00207 /*      for (int j=0;j<64;j++)
00208         {       name  = Form("hscurve_%i",j);
00209                 TH1I * h = (TH1I*)tf.Get(name);
00210                 hscurve[j] = h;
00211         }
00212 */
00213         TString func = "[0]/(1+exp((x-[1])/[2]))";
00214         TF1 * fermi[64];
00215         float scurve_max = 0;
00216         float scurve_inflex = 0;
00217         float inflex = 0;
00218         float sigma = 0;
00219         float inflex_err = 0;
00220         float sigma_err = 0;
00221         float content0,content1,content2;
00222         float chi2 = 0;
00223         float ndf = 0;
00224 
00225 
00226 #ifdef PLOT
00227                 set_style_myplain() ;
00228    TApplication theApp("App", &argc, argv);
00229 
00230 //   TCanvas *c = new TCanvas("c", "The Hello Canvas", 400, 400);
00231 
00232 //   TPaveLabel *hello = new TPaveLabel(0.2,0.4,0.8,0.6,"Hello World");
00233 //   hello->Draw();
00234 //   TPaveLabel *quit = new TPaveLabel(0.2,0.2,0.8,0.3,"Close via menu File/Quit");
00235 //   quit->Draw();
00236 //   c->Update();
00237 
00238    // Enter event loop, one can now interact with the objects in
00239    // the canvas. Select "Exit ROOT" from Canvas "File" menu to exit
00240    // the event loop and execute the next statements.
00241 //   theApp.Run(kTRUE);
00242 
00243 //   TLine *l = new TLine(0.1,0.2,0.5,0.9);
00244 //   l->Draw();
00245 //   quit->SetLabel("Select File/Quit again to exit app");
00246 //   c->Update();
00247 
00248    // Here we don't return from the eventloop. "Exit ROOT" will quit the app.
00249 //   theApp.Run();
00250 
00251 
00252 
00253 
00254 
00255 
00256 //              set_style_myplain() ;
00257                 TCanvas *c1 = new TCanvas("c1","Test MICROROC");//,0,0,1024,1024);
00258 //              TApplication theApp("App", &argc, argv);
00259 //              hscurve->Draw("AP");
00260 //              theApp.Run() ;
00261 #endif
00262 //      for (int j=0;j<64;j++)
00263 //      {       
00264                 int j = chan ;
00265                 name  = Form("fermi_%i",j);
00266                 fermi[j] = new TF1(name,func,0,1024);
00267                 //fermi[j]->SetLineColor(j+1);
00268                 fermi[j]->SetLineColor(kRed);
00269                 fermi[j]->SetLineWidth(2);
00270                 /*If SCurve exist*/
00271                 if (hscurve->GetEntries()!=0)
00272                 {       /*Find threshold min tmin*/
00273                         for (int b=0;b<1023;b++)
00274                         {       content = hscurve->GetBinContent(b+1);
00275                                 if (content != 0)
00276                                 {       tmin = b ;
00277                                         b = 1023 ;
00278                                         scurve_max = content;
00279                                 }
00280                         }
00281                         /*Find threshold at which scurve = 0 tmax*/
00282                         for (int b=1024;b>tmin;b--)
00283                         {       if (hscurve->GetBinContent(b+1) != 0)
00284                                 {       tmax = b+1 ;
00285                                         b = tmin ;
00286                                 }
00287                         }
00288                         /*Find threshold at which scurve reaches the max*/
00289                         for (int b=tmax;b>=0;b--)
00290                         {       if (hscurve->GetBinContent(b+1) == scurve_max)
00291                                 {       tmed = b+1 ;
00292                                         b = -1 ;
00293                                 }
00294                         }
00295                         scurve_inflex = 0.5*(tmax+tmed);
00296                         delta = tmax - tmed;
00297   
00298                         /*If SCurve is a step function, do not fit parameters*/
00299                         if (delta==0)
00300                         {       inflex = scurve_inflex;
00301                                 sigma = 0.;
00302                                 inflex_err = 0.;
00303                                 sigma_err = 0.;
00304                         }
00305                         else
00306                         {       /*Set parameters*/
00307                                 fermi[j]->FixParameter(0,scurve_max);
00308                                 fermi[j]->SetParameter(1,scurve_inflex);
00309                                 fermi[j]->SetParameter(2,delta/10.);
00310                                 /*Perform fit*/
00311                                 hscurve->Fit(fermi[j],"q","",0,tmax+50);
00312                                 /*Get parameters*/
00313                                 inflex = fermi[j]->GetParameter(1);
00314                                 sigma = sqrt(pow(fermi[j]->GetParameter(2),2));
00315                                 inflex_err = sqrt(pow(fermi[j]->GetParError(1),2));
00316                                 sigma_err = sqrt(pow(fermi[j]->GetParError(2),2));
00317                                 chi2 = fermi[j]->GetChisquare();
00318                                 ndf = fermi[j]->GetNDF();
00319                         }
00320 
00321                         hc_inflex_sigma->SetBinContent(j+1,inflex);
00322                         hc_inflex_sigma->SetBinError(j+1,sigma);
00323 
00324                         hc_inflex->SetBinContent(j+1,inflex);
00325                         hc_sigma->SetBinContent(j+1,sigma);
00326 
00327                         hc_inflex_err->SetBinContent(j+1,inflex_err);
00328                         hc_sigma_err->SetBinContent(j+1,sigma_err);
00329 
00330                         hinflex->Fill(inflex);
00331                         hsigma->Fill(sigma);
00332 
00333                         #ifdef DISPLAY
00334                                 cout<<" Channel\t" << j << "\t" << inflex << "\t width\t" << sigma << endl ;
00335                         #endif
00336                 }
00337 //      }
00338         #ifdef WRITE
00339                 root_filename.ReplaceAll(".root", "");
00340                 TString outfile = root_directory + "../RESULTS/" +  "parameters_" + root_filename + ".txt";
00341                 ofstream out ;
00342                 out.open(outfile) ;
00343                 for (int j=0;j<64;j++)
00344                 {       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;
00345                 }
00346         #endif
00347         #ifdef PLOT
00348                 #ifdef INFORENO
00349                         cout << "Plotting" << endl ;
00350                 #endif
00351 
00352                 TString adj;
00353                 if (dac==0){adj = "Low";}
00354                 if (dac==1){adj = "Medium";}
00355                 if (dac==2){adj = "High";}
00356                 c1->Connect("Closed()", "TApplication", &theApp, "Terminate()");
00357                 //c1->Divide(8,8);
00358                 c1->cd();
00359                 hscurve->Draw("");
00360                 TString htitle = adj + " threshold";
00361 /*              for (int j=0;j<64;j++)
00362                 {       //c0->cd(j/8,j%8);
00363                         c1->cd(j+1);
00364                         hscurve->SetLineWidth(2);
00365                         hscurve->Draw("");
00366                         c1->Update() ;
00367                 }               
00368 */              
00369                 hscurve->GetXaxis()->SetTitle("Threshold [DAC]");
00370                 hscurve->GetYaxis()->SetTitle("N");
00371                 TPaveStats *p1 = (TPaveStats*) hscurve->GetListOfFunctions()->FindObject("stats");
00372                 
00373                 
00374                 c1->Update() ;
00375                 #ifdef INFORENO
00376                         cout << "Update..." << endl ;
00377                 #endif
00378                 theApp.Run() ;
00379                 #ifdef INFORENO
00380                         cout << "Do you see something ?" << endl ;
00381                 #endif
00382         #endif
00383     return 0 ;
00384 }
00385 
00386 
00387 
00388 //**********************************************************************************************
00389 //**********************************************************************************************
00390 
00391 
00392 void usage(char *progname)
00393 {
00394         printf("Usage : %s directory rootfile [0|1|2] channel\n", progname);
00395         printf("        Where [0|1|2] is the dac number\n");
00396         printf("               channel is the channel number (0 to 63)\n");
00397 }
00398 
00399 //**********************************************************************************************
00400 //**********************************************************************************************
00401 #ifdef PLOT
00402 void set_style_myplain()
00403 {
00404         style_myplain->SetCanvasBorderMode(0);
00405         style_myplain->SetCanvasColor(0);
00406         style_myplain->SetDrawBorder(0);
00407         style_myplain->SetPadBorderMode(0);
00408         style_myplain->SetPadColor(10);
00409         style_myplain->SetFrameLineColor(1);
00410         style_myplain->SetFrameFillColor(5);
00411         style_myplain->SetFrameFillStyle(0);
00412         style_myplain->SetFrameBorderMode(0);
00413 
00414         //style_myplain->SetFrameLineWidth(5); //test for setframelinecolor !
00415 
00416         style_myplain->SetLegendBorderSize(1);
00417 
00418         style_myplain->SetTextColor(231);
00419 
00420 //      style_myplain->SetFillColor(0); // If enable : impossible to fill histo with color !!! see note in the end !
00421         style_myplain->SetLineColor(231);
00422 
00423         style_myplain->SetStatColor(10);
00424         style_myplain->SetStatTextColor(1);
00425         style_myplain->SetStatBorderSize(1) ;
00426         style_myplain->SetStatX(0.9);
00427         style_myplain->SetStatY(0.9);
00428         style_myplain->SetStatW(0.1);
00429         style_myplain->SetStatH(0.1);
00430         style_myplain->SetOptStat("");//"emr" by default
00431         style_myplain->SetOptFit(0111);
00432 
00433         style_myplain->SetTitleTextColor(1);
00434         style_myplain->SetTitleBorderSize(0);
00435         style_myplain->SetTitleAlign(23);
00436         style_myplain->SetTitleX(0.5);
00437 
00438         style_myplain->SetAxisColor(1,"x");     // bleu
00439         style_myplain->SetAxisColor(1,"y");
00440         style_myplain->SetAxisColor(1,"z");
00441         style_myplain->SetLabelColor(1,"x");
00442         style_myplain->SetLabelColor(1,"x");
00443         style_myplain->SetLabelColor(1,"y");
00444         style_myplain->SetLabelColor(1,"z");
00445         style_myplain->SetLabelSize(0.03,"x");
00446         style_myplain->SetLabelSize(0.03,"y");
00447         style_myplain->SetLabelSize(0.03,"z");
00448         style_myplain->SetTitleColor(1,"x");
00449         style_myplain->SetTitleColor(1,"y");
00450         style_myplain->SetTitleColor(1,"z");
00451         style_myplain->SetTitleFillColor(10);
00452 
00453         style_myplain->SetHistFillColor(10);
00454         //style_myplain->SetHistFillStyle(0);
00455         style_myplain->SetHistFillStyle(1001);
00456         style_myplain->SetHistLineColor(1);                     // 231 !!!
00457 
00458         style_myplain->SetFuncColor(kOrange+10);        // dark orange
00459         style_myplain->SetPalette(1);                           // temp�rature like
00460         
00461 //      style_myplain->SetPaperSize(100,100);  // for Hi Res plots only !
00462 
00463         style_myplain->cd();
00464 
00465         gROOT->SetStyle("myplain"); 
00466 
00467         gROOT->ForceStyle(true);
00468 
00469 
00470 
00471 // frame style not applied ! why ????????
00472 // do it manually : right click on frame then "use current style"
00473 }
00474 #endif
/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/Renaud/fit_scurve_chip.cpp
