/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/root/DiracScurve.cpp

Go to the documentation of this file.

/* @version $Revision: 1566 $ * @modifiedby $Author: jacquem $ * @lastmodified $Date: 2012-03-08 14:00:10 +0100 (Thu, 08 Mar 2012) $ */
#include <TROOT.h>
#include <TCanvas.h>
#include <TTree.h>
#include <TFile.h>
#include <TSystem.h>
#include <TColor.h>
#include <TStyle.h>
#include <TAxis.h>
#include <TString.h>
#include <TGraph.h>
#include <TLegend.h>
#include <TCanvas.h>
#include <TApplication.h>
#include <TList.h>
#include <TH2F.h>

#include <sstream>
#include <fstream>

#include "Log.hh"
#include "Run.hh"
#include "Chamber.hh"
#include "Event.hh"
#include "Detector.hh"
#include "Centaure.hh"
#include "DiracLabview.hh"
#include "DiracChamber1.hh"
#include "GassiplexChamber1.hh"
#include "GassiplexChamber4.hh"
#include "Hardroc1Chamber.hh"
#include "XMLTool.hh"



string itos(int i)
{
        stringstream s;
        s << i;
        return s.str();
}


typedef struct{
        Short_t fChipID;
        Short_t fRunNb;
        UInt_t fChannelMin, fChannelStep, fChannelMax;
        Double_t fAmplitudeMin, fAmplitudeStep, fAmplitudeMax;
        UInt_t fDacMinA[3],fDacStepA[3],fDacMaxA[3];
        UInt_t fDacMinB[3],fDacStepB[3],fDacMaxB[3];
        UInt_t fNbDump;
}THeader;

using namespace std;


