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

Go to the documentation of this file.

/* @version $Revision: 1328 $ * @modifiedby $Author: jacquem $ * @lastmodified $Date: 2011-10-03 17:04:17 +0200 (Mon, 03 Oct 2011) $ */
#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/MTTrack.hh"
#include "root/MTChamber.hh"
#include "root/MTDetector.hh"
#include "root/MTBoard.hh"
#include "root/MTChip.hh"
#include "root/MTChannel.hh"
#include "root/MTEvent.hh"
#include "root/MTDif.hh"
#include "root/MTHardroc1Chip.hh"
#include "root/MTHardroc2Chip.hh"
#include "root/MTDiracChip.hh"
#include "root/MTMicrorocChip.hh"

#include "MicroException.hh"

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

//#define CHIP_ID 32
#define CHIP_ID 33

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



// src/analyse/root/example.cpp

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]);

        #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;
        int nbChannel = 0;
        float nbEvt = 0;
        int threshold = 0;
        //int chipid = 0;
        int hardid = 0;
        int value = 0;
        int serial = 0;
        //int content = 0;
        Double_t content = 0;
        UInt_t chipId ;
    int dac_min = 1025;
        int dac_max = 0;
        
        TH1I * hvalue = new TH1I("hvalue","",5,0,5);
        TH1I * hnpulse_threshold = new TH1I("hnpulse_threshold","",1024,0,1024);

        TH1I * hscurve[64] ;
        TH1I * hscurve_deriv[64] ;
        TString name, title, name_deriv, title_deriv ;
        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);
                name_deriv  = Form("deriv_%i",j);
                title_deriv = Form("Derivative from channel %i",j);
                hscurve_deriv[j] = new TH1I(name_deriv,title_deriv,1024,0,1024);
        }
  
        TFile f(root_file);
        
        int ttiter=0 ;
        
        // from example file !
        TIter nextkey(f.GetListOfKeys());
        TKey *key;
  while (key = (TKey*)nextkey())
  {
  TTree *tree = (TTree*)key->ReadObj();                
  MTRun *run=NULL;
  try 
        {
    run=dynamic_cast<MTRun *>(tree->GetUserInfo()->First());
//    cout << "---------------------------- NEW TTREE ---------------------------"<< endl;
//    ttiter++ ;
//    cout << ttiter << endl ;
//    nbEvt =0;
//    nbHit = 0;

    bool showConfig = true;
    if ( run != NULL && showConfig) 
    {
        //run->Info(); 
                        MTDetector* det = run->GetDetector();
                        // loop over chamber config
                                
                        const std::map<UInt_t,MTChamber*>& chambers = det->GetChambers();
                        for( std::map<UInt_t,MTChamber*>::const_iterator ii=chambers.begin(); ii!=chambers.end(); ++ii)
                        {
                                const MTChamber& chamber= *((*ii).second);
                                //cout << "chamber id: " << chamber.GetId() << endl;
                                //cout << "chamber type: " << chamber.GetType() << endl;

                                const std::map<UInt_t,MTBoard*>& boards = chamber.GetBoards();
                                for( std::map<UInt_t,MTBoard*>::const_iterator ii=boards.begin(); ii!=boards.end(); ++ii)
                                {
                                        const MTBoard& board= *((*ii).second);
                                        //board.Info();

                                        //cout << "\t\t\t**** Chips configuration ****" << endl ;

                                        map<UInt_t,MTChip*> chips = board.GetChips();
                                        for( map<UInt_t,MTChip*>::const_iterator iiChip=chips.begin(); iiChip!=chips.end(); ++iiChip)
                                        {
                                                const MTChip& chip  = *((*iiChip).second);
                                                if ( chip.IsConfigured() )
                                                {
                                                //cout <<  "   -- New Chip configuration:" <<  chip.GetSoftId().toString() << " softId[" << chip.GetSoftId().GetValue() <<"], serial Number[" << chip.GetSerialNumber() << "]"<< endl <<endl  <<endl ;
                                                        if (chip.GetSerialNumber() == 14) { // serial 14 = id 32 -> see xml file
                                                        try{
                                                        
                                                        const MTMicrorocChip &microChip = dynamic_cast<const MTMicrorocChip &> (chip);

                                                        if (dac==0){threshold = microChip.GetThresholdDac_0();}
                                                        if (dac==1){threshold = microChip.GetThresholdDac_1();}
                                                        if (dac==2){threshold = microChip.GetThresholdDac_2();}
                                                        //cout << "Chip: " << chip.GetSerialNumber() << " SoftId: " << chip.GetSoftId().toString() << " Threshold:" << threshold << " SoftId:" << endl ;
                                                        }
                                                        catch (...) {}
                                                        
                                                if (threshold<dac_min){dac_min = threshold;}
                                        if (threshold>dac_max){dac_max = threshold;}

                                                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);
//                                                      if (chip.GetSerialNumber() == 14) { 
                                                        //ttiter++ ;
                                                        //cout << "procees " << ttiter << endl ;
                                                        nbChannel = evt->GetNchannel();
                                                        for(int i=0;i<nbChannel ;i++)
                                                        {       channel = (MTChannel*)evt->GetChannels()->UncheckedAt(i);
                                                                //if (channel->getChip()->GetSerialNumber() == 14)
                                                                //cout << channel->>GetSerialNumber() << endl ;
                                                                if (chipId = channel->GetChipId() == CHIP_ID) // numéro du chip
                                                                {hardid = channel->GetHardId(); // numéro de la voie
                                                                value = channel->GetDigitalValue();
                                                                hvalue->Fill(value);
                                                                if (value==(dac+1)) // ie 1 for low, 2 for med and 3 for high ....
                                                                {       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
                                                                }}
                                                        }
                                                }
        
                                                        
                                                        
                                                        
                                                
} // fin serial == 14   
                                        }
                                }
                        }
                }
        }
        
        

        
        
        
        
        }
        catch (MicroException e) {}             
        }       
                                                

