/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/Renaud/test_microroc.cpp

Go to the documentation of this file.

/* @version $Revision: 1585 $ * @modifiedby $Author: jacquem $ * @lastmodified $Date: 2012-03-13 11:03:35 +0100 (Tue, 13 Mar 2012) $ */
#include <map>
#include "TKey.h"
#include <TROOT.h>
#include <TRint.h>
#include <TKey.h>
#include <TTree.h>
#include <TFile.h>
#include <TSystem.h>
#include <iostream>
#include <TF1.h>
#include <TH2.h>
#include <TGraph.h>
#include <sstream>
#include "Log.hh"
#include <fstream>
#include "TCanvas.h"
#include "TApplication.h"
#include "TColor.h"
#include "TStyle.h"
#include "TAxis.h"
#include "TString.h"
#include "TLegend.h"
#include "Bytes.h"



#include "root/MTRun.hh"
#include "root/MTEvent.hh"
#include "root/MTChannel.hh"
#include "root/MTMicrorocChip.hh"

#include "MicroException.hh"

//#define DISPLAY
#undef DISPLAY
//#define R_DEBUG
#undef R_DEBUG
//#define PLOT
#undef PLOT
#define WRITE

using namespace std;
void usage(char *progname) ;
#ifdef PLOT
        void set_style_myplain() ;
        TStyle *style_myplain=new TStyle("myplain","classic style");
#endif
////////////////////////////////////////////////////////////////////////////////////////////

