/[PAMELA software]/quicklook/tracker/flight/macros/fillpedsig.cpp
ViewVC logotype

Contents of /quicklook/tracker/flight/macros/fillpedsig.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (show annotations) (download) (vendor branch)
Tue Mar 14 16:09:30 2006 UTC (18 years, 8 months ago) by pam-fi
Branch: trk-QLook
CVS Tags: R1v00, R1v01
Changes since 1.1: +0 -0 lines
First CVS release of tracker flight quick-look 

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 }

  ViewVC Help
Powered by ViewVC 1.1.23