//      cout << "ttree iter " << ttiter << endl ;

//cout << "           " << endl ;

        
        //===========================================================================*
        // 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;
        }*/

        float scurve_max = 0;
        float scurve_inflex = 0;
        float inflex = 0;
        float sigma = 0;
        float inflex_err = 0;
        float sigma_err = 0;
        float chi2 = 0;
        float ndf = 0;
        Double_t deriv = 0;
        
        bool start_100 ; // test if a channel starts at 100% else warning
        bool has_flipped ; // test if each channel flips from 100% to 0% else, warning message
        Int_t noisy;    // test if each channel is noisy, ie has entry beyond flip point
        Int_t noise_limit = 2;  // number of entry accepted
        Int_t noise_window = 10 ;// number of bin defore dac_max to verify noise
        Int_t flip_limit = 30 ;
          
        for (int j=0;j<64;j++)
        {       noisy = 0 ;
                for (int k=dac_min+1;k<=dac_max;k++) // ATTN : bin = dac + 1 !!
                {       deriv = ( hscurve[j]->GetBinContent(k)) - (hscurve[j]->GetBinContent(k+1));
                        //cout << j << " " << k << " " << hscurve[j]->GetBinContent(k) << " " << hscurve[j]->GetBinContent(k+1) << " " << deriv << endl ;
                        // test if each channel has flipped
                        if  (k==dac_min+1) // first loop
                        {       if (hscurve[j]->GetBinContent(k)>=100)
                                        start_100 = 1 ;
                                else start_100 = 0 ;
                        }
                        
                        if  (k>(dac_max - noise_window)) // test noisy channel
                        {       if (hscurve[j]->GetBinContent(k)>=noise_limit)
                                        noisy++ ;
                                //cout << "Channel " << j << " dac " << k << " noise " << noisy << endl ;
                        }
                        
                        
                        
                        if  (k==dac_max-1) // last loop has flipped ?
                        {       if (hscurve[j]->GetBinContent(k)<=flip_limit && start_100 == 1 )
                                {       has_flipped = 1 ;
                                        //if (hscurve[j]->GetBinContent(k)>0)
                                        //      noisy++;// = 0 ;
                                        //else noisy = 1 ;
                                }
                                else has_flipped = 0 ;
                        }                       
                        if (deriv > 0.)
                        {       hscurve_deriv[j]->SetBinContent(k,deriv);
                                hscurve_deriv[j]->SetEntries((hscurve_deriv[j]->GetEntries())+deriv);
                                //hscurve_deriv[j]->Fill(k,deriv);
                        }
                }
                if (has_flipped == 0)
                {       cout << "****** " << argv[0] << ": Warning! Chip " << serial << " Channel " << j << " has not flipped" ;
                        if (start_100 == 1) cout << " (stuck to 100 %)" << endl ;
                        else cout << " (stuck to 0 %)" << endl ;
                }
                if (noisy==noise_window) cout << "****** " << argv[0] << ": Warning! Chip " << serial << " Channel " << j << " may be noisy, check root scurve file !" << endl ;
                
         }