int main(int argc, char**argv){

  /*****************************************************************/
  /* ARGUMENTS TO THE PROGRAM **************************************/
  /* 1 - DIRECTORY OF THE ROOT FILE PRODUCED BY THE RECONSTRUCTION */
  /* 2 - NAME OF THE ROOT FILE *************************************/
  /* 3 - THRESHOLD THAT IS MOVED DURING THE CALIBRATION ************/
  /*    0 -> DAC_0 (LOW THRESHOLD) *********************************/
  /*    1 -> DAC_1 (MEDIUM THRESHOLD) ******************************/
  /*    2 -> DAC_2 (HIGH THRESHOLD) ********************************/
  /*****************************************************************/
        if(argc != 4)
        {       usage(argv[0]) ;
                return 1 ;
        }
        TString root_directory = argv[1];
        TString root_filename = argv[2];
        int dac = atoi(argv[3]);

        //bool display = 0;
        #ifdef PLOT
                set_style_myplain() ;
        #endif
        TCanvas *c1 = new TCanvas("c1","Test MICROROC",0,0,1000,1000);
        TString root_file = root_directory + root_filename;
        TString scurve_file = root_directory + "scurve_" + root_filename;

//  cout<<"Wish to write "<<scurve_file<<endl;

        int nbChannel = 0;
        float nbEvt = 0;
        int threshold = 0;
        //int chipid = 0;
        int hardid = 0;
        int value = 0;
        //int content = 0;
        Double_t content = 0;
        UInt_t chipId ;
 
        TH1I * hvalue = new TH1I("hvalue","",5,0,5);
        TH1I * hnpulse_threshold = new TH1I("hnpulse_threshold","",1024,0,1024);

        TH1I * hscurve[64];
        TString name,title;
 
        for (int j=0;j<64;j++)
        {       name  = Form("hscurve_%i",j);
                title = Form("SCurve from channel %i",j);
                hscurve[j] = new TH1I(name,title,1024,0,1024);
        }
  
        TFile f(root_file);
        TIter nextkey(f.GetListOfKeys());
        TKey *key;
        while (key = (TKey*)nextkey())
        {       TTree *tree = (TTree*)key->ReadObj();                      
                MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");
                const MTChip&  chip = run->GetOneChip();
                //    cout << "------------ Chip[" << chip.GetId() << "], softId: [" << chip.GetSoftId() <<"] configuration ---------------" << endl ;

                try
                {       const MTMicrorocChip &microChip = dynamic_cast<const MTMicrorocChip &> (chip);
                        //const MTHardroc2Chip &microChip = dynamic_cast<const MTHardroc2Chip &> (chip);
                        //cout << "threshold dac 0:[" << microChip.GetThresholdDac_0()  << endl;
                        //cout << "threshold dac 1:[" << microChip.GetThresholdDac_1()  << endl;
                        //cout << "threshold dac 2:[" << microChip.GetThresholdDac_2()  << endl;
                        //chipId = chip.GetId() ;
                        //cout << "Chip identifier: " << chipId << endl ;
                        if (dac==0){threshold = microChip.GetThresholdDac_0();}
                        if (dac==1){threshold = microChip.GetThresholdDac_1();}
                        if (dac==2){threshold = microChip.GetThresholdDac_2();}
                }
                catch (...) {}

                MTEvent *evt =  new MTEvent();
                TBranch *branch= tree->GetBranch("MTEvent");
                branch->SetAddress(&evt);
                MTChannel* channel = NULL;
                nbEvt = tree->GetEntries();
//              cout << "Number of events in the tree: " << nbEvt << endl ;

                hnpulse_threshold->SetBinContent(threshold+1,nbEvt);

                for( int evtNum = 0; evtNum < nbEvt ; evtNum++)
                {       //if (evtNum%1000 == 0) {cout<<evtNum/nbEvt*100.<<" %"<<"\r"<<flush;}
                        tree->GetEntry(evtNum);
                        nbChannel = evt->GetNchannel();
                        for(int i=0;i<nbChannel ;i++)
                        {       channel = (MTChannel*)evt->GetChannels()->UncheckedAt(i);
                                chipId = channel->GetChipId();
                                hardid = channel->GetHardId();
                                value = channel->GetDigitalValue();
                                hvalue->Fill(value);
                                if (value==(dac+1))
                                {       hscurve[hardid]->Fill(threshold);
                                        #ifdef R_DEBUG
                                                cout << "---- digital value:" <<  value <<endl;
                                                cout << "---- hard id:" << hardid <<endl;
                                                //cout << "---- x:"<< channel->GetX() <<endl;
                                                //cout << "---- y:" << channel->GetY() <<endl;
                                                cout << "---- chip id:" << chipId <<endl;
                                        #endif
          
                                }
                        }
                }
        }
        
        


        //===========================================================================*
        // Now analyze scurves  

        TString t_scurve_file = root_directory + root_filename ;//+ ".root";
        //cout << t_scurve_file<< endl;

        int tmin = 0;
        int tmed = 0;
        int tmax = 0;
        int delta = 0;

//      TH1I * hscurve[64];

        TH1F * hinflex = new TH1F("hinflex","",1024,0,1024);
        TH1F * hsigma = new TH1F("hsigma","",50,-2,2);

        TH1F * hinflex_err = new TH1F("hinflex_err","",100,0,10);
        TH1F * hsigma_err = new TH1F("hsigma_err","",100,0,2);

        TH1F * hc_inflex = new TH1F("hc_inflex","",64,0,64);
        TH1F * hc_sigma = new TH1F("hc_sigma","",64,0,64);

        TH1F * hc_inflex_err = new TH1F("hc_inflex_err","",64,0,64);
        TH1F * hc_sigma_err = new TH1F("hc_sigma_err","",64,0,64);

        TH1F * hc_inflex_sigma = new TH1F("hc_inflex_sigma","",64,0,64);

//      TString name;

/*      for (int j=0;j<64;j++)
        {       name  = Form("hscurve_%i",j);
                TH1I * h = (TH1I*)tf.Get(name);
                hscurve[j] = h;
        }
*/
        TString func = "[0]/(1+exp((x-[1])/[2]))";
        TF1 * fermi[64];
        float scurve_max = 0;
        float scurve_inflex = 0;
        float inflex = 0;
        float sigma = 0;
        float inflex_err = 0;
        float sigma_err = 0;
        float content0,content1,content2;
        float chi2 = 0;
        float ndf = 0;
        
        for (int j=0;j<64;j++)
        {       name  = Form("fermi_%i",j);
                fermi[j] = new TF1(name,func,0,1024);
                //fermi[j]->SetLineColor(j+1);
                fermi[j]->SetLineColor(kRed);
                fermi[j]->SetLineWidth(2);
                /*If SCurve exist*/
                if (hscurve[j]->GetEntries()!=0)
                {       /*Find threshold min tmin*/
                        for (int b=0;b<1023;b++)
                        {       content = hscurve[j]->GetBinContent(b+1);
                                if (content != 0)
                                {       tmin = b ;
                                        b = 1023 ;
                                        scurve_max = content;
                                }
                        }
                        /*Find threshold at which scurve = 0 tmax*/
                        for (int b=1024;b>tmin;b--)
                        {       if (hscurve[j]->GetBinContent(b+1) != 0)
                                {       tmax = b+1 ;
                                        b = tmin ;
                                }
                        }
                        /*Find threshold at which scurve reaches the max*/
                        for (int b=tmax;b>=0;b--)
                        {       if (hscurve[j]->GetBinContent(b+1) == scurve_max)
                                {       tmed = b+1 ;
                                        b = -1 ;
                                }
                        }
                        scurve_inflex = 0.5*(tmax+tmed);
                        delta = tmax - tmed;
  
                        /*If SCurve is a step function, do not fit parameters*/
                        if (delta==0)
                        {       inflex = scurve_inflex;
                                sigma = 0.;
                                inflex_err = 0.;
                                sigma_err = 0.;
                        }
                        else
                        {       /*Set parameters*/
                                fermi[j]->FixParameter(0,scurve_max);
                                fermi[j]->SetParameter(1,scurve_inflex);
                                fermi[j]->SetParameter(2,delta/10.);
                                /*Perform fit*/
                                hscurve[j]->Fit(fermi[j],"q","",0,tmax+50);
                                /*Get parameters*/
                                inflex = fermi[j]->GetParameter(1);
                                sigma = sqrt(pow(fermi[j]->GetParameter(2),2));
                                inflex_err = sqrt(pow(fermi[j]->GetParError(1),2));
                                sigma_err = sqrt(pow(fermi[j]->GetParError(2),2));
                                chi2 = fermi[j]->GetChisquare();
                                ndf = fermi[j]->GetNDF();
                        }

                        hc_inflex_sigma->SetBinContent(j+1,inflex);
                        hc_inflex_sigma->SetBinError(j+1,sigma);

                        hc_inflex->SetBinContent(j+1,inflex);
                        hc_sigma->SetBinContent(j+1,sigma);

                        hc_inflex_err->SetBinContent(j+1,inflex_err);
                        hc_sigma_err->SetBinContent(j+1,sigma_err);

                        hinflex->Fill(inflex);
                        hsigma->Fill(sigma);

                        #ifdef DISPLAY
                                cout<<" Channel\t" << j << "\t" << inflex << "\t width\t" << sigma << endl ;
                        #endif
                }
        }
        
        ////////////////////////////////////////////////////////////////////////////////
        // Write Scurves
        TFile * tf = new TFile(scurve_file,"RECREATE");
        hvalue->Write();
        hnpulse_threshold->Write();
        for (int j=0;j<64;j++)
        {       hscurve[j]->Write();}
        tf->Close();
        ////////////////////////////////////////////////////////////////////////////////
        
        #ifdef WRITE
                root_filename.ReplaceAll(".root", "");
                TString outfile = root_directory + "../RESULTS/" +  "parameters_" + root_filename + ".txt";
                ofstream out ;
                out.open(outfile) ;
                for (int j=0;j<64;j++)
                {       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;
                }
        #endif
        #ifdef PLOT
                TApplication theApp("App", &argc, argv);
                TString adj;
                if (dac==0){adj = "Low";}
                if (dac==1){adj = "Medium";}
                if (dac==2){adj = "High";}
                
                c1->Divide(8,8);
                //c1->cd();
                TString htitle = adj + " threshold";
                for (int j=0;j<64;j++)
                {       //c0->cd(j/8,j%8);
                        c1->cd(j+1);
                        hscurve[j]->SetLineWidth(2);
                        hscurve[j]->Draw("");
                        c1->Update() ;
                }
                theApp.Run(kTRUE) ;
                c1->Update() ;
        #endif
    return 0 ;
}



