Hi Axel,
Thanks for the answer. Please see bellow the shapeanalysis.cpp file :
#include “TROOT.h”
#include “TClass.h”
#include “TGraph.h”
#include “TF1.h”
#include “TH1F.h”
#include “TH2F.h”
#include “TH3F.h”
#include “TH3D.h”
#include “TTree.h”
#include “TCanvas.h”
#include “TBranch.h”
#include “Riostream.h”
#include “TStyle.h”
#include “TFile.h”
#include “TString.h”
#include “TLegend.h”
#include “TRandom3.h”
#include “TMath.h”
#include “math.h”
#include “TColor.h”
#include “TWaveForm.hh”
#include
#include
#include
#define MAXLENGTH 10000
using namespace std;
typedef vector<Double_t> vect;
class sample{
public:
UInt_t SampleValue[128];
};
class WaveFormStd{
public:
Int_t nb_sample_wf;
Double_t SampleW_F[MAXLENGTH];
Double_t T_0;
};
void shape_analysis();
int main(){
shape_analysis();
return 0;
}
void shape_analysis()
{
//////
////////On instancie la classe des donnees d’enregistrement ////////
///
sample asample1;
sample asample2;
sample asample3;
WaveFormStd astdwaveform;
Double_t Nrj(0.);
Double_t Tim(0.);
Double_t Q(0.);
Double_t t_Origine(0.);
Double_t Chi2(0.);
//////
////////option d’affichage des stats des histogrammes ////////
///
gStyle -> SetStatW(0.28);
gStyle -> SetStatH(0.13);
gStyle -> SetStatColor(41);
gStyle -> SetStatX(0.87);
gStyle -> SetStatY(0.85);
gStyle -> SetStatFont(42);
gStyle->SetOptStat(1111111);
gStyle->SetOptFit(1111);
//////
////// Ouverture du fichier WaveForm.root //////
///
TFile* file=new TFile("WaveForm.root");
TTree * tree_wave_form = (TTree*)file->Get("Wave_Form");
Int_t nentries = (Int_t)tree_wave_form->GetEntries();
cout<<"nentries : "<<nentries<<endl;
TBranch* RF_Faisceau = tree_wave_form->GetBranch("RF");
RF_Faisceau->SetAddress(&asample1);
TBranch* PMT1 = tree_wave_form->GetBranch("PM1");
PMT1->SetAddress(&asample2);
TBranch* PMT2 = tree_wave_form->GetBranch("PM2");
PMT2->SetAddress(&asample3);
//////
////// Ouverture du fichier StandardWaveForm.root //////
///
TFile* file_std_wf=new TFile("StandardWaveForm.root");
TTree * Tree_Energy_Time = (TTree*)file_std_wf->Get("Tree_Energy_Time");
Int_t nentries_std_wf = (Int_t)Tree_Energy_Time->GetEntries();
cout<<"nentries_std_wf : "<<nentries_std_wf<<endl;
TBranch* Energie = Tree_Energy_Time->GetBranch("energy");
Energie->SetAddress(&Nrj);
TBranch* TempsRecal = Tree_Energy_Time->GetBranch("time");
TempsRecal->SetAddress(&Tim);
TBranch* Charge = Tree_Energy_Time->GetBranch("charge");
Charge->SetAddress(&Q);
TTree * Tree_OutPut = (TTree*)file_std_wf->Get("Tree_OutPut");
Int_t nentries_output = (Int_t)Tree_OutPut->GetEntries();
cout<<"nentries_output : "<<nentries_output<<endl;
TBranch* Nb_sample_std = Tree_OutPut->GetBranch("nb_sample_std");
Nb_sample_std->SetAddress(&astdwaveform.nb_sample_wf);
TBranch* Sample_STD = Tree_OutPut->GetBranch("std_waveform");
Sample_STD->SetAddress(astdwaveform.SampleW_F);
TBranch* TempsOrigine = Tree_OutPut->GetBranch("T_Origine");
TempsOrigine->SetAddress(&astdwaveform.T_0);
TFile* Chi2_File=new TFile("Chi2_File.root","recreate");
TTree* Tree_Chi2=new TTree("Tree_Chi2","A Root Tree");
Tree_Chi2->Branch("CHI2",&Chi2,"Chi2/D");
Int_t N_sample=128;
Double_t T_sampling=1.e-9;
Double_t T_samplingUp=T_sampling/5;
Tree_OutPut->GetEntry(0);
t_Origine=astdwaveform.T_0;
cout<<"t_Origine "<<astdwaveform.T_0<<endl;
Int_t Nb_Sample_std=astdwaveform.nb_sample_wf;
cout<<"Nb_Sample_std "<<Nb_Sample_std<<endl;
Int_t Chi2Min=0;
Int_t Chi2Max=10.;
Int_t NBinChi2=2000;
TH1F* Hist_Chi2=new TH1F(“Hist_Chi2”,"",NBinChi2,Chi2Min,Chi2Max);
TWaveForm* W_F_STD=new TWaveForm(Nb_Sample_std,t_Origine,T_samplingUp,astdwaveform.SampleW_F,"wave_form_ref");
TGraph* Graph_STD=W_F_STD->Graph();
TCanvas* c1 = new TCanvas(“c1”,“Control”,50,10,600,600);
c1->cd();
c1->Print(“c1.svg”);
Graph_STD->Draw(“AL*”);
Graph_STD->SetMarkerColor(2);
Graph_STD->SetLineColor(2);
Double_t* pm1;
TWaveForm** wf0=new TWaveForm*[nentries];
Double_t PiedEstal=172.;
TWaveForm* wf_trifouillage;
TWaveForm* wf_oversampled;
for (Int_t i = 0 ; i < nentries ;i++){
pm1=new Double_t[N_sample];
tree_wave_form->GetEntry(i);
Tree_Energy_Time->GetEntry(i);
for (Int_t j = 0 ; j < N_sample ; j++){
pm1[j]=asample2.SampleValue[j];
}
//Int_t Origine_des_temps=Tim/T_sampling+0.5; /// Arrondi et non troncature !
//wf0[i]=new TWaveForm(N_sample,-Origine_des_temps*T_sampling,T_sampling,pm1,"wave_form");
wf0[i]=new TWaveForm(N_sample,0.,T_sampling,pm1,"wave_form");
(*wf0[i])-=PiedEstal; /// soustraction du piedestal
(*wf0[i])/=Q; /// normalisation en charge
wf_oversampled=wf0[i]->UpSample(5);
Int_t Origine_des_temps=Tim/T_samplingUp+0.5; /// Arrondi et non troncature !
wf_oversampled->Set_t0(-Origine_des_temps*T_samplingUp);
wf_trifouillage=new TWaveForm((*wf_oversampled));
/*if (i%1000==0){
wf_trifouillage->Draw("LPsame");
c1->Update();
}*/
(*wf_trifouillage)-=(*W_F_STD); /// Difference point par point
/// entre le i eme evenement et la waveform de ref std
Chi2=wf_trifouillage->chi_square(); /// Retourne la moyenne de la diff quad, le Chi2 quoi !
Tree_Chi2->Fill();
Hist_Chi2->Fill(Chi2);
}
Int_t EMin=0;
Int_t EMax=2000;
Int_t NBinE=2000;
TH1F* EnergieDeposee=new TH1F("EnergieDeposee","",NBinE,EMin,EMax);
TH1F* EnergieDeposeeCut=new TH1F("EnergieDeposeeCut","",NBinE,EMin,EMax);
TH2F* Chi2_vs_Energy=new TH2F("Chi2_vs_Energy","",NBinChi2,Chi2Min,Chi2Max,NBinE,EMin,EMax);
TCanvas* c2 = new TCanvas("c2","Chi2",50,10,600,600);
c2->cd();
Hist_Chi2->Draw();
c2->Print("c2.svg");
for (Int_t i = 0 ; i < nentries; i++){
Tree_Chi2->GetEntry(i);
Tree_Energy_Time->GetEntry(i);
if (Nrj>1600 && Nrj<1700){
EnergieDeposeeCut->Fill(Nrj);
Chi2_vs_Energy->Fill(Chi2/1.e17,Nrj);
}
EnergieDeposee->Fill(Nrj);
//Chi2_vs_Energy->Fill(Chi2/1.e17,Nrj);
}
TCanvas* c3 = new TCanvas("c3","Nrj",50,10,600,600);
EnergieDeposee->Draw();
EnergieDeposee->SetLineColor(2);
EnergieDeposeeCut->Draw("same");
EnergieDeposeeCut->SetLineColor(4);
c3->Print("c3.svg");
TCanvas* c4 = new TCanvas("c4","Chi2_Nrj",50,10,600,600);
//c4->SetLogz(1);
gStyle->SetPalette(1);
Chi2_vs_Energy->Draw("colz");
Double_t E1=EnergieDeposee->Integral();
Double_t E2=EnergieDeposeeCut->Integral();
cout<<"E1 "<<E1<<" E2 "<<E2<<endl;
c4->Print("c4.svg");
delete EnergieDeposee;
delete EnergieDeposeeCut;
delete Chi2_vs_Energy;
Chi2_File->Write();
Chi2_File->Close();
}
Thanks a lot
Cheers,
Loïc