source: trunk/MagicSoft/Mars/mtemp/mifae/programs/calib.cc@ 6139

Last change on this file since 6139 was 5069, checked in by blanch, 21 years ago
*** empty log message ***
File size: 9.2 KB
Line 
1/*********************************/
2/* Compute the hillas parameters */
3/*********************************/
4
5#include "TString.h"
6#include "TArrayS.h"
7
8#include "MParList.h"
9#include "MTaskList.h"
10#include "MPedestalCam.h"
11#include "MBadPixelsCam.h"
12#include "MBadPixelsTreat.h"
13#include "MReadMarsFile.h"
14#include "MGeomApply.h"
15#include "MPedCalcPedRun.h"
16#include "MEvtLoop.h"
17#include "MGeomCamMagic.h"
18#include "MExtractedSignalCam.h"
19#include "MCalibrationChargeCam.h"
20#include "MCalibrationQECam.h"
21#include "MCalibrationQEPix.h"
22#include "MHCalibrationChargeCam.h"
23#include "MHCalibrationRelTimeCam.h"
24#include "MExtractor.h"
25#include "MExtractFixedWindow.h"
26#include "MExtractSlidingWindow.h"
27#include "MExtractSignal.h"
28#include "MCalibrationChargeCalc.h"
29#include "MFCosmics.h"
30#include "MContinue.h"
31#include "MLog.h"
32#include "MCerPhotEvt.h"
33#include "MPedPhotCam.h"
34#include "MCalibrate.h"
35#include "MPedPhotCalc.h"
36#include "MHillas.h"
37#include "MNewImagePar.h"
38#include "MRawRunHeader.h"
39#include "MSrcPosCam.h"
40#include "MBlindPixelCalc.h"
41#include "MImgCleanStd.h"
42#include "MHillasSrcCalc.h"
43#include "MHillasCalc.h"
44#include "MArrivalTimeCam.h"
45#include "MArrivalTimeCalc2.h"
46#include "MIslands.h"
47#include "MIslandCalc.h"
48#include "MIslandClean.h"
49#include "MWriteRootFile.h"
50#include "MArgs.h"
51#include "MRunIter.h"
52#include "MJPedestal.h"
53#include "MJCalibration.h"
54#include "MHillasDisplay.h"
55#include "MF.h"
56#include "MContinue.h"
57
58#include "TApplication.h"
59
60#include <iostream>
61#include <fstream>
62#include <stdlib.h>
63
64using namespace std;
65
66Bool_t readDatacards(TString& filename);
67void calib();
68
69// initial and final time slices to be used in signal extraction
70const Byte_t hifirst = 1;
71const Byte_t hilast = 14;
72const Byte_t lofirst = 3;
73const Byte_t lolast = 14;
74
75// declaration of variables read from datacards
76TString outname;
77TString idirname;
78TString filter;
79TString psfilename("makehillas.ps");
80MRunIter caliter;
81MRunIter pediter;
82MRunIter datiter;
83UInt_t display = 0;
84ULong_t nmaxevents= 999999999;
85Short_t calflag = 1;
86Float_t lcore = 3.0;
87Float_t ltail = 1.5;
88Int_t islflag = 0;
89Float_t lnew = 40;
90Int_t kmethod = 1;
91Int_t kalgorithm = 1;
92Int_t nfiles = 0;
93
94const TString defaultcard="makehillas.datacard";
95/*************************************************************/
96static void Usage()
97{
98 gLog <<endl;
99 gLog << "Usage is:" << endl;
100 gLog << " calib [-h] [-?] <datacards>" << endl << endl;
101 gLog << " <datacards>: datacards file name (dafault input.datacards)" << endl;
102 gLog << " -?/-h: This help" << endl << endl;
103}
104/*************************************************************/
105int main(int argc, char **argv)
106{
107 // create a TApplication to be able to
108 TApplication app("Application",0,0);
109
110 // evaluate arguments
111 MArgs arg(argc, argv);
112 if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help"))
113 {
114 Usage();
115 return -1;
116 }
117
118 TString datacard = arg.GetArgumentStr(0);
119 if(!datacard.Length())
120 datacard = defaultcard;
121
122 if(!readDatacards(datacard))
123 {
124 cout << "Error reading datacards. Stoping" << endl;
125 return -1;
126 }
127
128 calib();
129}
130
131/*************************************************************/
132void calib()
133{
134 // Set the general tasks/containers
135 MExtractFixedWindow extractor;
136 extractor.SetRange(hifirst,hilast,lofirst,lolast);
137
138 // MExtractSlidingWindow extractor;
139 // extractor.SetRange(hifirst,hilast,lofirst,lolast);
140 // extractor.SetWindowSize(2,2);
141
142 MGeomCamMagic geomcam;
143 MGeomApply geomapl;
144
145 /************************************/
146 /* FIRST LOOP: PEDESTAL COMPUTATION */
147 /************************************/
148
149 // If you want to exclude pixels from the beginning, read
150 // an ascii-file with the corr. pixel numbers (see MBadPixelsCam)
151 //badcam.AsciiRead("badpixels.dat");
152
153 MJPedestal pedloop;
154 pedloop.SetInput(&pediter);
155 pedloop.SetExtractor(&extractor);
156 // pedloop.SetBadPixels(badcam);
157
158 if (!pedloop.Process())
159 return;
160
161 /*****************************/
162 /* SECOND LOOP: CALIBRATION */
163 /*****************************/
164
165 MJCalibration calloop;
166
167 calloop.SetExtractor(&extractor);
168 calloop.SetInput(&caliter);
169 calloop.SetBadPixels(pedloop.GetBadPixels());
170 if(calflag==2)
171 calloop.SetUseBlindPixel();
172
173 if(calflag>0)
174 if (!calloop.Process(pedloop.GetPedestalCam()))
175 return;
176
177 Float_t meanCharge[577];
178 Float_t meanSigma[577];
179 Float_t meanFADC2Phe[577];
180 Float_t meanFADC2Phe_error[577];
181 Float_t meanFADCtoPh[577];
182 Float_t prob[577];
183 ofstream fout;
184
185 fout.open(outname);
186
187 for (Int_t i=0; i<577; i++)
188 {
189 meanCharge[i] = 0;
190 MCalibrationChargePix calpix = (MCalibrationChargePix&)(calloop.GetCalibrationCam())[i];
191 MCalibrationQEPix qepix = ( MCalibrationQEPix&)(calloop.GetQECam())[i];
192 meanCharge[i] = calpix.GetMean();
193 meanSigma[i] = calpix.GetSigma();
194 meanFADC2Phe[i]=calpix.GetMeanConvFADC2Phe();
195 meanFADC2Phe_error[i]=calpix.GetMeanConvFADC2PheErr();
196 meanFADCtoPh[i]=calpix.GetMeanConvFADC2Phe()/qepix.GetQECascadesFFactor(0);
197 prob[i]=calpix.GetProb();
198 fout << i << '\t' << meanCharge[i] << '\t' << meanSigma[i]<< '\t' << meanFADC2Phe[i]<< '\t' << meanFADC2Phe_error[i]<< '\t' <<meanFADCtoPh[i]<< '\t' <<prob[i]<< '\t' <<(Int_t)calpix.IsHiGainSaturation()<< '\t' <<(Int_t)calpix.IsFFactorMethodValid()<<endl;
199 }
200
201 fout.close();
202
203 /*****************************************************/
204 /* THIRD LOOP: Q in fadc event by event CALIBRATION */
205 /*****************************************************/
206
207 MParList plist4;
208 MTaskList tlist4;
209 plist4.AddToList(&tlist4);
210
211 // containers
212
213 MCerPhotEvt nphot;
214 plist4.AddToList(&geomcam);
215 plist4.AddToList(&pedloop.GetPedestalCam());
216 plist4.AddToList(&nphot);
217
218 //tasks
219 MReadMarsFile read4("Events");
220 static_cast<MRead&>(read4).AddFiles(caliter);
221 read4.DisableAutoScheme();
222
223 MCalibrate::CalibrationMode_t calMode=MCalibrate::kNone;
224
225 MCalibrate photcalc(calMode);
226
227 tlist4.AddToList(&read4);
228 tlist4.AddToList(&geomapl);
229 tlist4.AddToList(&extractor);
230 tlist4.AddToList(&photcalc);
231
232 // Create and setup the eventloop
233 MEvtLoop datloop;
234 datloop.SetParList(&plist4);
235
236 // if (!datloop.Eventloop(nmaxevents))
237 // return;
238 outname+="_evt";
239
240 datloop.PreProcess();
241
242 fout.open(outname);
243 while(tlist4.Process())
244 {
245 for(int i=0;i<577;i++)
246 fout<<nphot[i].GetNumPhotons()<<" ";
247 fout<<endl;
248 }
249 fout.close();
250 datloop.PostProcess();
251 return;
252
253}
254//-------------------------------------------------------------------------------
255
256Bool_t readDatacards(TString& filename)
257{
258 ifstream ifun(filename.Data());
259 if(!ifun)
260 {
261 cout << "File " << filename << " not found" << endl;
262 return kFALSE;
263 }
264
265 TString word;
266
267 while(ifun >> word)
268 {
269 // skip comments
270 if(word[0]=='/' && word[1]=='/')
271 {
272 while(ifun.get()!='\n'); // skip line
273 continue;
274 }
275
276 // input file directory
277 if(strcmp(word.Data(),"IDIR")==0)
278 {
279 if(idirname.Length())
280 cout << "readDataCards Warning: overriding input directory file name" << endl;
281 ifun >> idirname;
282 }
283
284 // pedestal runs
285 if(strcmp(word.Data(),"PRUNS")==0)
286 {
287 if(pediter.GetNumRuns())
288 cout << "readDataCards Warning: adding pedestal runs to the existing list" << endl;
289 ifun >> word;
290 pediter.AddRuns(word.Data(),idirname.Data());
291 }
292
293 // calibration runs
294 if(strcmp(word.Data(),"CRUNS")==0)
295 {
296 if(caliter.GetNumRuns())
297 cout << "readDataCards Warning: adding calibration runs to the existing list" << endl;
298 ifun >> word;
299 caliter.AddRuns(word.Data(),idirname.Data());
300 }
301
302 // output file name
303 if(strcmp(word.Data(),"OUTFILE")==0)
304 {
305 if(outname.Length())
306 cout << "readDataCards Warning: overriding output file name" << endl;
307 ifun >> outname;
308 }
309
310 // calibration flag
311 if(strcmp(word.Data(),"CALFLAG")==0)
312 ifun >> calflag;
313
314
315 }
316
317 pediter.Reset();
318 caliter.Reset();
319 datiter.Reset();
320 TString pfile;
321
322 // Dump read values
323 cout << "************************************************" << endl;
324 cout << "* Datacards read from file " << filename << endl;
325 cout << "************************************************" << endl;
326 cout << "Pedestal file (s): " << endl;
327 while(!(pfile=pediter.Next()).IsNull())
328 cout << pfile << endl;
329 cout << "Calibration file (s): " << endl;
330 while(!(pfile=caliter.Next()).IsNull())
331 cout << pfile << endl;
332 cout << "Output file name: " << outname << endl;
333 cout << "Calibration: ";
334 if(calflag==0)
335 cout << "Pixel area proportional intercalibration" << endl;
336 else if(calflag==-1)
337 cout << "No calibration whatsoever" << endl;
338 else if(calflag==1)
339 cout << "Default calibration" << endl;
340 else if(calflag==11)
341 cout << "Default calibration + bad pixels interpolation" << endl;
342
343 if(!pediter.GetNumEntries())
344 {
345 cout << "No pedestal file name specified" << endl;
346 return kFALSE;
347 }
348 if(!caliter.GetNumEntries() && calflag>0)
349 {
350 cout << "No calibration file name specified" << endl;
351 return kFALSE;
352 }
353 if(!outname.Length())
354 {
355 cout << "No output file name specified" << endl;
356 return kFALSE;
357 }
358
359
360 return kTRUE;
361}
362
Note: See TracBrowser for help on using the repository browser.