/*********************************/ /* Compute the hillas parameters */ /*********************************/ #include "TString.h" #include "TArrayS.h" #include "MParList.h" #include "MTaskList.h" #include "MPedestalCam.h" #include "MBadPixelsCam.h" #include "MReadMarsFile.h" #include "MGeomApply.h" #include "MPedCalcPedRun.h" #include "MEvtLoop.h" #include "MGeomCamMagic.h" #include "MExtractedSignalCam.h" #include "MCalibrationChargeCam.h" #include "MHCalibrationChargeCam.h" #include "MHCalibrationRelTimeCam.h" #include "MExtractFixedWindow.h" #include "MCalibrationChargeCalc.h" #include "MFCosmics.h" #include "MContinue.h" #include "MFillH.h" #include "MLog.h" #include "MCerPhotEvt.h" #include "MPedPhotCam.h" #include "MCalibrate.h" #include "MPedPhotCalc.h" #include "MHillas.h" #include "MRawRunHeader.h" #include "MSrcPosCam.h" #include "MBlindPixelCalc.h" #include "MImgCleanStd.h" #include "MHillasSrcCalc.h" #include "MHillasCalc.h" #include "MWriteRootFile.h" #include "MProgressBar.h" #include "MArgs.h" #include "MRunIter.h" #include "MJPedestal.h" #include "MJCalibration.h" #include #include #include using namespace std; Bool_t readDatacards(TString& filename); void makeHillas(); // initial and final time slices to be used in signal extraction const Byte_t hifirst = 1; const Byte_t hilast = 14; const Byte_t lofirst = 4; const Byte_t lolast = 13; // declaration of variables read from datacards TString outname; TString idirname; MRunIter caliter; MRunIter pediter; MRunIter datiter; ULong_t nmaxevents=999999999; Short_t calflag=1; Float_t lcore = 3.0; Float_t ltail = 1.5; Int_t nfiles = 0; const TString defaultcard="input.datacard"; /*************************************************************/ // makeHillas usage output static void Usage() { gLog <" << endl << endl; gLog << " : datacards file name (dafault input.datacards)" << endl; gLog << " -?/-h: This help" << endl << endl; } /*************************************************************/ // main program int main(int argc, char **argv) { // evaluate arguments MArgs arg(argc, argv); if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help")) { Usage(); return -1; } // get name of input datacard file TString datacard = arg.GetArgumentStr(0); if(!datacard.Length()) datacard = defaultcard; // read the datacards if(!readDatacards(datacard)) { cout << "Error reading datacards. Stoping" << endl; return -1; } // make the hillas file makeHillas(); } /*************************************************************/ /*************************************************************/ /*************************************************************/ void makeHillas() { // set the signal extractor and calibration mode MExtractFixedWindow extractor; extractor.SetRange(hifirst,hilast,lofirst,lolast); MCalibrate::CalibrationMode_t calMode=MCalibrate::kDefault; if(calflag==0) calMode=MCalibrate::kNone; if(calflag==-1) calMode=MCalibrate::kDummy; MCalibrate calibrate(calMode); // general containers MBadPixelsCam badcam; MGeomCamMagic geomcam; MGeomApply geomapl; // If you want to exclude pixels from the beginning, read // an ascii-file with the corr. pixel numbers (see MBadPixelsCam) // // badcam.AsciiRead("badpixels.dat"); /************************************/ /* FIRST LOOP: PEDESTAL COMPUTATION */ /************************************/ MJPedestal pedloop; pedloop.SetInput(&pediter); pedloop.SetBadPixels(badcam); cout << "*************************" << endl; cout << "** COMPUTING PEDESTALS **" << endl; cout << "*************************" << endl; if (!pedloop.Process()) return; /*****************************/ /* SECOND LOOP: CALIBRATION */ /*****************************/ MCalibrationQECam qecam; MJCalibration calloop; calloop.SetInput(&caliter); calloop.SetExtractor(&extractor); // // Set the corr. cams: // calloop.SetQECam(qecam); calloop.SetBadPixels(pedloop.GetBadPixels()); // // Apply rel. time calibration: // calloop.SetRelTimeCalibration(); // Use as arrival time extractor MArrivalTimeCalc2: // calloop.SetArrivalTimeLevel(2); cout << "***************************" << endl; cout << "** COMPUTING CALIBRATION **" << endl; cout << "***************************" << endl; if (!calloop.Process(pedloop.GetPedestalCam())) return; /************************************************************************/ /* THIRD LOOP: PEDESTAL CALIBRATION INTO PHOTONS */ /************************************************************************/ MParList plist3; MTaskList tlist3; plist3.AddToList(&tlist3); // containers MCerPhotEvt nphot; MPedPhotCam nphotrms; plist3.AddToList(&geomcam); plist3.AddToList(&pedloop.GetPedestalCam()); plist3.AddToList(&calloop.GetCalibrationCam()); plist3.AddToList(&calloop.GetQECam()); plist3.AddToList(&calloop.GetRelTimeCam()); plist3.AddToList(&calloop.GetBadPixels()); plist3.AddToList(&nphot); plist3.AddToList(&nphotrms); //tasks MReadMarsFile read3("Events"); read3.DisableAutoScheme(); static_cast(read3).AddFiles(pediter); MPedPhotCalc photrmscalc; tlist3.AddToList(&read3); tlist3.AddToList(&geomapl); tlist3.AddToList(&extractor); tlist3.AddToList(&calibrate); tlist3.AddToList(&photrmscalc); // Create and setup the eventloop MEvtLoop evtloop3; evtloop3.SetParList(&plist3); if (!evtloop3.Eventloop()) return; tlist3.PrintStatistics(); /************************************************************************/ /* FOURTH LOOP: DATA CALIBRATION INTO PHOTONS */ /************************************************************************/ MParList plist4; MTaskList tlist4; plist4.AddToList(&tlist4); // containers MHillas hillas; MSrcPosCam srcposcam; MRawRunHeader runhead; plist4.AddToList(&geomcam); plist4.AddToList(&pedloop.GetPedestalCam()); plist4.AddToList(&calloop.GetCalibrationCam()); plist4.AddToList(&calloop.GetQECam()); plist4.AddToList(&calloop.GetRelTimeCam()); plist4.AddToList(&calloop.GetBadPixels()); plist4.AddToList(&nphot); plist4.AddToList(&nphotrms); plist4.AddToList(&srcposcam); plist4.AddToList(&hillas); plist4.AddToList(&runhead); //tasks MReadMarsFile read4("Events"); read4.DisableAutoScheme(); static_cast(read4).AddFiles(datiter); MImgCleanStd clean(lcore,ltail); MHillasCalc hcalc; MHillasSrcCalc csrc1; MWriteRootFile write(outname,"RECREATE"); write.AddContainer("MHillas" , "Parameters"); write.AddContainer("MHillasSrc" , "Parameters"); write.AddContainer("MHillasExt" , "Parameters"); write.AddContainer("MNewImagePar" , "Parameters"); write.AddContainer("MRawEvtHeader" , "Parameters"); write.AddContainer("MRawRunHeader" , "Parameters"); write.AddContainer("MConcentration" , "Parameters"); write.AddContainer("MSrcPosCam", "Parameters"); tlist4.AddToList(&read4); tlist4.AddToList(&geomapl); tlist4.AddToList(&extractor); tlist4.AddToList(&calibrate); tlist4.AddToList(&clean); tlist4.AddToList(&hcalc); tlist4.AddToList(&csrc1); tlist4.AddToList(&write); // Create and setup the eventloop MEvtLoop datloop; datloop.SetParList(&plist4); cout << "*************************************************************" << endl; cout << "*** COMPUTING DATA USING EXTRACTED SIGNAL (IN PHOTONS) ***" << endl; cout << "*************************************************************" << endl; if (!datloop.Eventloop(nmaxevents)) return; tlist4.PrintStatistics(); } /******************************************************************************/ Bool_t readDatacards(TString& filename) { ifstream ifun(filename.Data()); if(!ifun) { cout << "File " << filename << " not found" << endl; return kFALSE; } TString word; while(ifun >> word) { // skip comments if(word[0]=='/' && word[1]=='/') { while(ifun.get()!='\n'); // skip line continue; } // number of events if(strcmp(word.Data(),"NEVENTS")==0) ifun >> nmaxevents; // input file directory if(strcmp(word.Data(),"IDIR")==0) { if(idirname.Length()) cout << "readDataCards Warning: overriding input directory file name" << endl; ifun >> idirname; } // pedestal runs if(strcmp(word.Data(),"PRUNS")==0) { if(pediter.GetNumRuns()) cout << "readDataCards Warning: adding pedestal runs to the existing list" << endl; ifun >> word; pediter.AddRuns(word.Data(),idirname.Data()); } // calibration runs if(strcmp(word.Data(),"CRUNS")==0) { if(caliter.GetNumRuns()) cout << "readDataCards Warning: adding calibration runs to the existing list" << endl; ifun >> word; caliter.AddRuns(word.Data(),idirname.Data()); } // data runs if(strcmp(word.Data(),"DRUNS")==0) { if(datiter.GetNumRuns()) cout << "readDataCards Warning: adding data runs to the existing list" << endl; ifun >> word; datiter.AddRuns(word.Data(),idirname.Data()); } // output file name if(strcmp(word.Data(),"OUTFILE")==0) { if(outname.Length()) cout << "readDataCards Warning: overriding output file name" << endl; ifun >> outname; } // calibration flag if(strcmp(word.Data(),"CALFLAG")==0) ifun >> calflag; // cleaning level if(strcmp(word.Data(),"CLEANLEVEL")==0) { ifun >> lcore; ifun >> ltail; } } pediter.Reset(); caliter.Reset(); datiter.Reset(); TString pfile; // Dump read values cout << "************************************************" << endl; cout << "* Datacards read from file " << filename << endl; cout << "************************************************" << endl; cout << "Pedestal file (s): " << endl; while(!(pfile=pediter.Next()).IsNull()) cout << pfile << endl; cout << "Calibration file (s): " << endl; while(!(pfile=caliter.Next()).IsNull()) cout << pfile << endl; cout << "Data file (s): " << endl; while(!(pfile=datiter.Next()).IsNull()) cout << pfile << endl; cout << "Maximum number of events: " << nmaxevents << endl; cout << "Output file name: " << outname << endl; cout << "Calibration flag: " << calflag << endl; cout << "Cleaning level: ("<