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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (hide annotations) (download)
Wed May 10 10:05:56 2006 UTC (18 years, 6 months ago) by pam-fi
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +0 -0 lines
FILE REMOVED
.cpp files of utilities removed, quicklook macros added

1 pam-fi 1.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