// FIT
        TF1 * gaussfunc[64];
        for (int j=0;j<64;j++)
        {       name = Form("g_%i",j);
                gaussfunc[j] = new TF1(name,"gaus",0,1024);
        }
        
        for (int j=0;j<64;j++)
        {       //gaussfunc[j]->SetParLimits(1, dac_min, dac_max) ;
                gaussfunc[j]->SetParameter(0, hscurve_deriv[j]->GetMaximumBin());
                gaussfunc[j]->SetParLimits(0, 1, 100);
                gaussfunc[j]->SetParameter(1, hscurve_deriv[j]->GetMean() ); //GetMaximumBin()
                //gaussfunc[j]->FixParameter(1, hscurve_deriv[j]->GetMean() );
                gaussfunc[j]->SetParLimits(1, dac_min, dac_max) ; //hscurve_deriv[j]->GetMaximumBin()-10, hscurve_deriv[j]->GetMaximumBin() + 10) ;
                gaussfunc[j]->SetParameter(2, hscurve_deriv[j]->GetRMS() );
                gaussfunc[j]->SetParLimits(2, 0.01, hscurve_deriv[j]->GetRMS()) ;
//              if (j==48)
//              hscurve_deriv[j]->Fit(gaussfunc[j],"BV", "",hscurve_deriv[j]->GetMean()-15, hscurve_deriv[j]->GetMean() + 15 );
//              else
                                hscurve_deriv[j]->Fit(gaussfunc[j],"WBQ", "",hscurve_deriv[j]->GetMean()-20, hscurve_deriv[j]->GetMean() + 20 );
                //hscurve_deriv[j]->Fit(gaussfunc[j],"Q");
                inflex = gaussfunc[j]->GetParameter(1);
                if (inflex<=dac_min || inflex>=dac_max)
                {       cout << "****** " << argv[0] << ": Warning! Chip " << serial << " Channel " << j << ", fit may be not correct ! Inflexion: " << inflex << " Dac min: " << dac_min << " Dac max: " << dac_max << endl ;
                }
                sigma = gaussfunc[j]->GetParameter(2);
                inflex_err = gaussfunc[j]->GetParError(1);
                sigma_err =gaussfunc[j]->GetParError(2);
                
                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);

                //cout << "Channel " << j << " init at " << hscurve_deriv[j]->GetMean() << " " << gaussfunc[j]->GetParameter(1) << endl ;
                // cout << j << " " << inflex << " " << sigma << endl ;
        }


        //////////////////////////////////////////////////////////////////////////////////////////
        // Write Scurves
        //////////////////////////////////////////////////////////////////////////////////////////
        TFile * tf = new TFile(scurve_file,"RECREATE");
        hvalue->Write();
        hnpulse_threshold->Write();
        for (int j=0;j<64;j++)
        {       hscurve[j]->Write();
                hscurve_deriv[j]->Write();
                //gaussfunc[j]->Write();
        }
        tf->Close();
        
        
        #ifdef WRITE
                root_filename.ReplaceAll(".root", "");
                //TString outfile = root_directory + "../RESULTS/" +  "parameters_" + root_filename + ".txt";
                TString outfile = root_directory + "../RESULTS/" +  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

---