#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 <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"Include dependency graph for DiracScurve.cpp:
Go to the source code of this file.
Classes | |
| struct | THeader |
Functions | |
| string | itos (int i) |
| int | main (int argc, char **argv) |
| string itos | ( | int | i | ) |
| int main | ( | int | argc, | |
| char ** | argv | |||
| ) |
Definition at line 60 of file DiracScurve.cpp.
00060 { 00061 FILELog::ReportingLevel() = FILELog::FromString(INFO); 00062 00063 //----- Control usage and option 00064 if ( argc != 4 ) { 00065 FILE_LOG(logERROR) << "usage: DiracScurve streeFile header rootFile" << endl; 00066 exit(1); 00067 } 00068 00069 string steerName; 00070 steerName.assign(argv[1]); 00071 00072 string headerName; 00073 headerName.assign(argv[2]); 00074 00075 string OutputName; 00076 OutputName.assign(argv[3]); 00077 00078 00079 00080 /***************************************************************************************** 00081 Ouverture des fichiers 00082 *****************************************************************************************/ 00083 // test if steerfile exist 00084 FILE *file = fopen(steerName.c_str(), "r"); 00085 if(file == NULL){ 00086 FILE_LOG(logERROR) << "Steer file ["<< steerName.c_str() << "] does not exist" << endl; 00087 exit(0); 00088 fclose(file); 00089 } 00090 00091 // test if header file exist 00092 fstream HeaderFile; 00093 HeaderFile.open(headerName.c_str(),ifstream::in) ; 00094 00095 // create Run, set detector and steerDesc from steer XML file 00096 SteerDesc steerDesc; 00097 XMLTool xml; 00098 xml.parse(steerDesc, steerName); 00099 Detector detector; 00100 Run run(detector); 00101 if ( !steerDesc.setRun(run)) 00102 { 00103 exit(-1); 00104 } 00105 fclose(file); 00106 string rootName = "default.root"; 00107 if ( steerDesc.outputFile.length() != 0) { rootName = steerDesc.outputFile; } 00108 /***************************************************************************************** 00109 Lecture du header 00110 *****************************************************************************************/ 00111 char strReadLine[256]; 00112 THeader MyHeader; 00113 Bool_t bFlagRun = true; 00114 while((HeaderFile.eof()!=1) && bFlagRun){ 00115 if (HeaderFile.bad()==1) cerr << "Error reading file !" << endl ; 00116 HeaderFile.getline(strReadLine,256) ; 00117 // Read Chip Identifier 00118 if (!strncmp(strReadLine,"Chip dirac N = ",15)){ 00119 sscanf(&strReadLine[15],"%hd",&(MyHeader.fChipID)); 00120 } 00121 // Read Run Number 00122 if (!strncmp(strReadLine,"nbrun = ",18)){ 00123 sscanf(&strReadLine[8],"%hd",&(MyHeader.fRunNb)); 00124 } 00125 // Read Channel Interval 00126 if (!strncmp(strReadLine,"CH_min:CH_step:CH_max = ",24)){ 00127 sscanf(&strReadLine[24],"%d %d %d",&(MyHeader.fChannelMin),&(MyHeader.fChannelStep),&(MyHeader.fChannelMax)); 00128 } 00129 // Read Amplitude Interval 00130 if (!strncmp(strReadLine,"Q_min:Q_step:Q_max = ",21)){ 00131 sscanf(&strReadLine[21],"%lE %lE %lE",&(MyHeader.fAmplitudeMin),&(MyHeader.fAmplitudeStep),&(MyHeader.fAmplitudeMax)); 00132 } 00133 // Read DAC_A Interval High 00134 if (!strncmp(strReadLine,"HI_A_min:HI_A_step:HI_A_max = ",30)){ 00135 sscanf(&strReadLine[30],"%d %d %d",&(MyHeader.fDacMinA[2]),&(MyHeader.fDacStepA[2]),&(MyHeader.fDacMaxA[2])); 00136 } 00137 // Read DAC_A Interval Mid 00138 if (!strncmp(strReadLine,"MID_A_min:MID_A_step:MID_A_max = ",33)){ 00139 sscanf(&strReadLine[33],"%d %d %d",&(MyHeader.fDacMinA[1]),&(MyHeader.fDacStepA[1]),&(MyHeader.fDacMaxA[1])); 00140 } 00141 // Read DAC_A Interval Low 00142 if (!strncmp(strReadLine,"LO_A_min:LO_A_step:LO_A_max = ",30)){ 00143 sscanf(&strReadLine[30],"%d %d %d",&(MyHeader.fDacMinA[0]),&(MyHeader.fDacStepA[0]),&(MyHeader.fDacMaxA[0])); 00144 } 00145 // Read DAC_B Interval High 00146 if (!strncmp(strReadLine,"HI_B_min:HI_B_step:HI_B_max = ",30)){ 00147 sscanf(&strReadLine[30],"%d %d %d",&(MyHeader.fDacMinB[2]),&(MyHeader.fDacStepB[2]),&(MyHeader.fDacMaxB[2])); 00148 } 00149 // Read DAC_B Interval Mid 00150 if (!strncmp(strReadLine,"MID_B_min:MID_B_step:MID_B_max = ",33)){ 00151 sscanf(&strReadLine[33],"%d %d %d",&(MyHeader.fDacMinB[1]),&(MyHeader.fDacStepB[1]),&(MyHeader.fDacMaxB[1])); 00152 } 00153 // Read DAC_B Interval Low 00154 if (!strncmp(strReadLine,"LO_B_min:LO_B_step:LO_B_max = ",30)){ 00155 sscanf(&strReadLine[30],"%d %d %d",&(MyHeader.fDacMinB[0]),&(MyHeader.fDacStepB[0]),&(MyHeader.fDacMaxB[0])); 00156 } 00157 // Read Number of Dump 00158 if (!strncmp(strReadLine,"nb acquisitions = ",18)){ 00159 sscanf(&strReadLine[18],"%d",&(MyHeader.fNbDump)); 00160 } 00161 } 00162 00163 /***************************************************************************************** 00164 Ouverture du fichier root 00165 *****************************************************************************************/ 00166 gSystem->Load("libMicro.so"); 00167 TFile *f = new TFile(rootName.c_str()); 00168 if ( f->IsZombie()) { 00169 FILE_LOG(logERROR) << "Root file ["<< rootName.c_str() << "] does not exist" << endl; 00170 exit(-1); 00171 } 00172 TTree *tree = (TTree *) gROOT->FindObject("t1"); 00173 Int_t iNbEntries=tree->GetEntries(); 00174 FILE_LOG(logINFO) << "Number of reoncstruct events: " << iNbEntries << endl; 00175 MTEvent *evt = new MTEvent(); 00176 TBranch *branch= tree->GetBranch("MTEvent"); 00177 branch->SetAddress(&evt); 00178 00179 /***************************************************************************************** 00180 Vérification du nombre d'évenements 00181 *****************************************************************************************/ 00182 Int_t iNbEvent=0; 00183 for (Double_t IndAmplitude = MyHeader.fAmplitudeMin;IndAmplitude <= (MyHeader.fAmplitudeMax+MyHeader.fAmplitudeStep/2);IndAmplitude += MyHeader.fAmplitudeStep){ 00184 FILE_LOG(logINFO) << "Amplitude : " << IndAmplitude << endl; 00185 00186 for (Int_t IndChannel = MyHeader.fChannelMin;IndChannel <= MyHeader.fChannelMax;IndChannel += MyHeader.fChannelStep){ 00187 FILE_LOG(logINFO) << "\tChannel : " << IndChannel << endl; 00188 for (Int_t IndDAC = 0;IndDAC < 3;IndDAC += 1){ // Low, Mid and High DAC 00189 Int_t fDacMin,fDacStep,fDacMax; 00190 if (IndChannel <32){ fDacMin=MyHeader.fDacMinA[IndDAC];fDacMax=MyHeader.fDacMaxA[IndDAC];fDacStep=MyHeader.fDacStepA[IndDAC]; } 00191 else { fDacMin=MyHeader.fDacMinB[IndDAC];fDacMax=MyHeader.fDacMaxB[IndDAC];fDacStep=MyHeader.fDacStepB[IndDAC]; } 00192 for (Int_t IndDac = fDacMin;IndDac <= fDacMax;IndDac += fDacStep){ 00193 iNbEvent += 8*MyHeader.fNbDump; 00194 } 00195 } 00196 } 00197 } 00198 FILE_LOG(logINFO) << "Nombre d'événements : " << iNbEvent << endl; 00199 00200 if (iNbEvent == iNbEntries){ 00201 00202 /************************************************************************************** 00203 PREPARE HISTOGRAMS 00204 *****************************************************************************************/ 00205 Double_t qMin = MyHeader.fAmplitudeMin,qMax = MyHeader.fAmplitudeMax,qStep = MyHeader.fAmplitudeStep; 00206 Int_t chMin = MyHeader.fChannelMin,chMax = MyHeader.fChannelMax,chStep = MyHeader.fChannelStep; 00207 Int_t dacMinL,dacMaxL,dacStepL,dacMinM,dacMaxM,dacStepM,dacMinH,dacMaxH,dacStepH; 00208 TClonesArray * hListLo = new TClonesArray ("TH2F", 4096) ; // we compute all channel responses ! 00209 TClonesArray * hListMi = new TClonesArray ("TH2F", 4096) ; // we compute all channel responses ! 00210 TClonesArray * hListHi = new TClonesArray ("TH2F", 4096) ; // we compute all channel responses ! 00211 00212 FILE_LOG(logINFO) << "Bins Min : " << qMin << " Max :" << qMax+qStep << " Nb of bins : " << fixed << (qMax+qStep-qMin)/qStep << endl ; 00213 00214 string hBaseName = "h_" ; 00215 string hName ; 00216 string hBaseTitle = "scurve_" ; 00217 string hTitle ; 00218 Int_t index ; 00219 // There are 64 histo for each stimd channel... 00220 // .... and 64(chmax-chmin)/chstep stimd channel 00221 for (Int_t stimdchan=chMin;stimdchan<=chMax;stimdchan+=chStep) 00222 { for (Int_t chan=1;chan<=64;chan++) 00223 { index = 64*(stimdchan)+(chan-1) ; 00224 00225 00226 if (chan ==1) { // Bank A 00227 dacMinL=MyHeader.fDacMinA[0];dacMaxL=MyHeader.fDacMaxA[0];dacStepL=MyHeader.fDacStepA[0]; 00228 dacMinM=MyHeader.fDacMinA[1];dacMaxM=MyHeader.fDacMaxA[1];dacStepM=MyHeader.fDacStepA[1]; 00229 dacMinH=MyHeader.fDacMinA[2];dacMaxH=MyHeader.fDacMaxA[2];dacStepH=MyHeader.fDacStepA[2]; 00230 } 00231 if (chan ==33) { // Bank B 00232 dacMinL=MyHeader.fDacMinB[0];dacMaxL=MyHeader.fDacMaxB[0];dacStepL=MyHeader.fDacStepB[0]; 00233 dacMinM=MyHeader.fDacMinB[1];dacMaxM=MyHeader.fDacMaxB[1];dacStepM=MyHeader.fDacStepB[1]; 00234 dacMinH=MyHeader.fDacMinB[2];dacMaxH=MyHeader.fDacMaxB[2];dacStepH=MyHeader.fDacStepB[2]; 00235 } 00236 00237 hName = hBaseName + "S" + itos(stimdchan+1) + "_H" + itos(chan) + "l"; 00238 hTitle = hBaseTitle + "S" + itos(stimdchan+1) + "_H" + itos(chan) + "l" ; 00239 hListLo->New(index) ; 00240 ( (TH2F *) hListLo->At(index) )->SetNameTitle(hName.c_str(), hTitle.c_str()) ; 00241 ( (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) ; 00242 00243 hName = hBaseName + "S" + itos(stimdchan+1) + "_H" + itos(chan) + "m"; 00244 hTitle = hBaseTitle + "S" + itos(stimdchan+1) + "_H" + itos(chan) + "m" ; 00245 hListMi->New(index) ; 00246 ((TH2F *) hListMi->At(index) )->SetNameTitle(hName.c_str(), hTitle.c_str()) ; 00247 ((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) ; 00248 00249 hName = hBaseName + "S" + itos(stimdchan+1) + "_H" + itos(chan) + "h"; 00250 hTitle = hBaseTitle + "S" + itos(stimdchan+1) + "_H" + itos(chan) + "h" ; 00251 hListHi->New(index) ; 00252 ((TH2F *) hListHi->At(index) )->SetNameTitle(hName.c_str(), hTitle.c_str()) ; 00253 ((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) ; 00254 } 00255 } 00256 FILE_LOG(logINFO) << "Init ok" << endl ; 00257 /***************************************************************************************** 00258 Calcul des Scurves 00259 *****************************************************************************************/ 00260 Int_t currHitValue, currHitChannel; 00261 Int_t iEvent=0; 00262 for (Double_t IndAmplitude = MyHeader.fAmplitudeMin;IndAmplitude <= (MyHeader.fAmplitudeMax+MyHeader.fAmplitudeStep/2);IndAmplitude += MyHeader.fAmplitudeStep){ 00263 00264 for (Int_t IndChannel = MyHeader.fChannelMin;IndChannel <= MyHeader.fChannelMax;IndChannel += MyHeader.fChannelStep){ 00265 for (Int_t IndDAC = 0;IndDAC < 3;IndDAC += 1){ // Low, Mid and High DAC 00266 Int_t fDacMin,fDacStep,fDacMax; 00267 if (IndChannel <32){ fDacMin=MyHeader.fDacMinA[IndDAC];fDacMax=MyHeader.fDacMaxA[IndDAC];fDacStep=MyHeader.fDacStepA[IndDAC]; } 00268 else { fDacMin=MyHeader.fDacMinB[IndDAC];fDacMax=MyHeader.fDacMaxB[IndDAC];fDacStep=MyHeader.fDacStepB[IndDAC]; } 00269 for (Int_t IndDac = fDacMin;IndDac <= fDacMax;IndDac += fDacStep){ 00270 for (Int_t IndEvent = 0;IndEvent < 8*MyHeader.fNbDump;IndEvent ++){ 00271 tree->GetEntry(iEvent); 00272 Int_t NC=evt->GetNchannel(); 00273 for(int i=0;i<NC;i++){ 00274 00275 MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(i); 00276 00277 currHitValue = channel->GetDigitalValue() ; // ...Xtract data from channel... 00278 currHitChannel = channel->GetHardId() ; // ... Xtract channel number... 00279 index = 64*IndChannel+currHitChannel ; 00280 switch(currHitValue){ 00281 case 3: 00282 if (IndDAC == 2){ 00283 ( (TH2F *) hListHi->At(index) )->Fill(IndDac, IndAmplitude) ; 00284 //FILE_LOG(logINFO) << "Amplitude 3 : " << IndChannel << " / " << currHitChannel << endl; 00285 } 00286 break ; 00287 case 2: 00288 if (IndDAC == 1){ 00289 ( (TH2F *) hListMi->At(index) )->Fill(IndDac, IndAmplitude) ; 00290 //FILE_LOG(logINFO) << "Amplitude 2 : " << IndChannel << " / " << currHitChannel << endl; 00291 } 00292 break ; 00293 case 1: 00294 if (IndDAC == 0){ 00295 ( (TH2F *) hListLo->At(index) )->Fill(IndDac, IndAmplitude) ; 00296 //FILE_LOG(logINFO) << "Amplitude 1 : " << IndChannel << " / " << currHitChannel << endl; 00297 } 00298 break ; 00299 default: 00300 cout << "ERROR in ROOT file !" << endl ; 00301 } 00302 } 00303 iEvent++; 00304 } 00305 } 00306 } 00307 } 00308 } 00309 /***************************************************************************************** 00310 NORMALIZE IN %, remove first and last bin 00311 *****************************************************************************************/ 00312 00313 // Normalize to 100% 00314 Float_t scale = 1./(8*MyHeader.fNbDump) ; 00315 for (Int_t stimdchan=chMin;stimdchan<=chMax;stimdchan+=chStep) 00316 { for (Int_t chan=1;chan<=64;chan++) 00317 { index = 64*(stimdchan)+(chan-1) ; 00318 //cout << scale ; 00319 ( (TH2F *) hListLo->At(index) )->Scale(scale) ; 00320 ( (TH2F *) hListMi->At(index) )->Scale(scale) ; 00321 ( (TH2F *) hListHi->At(index) )->Scale(scale) ; 00322 } 00323 } 00324 00325 cout << "Normalization ok" << endl ; 00326 /***************************************************************************************** 00327 Enregistrement des Scurves 00328 *****************************************************************************************/ 00329 //FILE_LOG(logINFO) << "Events : " << iEvent << endl ; 00330 TFile *histFile = new TFile(OutputName.c_str(), "RECREATE") ; 00331 if (histFile->IsOpen() == 0) 00332 { cerr << "Could not open data file "<< endl ; 00333 exit(1) ; 00334 } 00335 histFile->cd() ; 00336 hListLo->Write() ; 00337 hListMi->Write() ; 00338 00339 hListHi->Write() ; 00340 00341 00342 /***************************************************************************************** 00343 Affichage des Scurves 00344 *****************************************************************************************/ 00345 /* argc=0; 00346 TApplication theApp("App", &argc, argv); // see $ROOTSYS/test/hworld.cxx 00347 gStyle->SetPalette(1); 00348 //gStyle->SetCanvasColor(33); 00349 //gStyle->SetFrameFillColor(18); 00350 TCanvas *c1 = new TCanvas("c1","S-Curves") ; 00351 00352 c1->Connect("Closed()", "TApplication", &theApp, "Terminate()"); 00353 string name="h_S32_H32l"; 00354 TH2F* h=(TH2F*) gROOT->FindObject(name.c_str()) ; 00355 c1->cd() ; 00356 h->Draw("lego2") ; 00357 theApp.Run();*/ 00358 00359 hListLo->Delete() ; 00360 hListMi->Delete() ; 00361 hListHi->Delete() ; 00362 evt->Delete() ; 00363 //h->Delete(); 00364 } 00365 return 0; 00366 }
1.4.7