/*****************************************************************************************
  Lecture des arguments
*****************************************************************************************/
int main(int argc, char **argv){  
  FILELog::ReportingLevel() = FILELog::FromString(INFO);

  //----- Control usage and option
  if ( argc != 4  ) {
   FILE_LOG(logERROR)  << "usage: DiracScurve streeFile header rootFile" << endl;
   exit(1);
  }

  string steerName;
  steerName.assign(argv[1]);

  string headerName;
  headerName.assign(argv[2]);

  string OutputName;
  OutputName.assign(argv[3]);



/*****************************************************************************************
  Ouverture des fichiers
*****************************************************************************************/
// test if steerfile exist
        FILE *file = fopen(steerName.c_str(), "r");
        if(file == NULL){
                FILE_LOG(logERROR)  << "Steer file ["<< steerName.c_str() << "] does not exist" << endl;
                exit(0);
                fclose(file);
        }

// test if header file exist
        fstream HeaderFile;
        HeaderFile.open(headerName.c_str(),ifstream::in) ;

// create Run, set detector and steerDesc from steer XML file
        SteerDesc steerDesc;
        XMLTool xml;
        xml.parse(steerDesc, steerName);
        Detector detector;
        Run run(detector);
        if  ( !steerDesc.setRun(run))
        {
          exit(-1);
        }
        fclose(file);
        string rootName = "default.root";
        if ( steerDesc.outputFile.length() != 0)  { rootName = steerDesc.outputFile; }
/*****************************************************************************************
  Lecture du header
*****************************************************************************************/
        char strReadLine[256];
        THeader MyHeader;
        Bool_t bFlagRun = true;
        while((HeaderFile.eof()!=1) && bFlagRun){
                if (HeaderFile.bad()==1)        cerr << "Error reading file !" << endl ;        
                HeaderFile.getline(strReadLine,256) ;
                // Read Chip Identifier
                if (!strncmp(strReadLine,"Chip dirac N = ",15)){
                        sscanf(&strReadLine[15],"%hd",&(MyHeader.fChipID));
                }
                // Read Run Number
                if (!strncmp(strReadLine,"nbrun = ",18)){
                        sscanf(&strReadLine[8],"%hd",&(MyHeader.fRunNb));
                }
                // Read Channel Interval
                if (!strncmp(strReadLine,"CH_min:CH_step:CH_max = ",24)){
                        sscanf(&strReadLine[24],"%d %d %d",&(MyHeader.fChannelMin),&(MyHeader.fChannelStep),&(MyHeader.fChannelMax));
                }
                // Read Amplitude Interval
                if (!strncmp(strReadLine,"Q_min:Q_step:Q_max = ",21)){
                        sscanf(&strReadLine[21],"%lE %lE %lE",&(MyHeader.fAmplitudeMin),&(MyHeader.fAmplitudeStep),&(MyHeader.fAmplitudeMax));
                }
                // Read DAC_A Interval High
                if (!strncmp(strReadLine,"HI_A_min:HI_A_step:HI_A_max = ",30)){
                        sscanf(&strReadLine[30],"%d %d %d",&(MyHeader.fDacMinA[2]),&(MyHeader.fDacStepA[2]),&(MyHeader.fDacMaxA[2]));
                }
                // Read DAC_A Interval Mid
                if (!strncmp(strReadLine,"MID_A_min:MID_A_step:MID_A_max = ",33)){
                        sscanf(&strReadLine[33],"%d %d %d",&(MyHeader.fDacMinA[1]),&(MyHeader.fDacStepA[1]),&(MyHeader.fDacMaxA[1]));
                }
                // Read DAC_A Interval Low
                if (!strncmp(strReadLine,"LO_A_min:LO_A_step:LO_A_max = ",30)){
                        sscanf(&strReadLine[30],"%d %d %d",&(MyHeader.fDacMinA[0]),&(MyHeader.fDacStepA[0]),&(MyHeader.fDacMaxA[0]));
                }
                // Read DAC_B Interval High
                if (!strncmp(strReadLine,"HI_B_min:HI_B_step:HI_B_max = ",30)){
                        sscanf(&strReadLine[30],"%d %d %d",&(MyHeader.fDacMinB[2]),&(MyHeader.fDacStepB[2]),&(MyHeader.fDacMaxB[2]));
                }
                // Read DAC_B Interval Mid
                if (!strncmp(strReadLine,"MID_B_min:MID_B_step:MID_B_max = ",33)){
                        sscanf(&strReadLine[33],"%d %d %d",&(MyHeader.fDacMinB[1]),&(MyHeader.fDacStepB[1]),&(MyHeader.fDacMaxB[1]));
                }
                // Read DAC_B Interval Low
                if (!strncmp(strReadLine,"LO_B_min:LO_B_step:LO_B_max = ",30)){
                        sscanf(&strReadLine[30],"%d %d %d",&(MyHeader.fDacMinB[0]),&(MyHeader.fDacStepB[0]),&(MyHeader.fDacMaxB[0]));           
                }
                // Read Number of Dump
                if (!strncmp(strReadLine,"nb acquisitions = ",18)){
                        sscanf(&strReadLine[18],"%d",&(MyHeader.fNbDump));
                }
        }

/*****************************************************************************************
  Ouverture du fichier root
*****************************************************************************************/
        gSystem->Load("libMicro.so");
        TFile *f = new TFile(rootName.c_str());
        if ( f->IsZombie()) {
                FILE_LOG(logERROR) << "Root file ["<< rootName.c_str() << "] does not exist" << endl;
                exit(-1);
        }
        TTree *tree = (TTree *) gROOT->FindObject("t1");
        Int_t iNbEntries=tree->GetEntries();
        FILE_LOG(logINFO) << "Number of reoncstruct events: " << iNbEntries << endl;
        MTEvent *evt = new MTEvent();
        TBranch *branch= tree->GetBranch("MTEvent");
        branch->SetAddress(&evt);

/*****************************************************************************************
  Vérification du nombre d'évenements
*****************************************************************************************/
        Int_t iNbEvent=0;
        for (Double_t IndAmplitude = MyHeader.fAmplitudeMin;IndAmplitude <= (MyHeader.fAmplitudeMax+MyHeader.fAmplitudeStep/2);IndAmplitude += MyHeader.fAmplitudeStep){
                FILE_LOG(logINFO) << "Amplitude : " << IndAmplitude << endl;

                for (Int_t IndChannel = MyHeader.fChannelMin;IndChannel <= MyHeader.fChannelMax;IndChannel += MyHeader.fChannelStep){
                        FILE_LOG(logINFO) << "\tChannel : " << IndChannel << endl;
                        for (Int_t IndDAC = 0;IndDAC < 3;IndDAC += 1){ // Low, Mid and High DAC
                                Int_t fDacMin,fDacStep,fDacMax;
                                if (IndChannel <32){ fDacMin=MyHeader.fDacMinA[IndDAC];fDacMax=MyHeader.fDacMaxA[IndDAC];fDacStep=MyHeader.fDacStepA[IndDAC]; }
                                                                else   {        fDacMin=MyHeader.fDacMinB[IndDAC];fDacMax=MyHeader.fDacMaxB[IndDAC];fDacStep=MyHeader.fDacStepB[IndDAC]; }
                                for (Int_t IndDac = fDacMin;IndDac <= fDacMax;IndDac += fDacStep){
                                        iNbEvent += 8*MyHeader.fNbDump;
                                }
                        }
                }
        }
        FILE_LOG(logINFO) << "Nombre d'événements : " << iNbEvent << endl;

        if (iNbEvent == iNbEntries){

/**************************************************************************************
  PREPARE HISTOGRAMS
*****************************************************************************************/
                Double_t qMin = MyHeader.fAmplitudeMin,qMax = MyHeader.fAmplitudeMax,qStep = MyHeader.fAmplitudeStep;
                Int_t chMin = MyHeader.fChannelMin,chMax = MyHeader.fChannelMax,chStep = MyHeader.fChannelStep;
                Int_t dacMinL,dacMaxL,dacStepL,dacMinM,dacMaxM,dacStepM,dacMinH,dacMaxH,dacStepH;
                TClonesArray * hListLo = new TClonesArray ("TH2F", 4096) ;              // we compute all channel responses !
                TClonesArray * hListMi = new TClonesArray ("TH2F", 4096) ;              // we compute all channel responses !
                TClonesArray * hListHi = new TClonesArray ("TH2F", 4096) ;              // we compute all channel responses !
        
                FILE_LOG(logINFO) << "Bins Min : " << qMin << " Max :" << qMax+qStep << " Nb of bins : " << fixed << (qMax+qStep-qMin)/qStep << endl ;

                string hBaseName = "h_" ;
                string hName ;
                string hBaseTitle = "scurve_" ;
                string hTitle ;
                Int_t index ;
                // There are 64 histo for each stimd channel...
                // .... and 64(chmax-chmin)/chstep stimd channel
                for (Int_t stimdchan=chMin;stimdchan<=chMax;stimdchan+=chStep)
                {       for (Int_t chan=1;chan<=64;chan++)
                        {       index = 64*(stimdchan)+(chan-1) ;
                

                                if (chan ==1)   { // Bank A
                                        dacMinL=MyHeader.fDacMinA[0];dacMaxL=MyHeader.fDacMaxA[0];dacStepL=MyHeader.fDacStepA[0]; 
                                        dacMinM=MyHeader.fDacMinA[1];dacMaxM=MyHeader.fDacMaxA[1];dacStepM=MyHeader.fDacStepA[1];
                                        dacMinH=MyHeader.fDacMinA[2];dacMaxH=MyHeader.fDacMaxA[2];dacStepH=MyHeader.fDacStepA[2];
                                }
                                if (chan ==33)  { // Bank B
                                        dacMinL=MyHeader.fDacMinB[0];dacMaxL=MyHeader.fDacMaxB[0];dacStepL=MyHeader.fDacStepB[0]; 
                                        dacMinM=MyHeader.fDacMinB[1];dacMaxM=MyHeader.fDacMaxB[1];dacStepM=MyHeader.fDacStepB[1];
                                        dacMinH=MyHeader.fDacMinB[2];dacMaxH=MyHeader.fDacMaxB[2];dacStepH=MyHeader.fDacStepB[2];
                                }

                                hName = hBaseName + "S" + itos(stimdchan+1) + "_H" + itos(chan) + "l";
                                hTitle = hBaseTitle + "S" + itos(stimdchan+1) + "_H" + itos(chan) + "l" ;
                                hListLo->New(index) ;
                                ( (TH2F *) hListLo->At(index) )->SetNameTitle(hName.c_str(), hTitle.c_str()) ;
                                ( (TH2F *) hListLo->At(index) )->SetBins( (dacMaxL+1-dacMinL)/dacStepL, dacMinL, dacMaxL+1, static_cast <Int_t> ( ((qMax+qStep-qMin)/qStep) + 0.5 ), qMin, qMax+qStep) ;

                                hName = hBaseName + "S" + itos(stimdchan+1) + "_H" + itos(chan) + "m";
                                hTitle = hBaseTitle + "S" + itos(stimdchan+1) + "_H" + itos(chan) + "m" ;
                                hListMi->New(index) ;
                                ((TH2F *) hListMi->At(index) )->SetNameTitle(hName.c_str(), hTitle.c_str()) ;
                                ((TH2F *) hListMi->At(index) )->SetBins( (dacMaxM+1-dacMinM)/dacStepM, dacMinM, dacMaxM+1, static_cast <Int_t> ( ((qMax+qStep-qMin)/qStep) + 0.5 ), qMin, qMax+qStep) ;

                                hName = hBaseName + "S" + itos(stimdchan+1) + "_H" + itos(chan) + "h";
                                hTitle = hBaseTitle + "S" + itos(stimdchan+1) + "_H" + itos(chan) + "h" ;
                                hListHi->New(index) ;
                                ((TH2F *) hListHi->At(index) )->SetNameTitle(hName.c_str(), hTitle.c_str()) ;
                                ((TH2F *) hListHi->At(index) )->SetBins( (dacMaxH+1-dacMinH)/dacStepH, dacMinH, dacMaxH+1, static_cast <Int_t> ( ((qMax+qStep-qMin)/qStep) + 0.5 ), qMin, qMax+qStep) ;
                        }
                }
                FILE_LOG(logINFO) << "Init ok" << endl ;
/*****************************************************************************************
  Calcul des Scurves
*****************************************************************************************/
                Int_t currHitValue, currHitChannel;
                Int_t iEvent=0;
                for (Double_t IndAmplitude = MyHeader.fAmplitudeMin;IndAmplitude <= (MyHeader.fAmplitudeMax+MyHeader.fAmplitudeStep/2);IndAmplitude += MyHeader.fAmplitudeStep){

                        for (Int_t IndChannel = MyHeader.fChannelMin;IndChannel <= MyHeader.fChannelMax;IndChannel += MyHeader.fChannelStep){
                                for (Int_t IndDAC = 0;IndDAC < 3;IndDAC += 1){ // Low, Mid and High DAC
                                        Int_t fDacMin,fDacStep,fDacMax;
                                        if (IndChannel <32){ fDacMin=MyHeader.fDacMinA[IndDAC];fDacMax=MyHeader.fDacMaxA[IndDAC];fDacStep=MyHeader.fDacStepA[IndDAC]; }
                                                                        else   {        fDacMin=MyHeader.fDacMinB[IndDAC];fDacMax=MyHeader.fDacMaxB[IndDAC];fDacStep=MyHeader.fDacStepB[IndDAC]; }
                                        for (Int_t IndDac = fDacMin;IndDac <= fDacMax;IndDac += fDacStep){
                                                for (Int_t IndEvent = 0;IndEvent < 8*MyHeader.fNbDump;IndEvent ++){
                                                        tree->GetEntry(iEvent);
                                                        Int_t NC=evt->GetNchannel();
                                                        for(int i=0;i<NC;i++){

                                                                MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(i);
                                        
                                                                currHitValue = channel->GetDigitalValue() ;                     // ...Xtract data from channel...
                                                                currHitChannel = channel->GetHardId() ;         // ... Xtract channel number...
                                                                index = 64*IndChannel+currHitChannel ;          
                                                                switch(currHitValue){   
                                                                        case 3:
                                                                                if (IndDAC == 2){
                                                                                        ( (TH2F *) hListHi->At(index) )->Fill(IndDac, IndAmplitude) ;
                                                                                        //FILE_LOG(logINFO) << "Amplitude 3 : " << IndChannel << " / " << currHitChannel << endl;
                                                                                }       
                                                                                break ;
                                                                        case 2:
                                                                                if (IndDAC == 1){
                                                                                        ( (TH2F *) hListMi->At(index) )->Fill(IndDac, IndAmplitude) ;
                                                                                        //FILE_LOG(logINFO) << "Amplitude 2 : " << IndChannel << " / " << currHitChannel << endl;
                                                                                }
                                                                                break ;
                                                                        case 1:
                                                                                if (IndDAC == 0){
                                                                                        ( (TH2F *) hListLo->At(index) )->Fill(IndDac, IndAmplitude) ;
                                                                                        //FILE_LOG(logINFO) << "Amplitude 1 : " << IndChannel << " / " << currHitChannel << endl;
                                                                                }
                                                                                break ;
                                                                        default:
                                                                                cout << "ERROR in ROOT file !" << endl ;
                                                                }       
                                                        }
                                                        iEvent++;
                                                }
                                        }
                                }
                        }
                }  
/*****************************************************************************************
  NORMALIZE IN %, remove first and last bin
*****************************************************************************************/

        // Normalize to 100%
        Float_t scale = 1./(8*MyHeader.fNbDump) ;
        for (Int_t stimdchan=chMin;stimdchan<=chMax;stimdchan+=chStep)
        {       for (Int_t chan=1;chan<=64;chan++)
                {       index = 64*(stimdchan)+(chan-1) ;
                        //cout << scale ;
                        ( (TH2F *) hListLo->At(index) )->Scale(scale) ;
                        ( (TH2F *) hListMi->At(index) )->Scale(scale) ;
                        ( (TH2F *) hListHi->At(index) )->Scale(scale) ;
                }
        }

        cout << "Normalization ok" << endl ;
/*****************************************************************************************
  Enregistrement des Scurves
*****************************************************************************************/
                //FILE_LOG(logINFO) << "Events : " << iEvent << endl ;                          
                TFile *histFile = new TFile(OutputName.c_str(), "RECREATE") ;
                if (histFile->IsOpen() == 0)
                {       cerr << "Could not open data file "<< endl ;
                        exit(1) ;
                }
                histFile->cd() ;
                hListLo->Write() ;
                hListMi->Write() ;      

                hListHi->Write() ;


/*****************************************************************************************
  Affichage des Scurves
*****************************************************************************************/
        /*      argc=0;
                TApplication theApp("App", &argc, argv);                                        // see $ROOTSYS/test/hworld.cxx
                gStyle->SetPalette(1);
        //gStyle->SetCanvasColor(33);
        //gStyle->SetFrameFillColor(18);
                TCanvas *c1 = new TCanvas("c1","S-Curves") ;
                
                c1->Connect("Closed()", "TApplication", &theApp, "Terminate()");
                string name="h_S32_H32l";
                TH2F* h=(TH2F*) gROOT->FindObject(name.c_str()) ;
                c1->cd() ;
                h->Draw("lego2") ;
                theApp.Run();*/

                hListLo->Delete() ;
                hListMi->Delete() ;
                hListHi->Delete() ;
                evt->Delete() ;
                //h->Delete();
        }
   return 0;
}

---