1 |
#include <utils/yodaUtility.h> |
2 |
#include "TrkFunctions.cpp" |
3 |
#include <iostream> |
4 |
#include <fstream> |
5 |
#include <TPaveText.h> |
6 |
#include <TLatex.h> |
7 |
#include <TCanvas.h> |
8 |
#include <TGraph.h> |
9 |
#include <TStyle.h> |
10 |
#include <TObjString.h> |
11 |
#include <PscuHeader.h> |
12 |
#include <EventHeader.h> |
13 |
#include <CalibTrk1Event.h> |
14 |
#include <CalibTrk2Event.h> |
15 |
|
16 |
/** |
17 |
* fillpedsig --> is a function to calculate the mean value for |
18 |
* sigma, variance of sigma, pedestal and variance of pedestal |
19 |
* and fill a root.file whit these values |
20 |
* |
21 |
* autor: D.Fedele |
22 |
* version 1.0 |
23 |
* Parameters: |
24 |
* file - the path to the root file to analyse |
25 |
* outdir - total path of output file |
26 |
* |
27 |
*/ |
28 |
|
29 |
void fillpedsig(TString file,TString outdir) |
30 |
{ |
31 |
// |
32 |
// obtain information about the data file and select the output file |
33 |
const string filepath=file.Data(); |
34 |
Int_t dwpos = filepath.find("DW_"); |
35 |
Int_t dwpos1 = filepath.find(".root"); |
36 |
TString fpath=(filepath.c_str()); |
37 |
TString base,ffile,date1,date2,date; |
38 |
stringcopy(base,fpath,0,dwpos); |
39 |
stringcopy(ffile,fpath,dwpos,dwpos1); |
40 |
|
41 |
stringcopy(date1,fpath,dwpos+3,dwpos+9); |
42 |
stringcopy(date2,fpath,dwpos+10,dwpos1); |
43 |
|
44 |
TString out; |
45 |
stringstream outfile; |
46 |
if(outdir.Length()==0){ |
47 |
out = base; |
48 |
}else{ |
49 |
out = outdir; |
50 |
}; |
51 |
outfile<<out.Data()<<"pedsig.root"; |
52 |
|
53 |
// |
54 |
// initialise structure where data will store |
55 |
struct caltrk_def ctrk; |
56 |
Int_t nevents=0; |
57 |
Float_t OBT[2]; |
58 |
|
59 |
OBT[0]=0; |
60 |
OBT[1]=0; |
61 |
ctrk.good0[0]=0; |
62 |
ctrk.good0[1]=0; |
63 |
for(Int_t i=0;i<12;i++){ |
64 |
ctrk.daqmode[i]=0; |
65 |
ctrk.dspnum[i]=0; |
66 |
ctrk.calibnum[i]=0; |
67 |
ctrk.ncalev[i]=0; |
68 |
ctrk.calfl[i]=0; |
69 |
ctrk.ped1[i]=0; |
70 |
ctrk.ped2[i]=0; |
71 |
ctrk.ped3[i]=0; |
72 |
ctrk.sig1[i]=0; |
73 |
ctrk.sig2[i]=0; |
74 |
ctrk.sig3[i]=0; |
75 |
ctrk.nbad1[i]=0; |
76 |
ctrk.nbad2[i]=0; |
77 |
ctrk.nbad3[i]=0; |
78 |
ctrk.crc_hc[i]=0; |
79 |
ctrk.crc_c[i][0]=0; |
80 |
ctrk.crc_c[i][1]=0; |
81 |
ctrk.crc_c[i][2]=0; |
82 |
for(Int_t iii=0;iii<3072;iii++){ |
83 |
ctrk.dspped[i][iii]=0; |
84 |
ctrk.dspsig[i][iii]=0; |
85 |
ctrk.dspbad[i][iii]=0; |
86 |
} |
87 |
} |
88 |
|
89 |
// |
90 |
// open data files |
91 |
pamela::EventHeader *eh1=0,*eh2=0; |
92 |
pamela::PscuHeader *ph1=0,*ph2=0; |
93 |
pamela::CalibTrk1Event *trk1 = 0; |
94 |
pamela::CalibTrk2Event *trk2 = 0; |
95 |
|
96 |
TFile *datafile = new TFile(file); |
97 |
if ( !datafile ){ |
98 |
// printf("No data file, exiting...\n"); |
99 |
return; |
100 |
}; |
101 |
TTree *otr1,*otr2; |
102 |
|
103 |
otr1 = (TTree*)datafile->Get("CalibTrk1"); |
104 |
otr1->SetBranchAddress("CalibTrk1", &trk1); |
105 |
otr1->SetBranchAddress("Header",&eh1); |
106 |
otr2 = (TTree*)datafile->Get("CalibTrk2"); |
107 |
otr2->SetBranchAddress("CalibTrk2", &trk2); |
108 |
otr2->SetBranchAddress("Header",&eh2); |
109 |
|
110 |
|
111 |
if(otr1->GetEntries()==otr2->GetEntries()) |
112 |
nevents = otr1->GetEntries(); |
113 |
else{ |
114 |
// printf("WARNING: CalibTrk1 entries is different from CalibTrk2 entries"); |
115 |
return;} |
116 |
|
117 |
if (nevents<=0) { |
118 |
datafile->Close(); |
119 |
// printf("No calibration packets found, exiting...\n"); |
120 |
return; |
121 |
}; |
122 |
|
123 |
// |
124 |
// Define variables |
125 |
Int_t ndsp =0,numbad[12][12]; |
126 |
Float_t pedav[12][12],pedavtemp[12][12],sigav[12][12],sigavtemp[12][12]; |
127 |
Float_t pedvar[12][12],pedvartemp[12][12],sigvar[12][12],sigvartemp[12][12]; |
128 |
Long64_t datefile=0; |
129 |
|
130 |
|
131 |
// |
132 |
// open output files |
133 |
TFile *pedsig; |
134 |
TTree *pstree; |
135 |
ifstream pedsigtmp(outfile.str().c_str()); |
136 |
if(pedsigtmp.fail()){ |
137 |
pedsigtmp.close(); |
138 |
pedsig = new TFile(outfile.str().c_str(),"new"); |
139 |
pstree=new TTree("pstree","pedsig"); |
140 |
pstree->Branch("pedav",pedav,"pedav[12][12]/F"); |
141 |
pstree->Branch("pedvar",pedvar,"pedvar[12][12]/F"); |
142 |
pstree->Branch("sigav",sigav,"sigav[12][12]/F"); |
143 |
pstree->Branch("sigvar",sigvar,"sigvar[12][12]/F"); |
144 |
pstree->Branch("datefile",&datefile,"datefile/L"); |
145 |
pstree->Branch("numbad",numbad,"numbad[12][12]/I"); |
146 |
} |
147 |
if(!pedsigtmp.fail()){ |
148 |
pedsigtmp.close(); |
149 |
pedsig = new TFile(outfile.str().c_str(),"update"); |
150 |
pstree = (TTree*)pedsig->Get("pstree"); |
151 |
pstree->SetBranchAddress("pedav",pedav); |
152 |
pstree->SetBranchAddress("pedvar",pedvar); |
153 |
pstree->SetBranchAddress("sigav",sigav); |
154 |
pstree->SetBranchAddress("sigvar",sigvar); |
155 |
pstree->SetBranchAddress("datefile",&datefile); |
156 |
pstree->SetBranchAddress("numbad",numbad); |
157 |
} |
158 |
|
159 |
datefile=atoll(date1+date2); |
160 |
//********************************************************************** |
161 |
// |
162 |
// LOOP OVER EVENTS |
163 |
// |
164 |
//********************************************************************** |
165 |
|
166 |
for (Int_t i = 0; i < nevents; i++){ |
167 |
|
168 |
otr1->GetEntry(i); |
169 |
otr2->GetEntry(i); |
170 |
|
171 |
// |
172 |
// fill data storage structure |
173 |
ctrk.good0[0]=trk1->good0; |
174 |
ctrk.good0[1]=trk2->good0; |
175 |
for (Int_t m = 0; m < 6; m++){ |
176 |
Int_t plane=m; |
177 |
ph1 = eh1->GetPscuHeader(); |
178 |
OBT[0]= ph1->GetOrbitalTime(); |
179 |
ctrk.daqmode[trk1->DSPnumber[m]-1]=trk1->DAQmode[m]; |
180 |
ctrk.dspnum[trk1->DSPnumber[m]-1]=trk1->DSPnumber[m]; |
181 |
ctrk.calibnum[trk1->DSPnumber[m]-1]=trk1->calibnumber[m]; |
182 |
ctrk.ncalev[trk1->DSPnumber[m]-1]=trk1->ncalib_event[m]; |
183 |
ctrk.ped1[trk1->DSPnumber[m]-1]=trk1->ped_l1[m]; |
184 |
ctrk.ped2[trk1->DSPnumber[m]-1]=trk1->ped_l2[m]; |
185 |
ctrk.ped3[trk1->DSPnumber[m]-1]=trk1->ped_l3[m]; |
186 |
ctrk.sig1[trk1->DSPnumber[m]-1]=trk1->sig_l1[m]; |
187 |
ctrk.sig2[trk1->DSPnumber[m]-1]=trk1->sig_l2[m]; |
188 |
ctrk.sig3[trk1->DSPnumber[m]-1]=trk1->sig_l3[m]; |
189 |
ctrk.nbad1[trk1->DSPnumber[m]-1]=trk1->nbad_l1[m]; |
190 |
ctrk.nbad2[trk1->DSPnumber[m]-1]=trk1->nbad_l2[m]; |
191 |
ctrk.nbad3[trk1->DSPnumber[m]-1]=trk1->nbad_l3[m]; |
192 |
ctrk.calfl[trk1->DSPnumber[m]-1]=trk1->cal_flag[m]; |
193 |
ctrk.crc_c[trk1->DSPnumber[m]-1][0]=trk1->crc_cal[m][0]; |
194 |
ctrk.crc_c[trk1->DSPnumber[m]-1][1]=trk1->crc_cal[m][1]; |
195 |
ctrk.crc_c[trk1->DSPnumber[m]-1][2]=trk1->crc_cal[m][2]; |
196 |
ctrk.crc_hc[trk1->DSPnumber[m]-1]=trk1->crc_hcal[m]; |
197 |
for (Int_t j = 0; j < 3072; j++){ |
198 |
ctrk.dspped[trk1->DSPnumber[m]-1][j]=trk1->DSPped_par[m][j]; |
199 |
ctrk.dspsig[trk1->DSPnumber[m]-1][j]=trk1->DSPsig_par[m][j]; |
200 |
ctrk.dspbad[trk1->DSPnumber[m]-1][j]=trk1->DSPbad_par[m][j]; |
201 |
} |
202 |
ph2 = eh2->GetPscuHeader(); |
203 |
OBT[1]= ph2->GetOrbitalTime(); |
204 |
ctrk.daqmode[trk2->DSPnumber[m]-1]=trk2->DAQmode[m]; |
205 |
ctrk.dspnum[trk2->DSPnumber[m]-1]=trk2->DSPnumber[m]; |
206 |
ctrk.calibnum[trk2->DSPnumber[m]-1]=trk2->calibnumber[m]; |
207 |
ctrk.ncalev[trk2->DSPnumber[m]-1]=trk2->ncalib_event[m]; |
208 |
ctrk.ped1[trk2->DSPnumber[m]-1]=trk2->ped_l1[m]; |
209 |
ctrk.ped2[trk2->DSPnumber[m]-1]=trk2->ped_l2[m]; |
210 |
ctrk.ped3[trk2->DSPnumber[m]-1]=trk2->ped_l3[m]; |
211 |
ctrk.sig1[trk2->DSPnumber[m]-1]=trk2->sig_l1[m]; |
212 |
ctrk.sig2[trk2->DSPnumber[m]-1]=trk2->sig_l2[m]; |
213 |
ctrk.sig3[trk2->DSPnumber[m]-1]=trk2->sig_l3[m]; |
214 |
ctrk.nbad1[trk2->DSPnumber[m]-1]=trk2->nbad_l1[m]; |
215 |
ctrk.nbad2[trk2->DSPnumber[m]-1]=trk2->nbad_l2[m]; |
216 |
ctrk.nbad3[trk2->DSPnumber[m]-1]=trk2->nbad_l3[m]; |
217 |
ctrk.calfl[trk2->DSPnumber[m]-1]=trk2->cal_flag[m]; |
218 |
ctrk.crc_c[trk1->DSPnumber[m]-1][0]=trk2->crc_cal[m][0]; |
219 |
ctrk.crc_c[trk1->DSPnumber[m]-1][1]=trk2->crc_cal[m][1]; |
220 |
ctrk.crc_c[trk1->DSPnumber[m]-1][2]=trk2->crc_cal[m][2]; |
221 |
ctrk.crc_hc[trk1->DSPnumber[m]-1]=trk2->crc_hcal[m]; |
222 |
for (Int_t j = 0; j < 3072; j++){ |
223 |
ctrk.dspped[trk2->DSPnumber[m]-1][j]=trk2->DSPped_par[m][j]; |
224 |
ctrk.dspsig[trk2->DSPnumber[m]-1][j]=trk2->DSPsig_par[m][j]; |
225 |
ctrk.dspbad[trk2->DSPnumber[m]-1][j]=trk2->DSPbad_par[m][j]; |
226 |
} |
227 |
} |
228 |
|
229 |
// |
230 |
// inizialise variables |
231 |
for(Int_t n = 0; n<12; n++){ |
232 |
for(Int_t nm = 0; nm<12; nm++){ |
233 |
pedav[n][nm]=0; |
234 |
pedavtemp[n][nm]=0; |
235 |
sigav[n][nm]=0; |
236 |
sigavtemp[n][nm]=0; |
237 |
pedvar[n][nm]=0; |
238 |
pedvartemp[n][nm]=0; |
239 |
sigvar[n][nm]=0; |
240 |
sigvartemp[n][nm]=0; |
241 |
numbad[n][nm]=0; |
242 |
} |
243 |
} |
244 |
|
245 |
// |
246 |
// calculate the mean value for sigma, variance of sigma, pedestal, variance of pedestal |
247 |
Int_t nn=0; |
248 |
for(Int_t n = 0; n<12; n++){ |
249 |
|
250 |
ndsp = ctrk.dspnum[n]; |
251 |
nn = ndsp-1; |
252 |
|
253 |
for(Int_t j = 0; j < 3072; j++){ |
254 |
if(ctrk.dspbad[nn][j]==0){ |
255 |
sigavtemp[nn][j/256]+=ctrk.dspsig[nn][j]; |
256 |
pedavtemp[nn][j/256]+=ctrk.dspped[nn][j]; |
257 |
} |
258 |
else |
259 |
numbad[nn][j/256]+=1; |
260 |
} |
261 |
for(Int_t ii=0;ii<12;ii++){ |
262 |
pedav[nn][ii]=pedavtemp[nn][ii]/256; |
263 |
sigav[nn][ii]=sigavtemp[nn][ii]/256; |
264 |
} |
265 |
for(Int_t j = 0; j < 3072; j++){ |
266 |
if(ctrk.dspbad[nn][j]==0){ |
267 |
sigvartemp[nn][j/256]+=pow((ctrk.dspsig[nn][j]-sigav[nn][j/256]),2); |
268 |
pedvartemp[nn][j/256]+=pow((ctrk.dspped[nn][j]-pedav[nn][j/256]),2); |
269 |
}; |
270 |
} |
271 |
for(Int_t ii=0;ii<12;ii++){ |
272 |
pedvar[nn][ii]=sqrt(pedvartemp[nn][ii]/255); |
273 |
sigvar[nn][ii]=sqrt(sigvartemp[nn][ii]/255); |
274 |
} |
275 |
} |
276 |
|
277 |
// |
278 |
// Fill the output Tree |
279 |
pstree->Fill(); |
280 |
|
281 |
};//end loop on events |
282 |
|
283 |
// |
284 |
//close all files |
285 |
datafile->Close(); |
286 |
pedsig->Write(); |
287 |
pedsig->Close(); |
288 |
|
289 |
return; |
290 |
} |