//**********************************************************************************************
//**********************************************************************************************


void usage(char *progname)
{
        printf("Usage : %s directory rootfile [0|1|2]", progname);
        printf("        Where [0|1|2] is the dac number\n");
        printf("        scurve_rootfile.root is created\n");
}

//**********************************************************************************************
//**********************************************************************************************
#ifdef PLOT
void set_style_myplain()
{
        style_myplain->SetCanvasBorderMode(0);
        style_myplain->SetCanvasColor(0);
        style_myplain->SetDrawBorder(0);
        style_myplain->SetPadBorderMode(0);
        style_myplain->SetPadColor(10);
        style_myplain->SetFrameLineColor(1);
        style_myplain->SetFrameFillColor(5);
        style_myplain->SetFrameFillStyle(0);
        style_myplain->SetFrameBorderMode(0);

        //style_myplain->SetFrameLineWidth(5); //test for setframelinecolor !

        style_myplain->SetLegendBorderSize(1);

        style_myplain->SetTextColor(231);

//      style_myplain->SetFillColor(0); // If enable : impossible to fill histo with color !!! see note in the end !
        style_myplain->SetLineColor(231);

        style_myplain->SetStatColor(10);
        style_myplain->SetStatTextColor(1);
        style_myplain->SetStatBorderSize(1) ;
        style_myplain->SetStatX(0.9);
        style_myplain->SetStatY(0.9);
        style_myplain->SetOptStat("");//"emr" by default
        style_myplain->SetOptFit(0111);

        style_myplain->SetTitleTextColor(1);
        style_myplain->SetTitleBorderSize(0);
        style_myplain->SetTitleAlign(23);
        style_myplain->SetTitleX(0.5);

        style_myplain->SetAxisColor(1,"x");     // bleu
        style_myplain->SetAxisColor(1,"y");
        style_myplain->SetAxisColor(1,"z");
        style_myplain->SetLabelColor(1,"x");
        style_myplain->SetLabelColor(1,"x");
        style_myplain->SetLabelColor(1,"y");
        style_myplain->SetLabelColor(1,"z");
        style_myplain->SetLabelSize(0.03,"x");
        style_myplain->SetLabelSize(0.03,"y");
        style_myplain->SetLabelSize(0.03,"z");
        style_myplain->SetTitleColor(1,"x");
        style_myplain->SetTitleColor(1,"y");
        style_myplain->SetTitleColor(1,"z");
        style_myplain->SetTitleFillColor(10);

        style_myplain->SetHistFillColor(10);
        //style_myplain->SetHistFillStyle(0);
        style_myplain->SetHistFillStyle(1001);
        style_myplain->SetHistLineColor(1);                     // 231 !!!

        style_myplain->SetFuncColor(kOrange+10);        // dark orange
        style_myplain->SetPalette(1);                           // temp�rature like
        
//      style_myplain->SetPaperSize(100,100);  // for Hi Res plots only !

        style_myplain->cd();

        gROOT->SetStyle("myplain"); 

        gROOT->ForceStyle(true);



// frame style not applied ! why ????????
// do it manually : right click on frame then "use current style"
}
#endif

---