/*********************************/ /* Compute the hillas parameters */ /*********************************/ #include "TString.h" #include "TArrayS.h" #include "TFile.h" #include "TH1.h" #include "TH2.h" #include "TProfile.h" #include "MArray.h" #include "MParContainer.h" #include "MParList.h" #include "MTaskList.h" #include "MPedestalCam.h" #include "MBadPixelsCam.h" #include "MBadPixelsTreat.h" #include "MReadMarsFile.h" #include "MReadReports.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 "MExtractor.h" #include "MExtractFixedWindow.h" #include "MExtractFixedWindowPeakSearch.h" #include "MExtractSlidingWindow.h" #include "MExtractTimeAndChargeSpline.h" #include "MPedCalcFromLoGain.h" #include "MExtractSignal.h" #include "MCalibrationChargeCalc.h" #include "MFCosmics.h" #include "MContinue.h" #include "MLog.h" #include "MCerPhotEvt.h" #include "MPedPhotCam.h" #include "MArrivalTime.h" #include "MCalibrateData.h" #include "MPedPhotCalc.h" #include "MHillas.h" #include "MNewImagePar.h" #include "MRawRunHeader.h" #include "MSrcPosCam.h" #include "MImgCleanStd.h" #include "MHillasSrcCalc.h" #include "MHillasCalc.h" #include "MArrivalTimeCam.h" #include "MArrivalTimeCalc2.h" #include "MIslands.h" #include "MImgIsland.h" #include "MIslandsCalc.h" #include "MIslandsClean.h" #include "MWriteRootFile.h" #include "MArgs.h" #include "MRunIter.h" #include "MJPedestal.h" #include "MJCalibration.h" #include "MHillasDisplay.h" #include "MF.h" #include "MContinue.h" #include "MReportDrive.h" #include "MTopology.h" #include "MTopologyCalc.h" #include "MRanForest.h" #include "MRanTree.h" #include "MRanForestCalc.h" #include "MRanForestFill.h" #include "MPointingPos.h" #include "MPointingPosCalc.h" #include "MImageParDisp.h" #include "MDispParameters.h" #include "MDispCalc.h" #include "MHDisp.h" #include "MFillH.h" #include "TApplication.h" #include "TClass.h" #include #include #include using namespace std; Bool_t readDatacards(TString& filename); void makeHillas(); // declaration of variables read from datacards TString outname; const TString outpath = "./"; TString idirname; TString filter; TString psfilename("makehillas.ps"); TString bpfilename; MRunIter pedcaliter; MRunIter caliter; MRunIter datiter; UInt_t display = 0; ULong_t nmaxevents= 999999999; Short_t calflag = 1; Bool_t caltimeflag = kFALSE; Short_t cleanflag = 1; UShort_t lrings = 1; Float_t lcore = 3.0; Float_t ltail = 1.5; Int_t islflag = 0; Int_t topflag = 0; Float_t lnew = 40; Int_t kmethod = 1; Int_t kalgorithm = 1; Int_t nfiles = 0; Int_t hifirst = 0; Int_t hilast = 14; Int_t lofirst = 2; Int_t lolast = 14; Int_t wsize = 6; Int_t sext = 0; TString rffilename; TString dispfilename; const TString defaultcard="makehillas.datacard"; char* chext[4]={"Fixed window","Sliding window","Peak Search","Spline"}; /*************************************************************/ static void Usage() { gLog <" << endl << endl; gLog << " : datacards file name (dafault makehillas.datacards)" << endl; gLog << " -?/-h: This help" << endl << endl; } /*************************************************************/ int main(int argc, char **argv) { // create a TApplication to be able to TApplication app("Application",0,0); // to deal correctly with the streamer MArray::Class()->IgnoreTObjectStreamer(); MParContainer::Class()->IgnoreTObjectStreamer(); // evaluate arguments MArgs arg(argc, argv); if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help")) { Usage(); return -1; } TString datacard = arg.GetArgumentStr(0); if(!datacard.Length()) datacard = defaultcard; if(!readDatacards(datacard)) { cout << "Error reading datacards. Stoping" << endl; return -1; } makeHillas(); } /*************************************************************/ void makeHillas() { // Set the signal extractor MExtractor* extractor; switch(sext) { case 0: extractor = new MExtractFixedWindow(); break; case 1: extractor = new MExtractSlidingWindow(); ((MExtractSlidingWindow*)extractor)->SetWindowSize(wsize,wsize); break; case 2: extractor = new MExtractFixedWindowPeakSearch(); ((MExtractFixedWindowPeakSearch*)extractor)->SetWindows(wsize,wsize,4); break; case 3: extractor = new MExtractTimeAndChargeSpline(); ((MExtractTimeAndChargeSpline*)extractor)->SetTimeType(MExtractTimeAndChargeSpline::kHalfMaximum); ((MExtractTimeAndChargeSpline*)extractor)->SetChargeType(MExtractTimeAndChargeSpline::kIntegral); ((MExtractTimeAndChargeSpline*)extractor)->SetRiseTime((Float_t)wsize*0.25); ((MExtractTimeAndChargeSpline*)extractor)->SetFallTime((Float_t)wsize*0.75); break; default: extractor = new MExtractFixedWindow(); break; } extractor->SetRange(hifirst,hilast,lofirst,lolast); /****************************************************/ /* FIRST LOOP: PEDESTAL COMPUTATION FOR CALIBRATION */ /****************************************************/ // If you want to exclude pixels from the beginning, read // an ascii-file with the corr. pixel numbers (see MBadPixelsCam) //badcam.AsciiRead("badpixels.dat"); MBadPixelsCam badcam; MGeomCamMagic geomcam; MGeomApply geomapl; badcam.InitSize(geomcam.GetNumPixels()); if(bpfilename.Length()>0) { ifstream fin("/home/emma/Mars/mtemp/mifae/programs/badpixels.dat"); badcam.AsciiRead((istream&)fin); badcam.Print(); } MJPedestal pedloop1; MJPedestal pedloop2; MJCalibration calloop; if(calflag>0) { pedloop1.SetInput(&pedcaliter); pedloop1.SetExtractor(extractor); pedloop1.SetBadPixels(badcam); pedloop1.SetNoStorage(); if (!pedloop1.Process()) return; if (extractor->InheritsFrom("MExtractTimeAndCharge")) { /***********************************************************/ /* NEEDED FOR SECOND LOOP: EXTRACTOR RESOLUTION COMPUTATION */ /***********************************************************/ pedloop2.SetNoStorage(); pedloop2.SetExtractor(extractor); pedloop2.SetBadPixels(pedloop1.GetBadPixels()); pedloop2.SetPedestals(pedloop1.GetPedestalCam()); pedloop2.SetInput(&pedcaliter); if (!pedloop2.Process()) return; calloop.SetExtractorCam(pedloop2.GetPedestalCam()); } } MPedestalCam& pedcammean = pedloop1.GetPedestalCam(); extractor->SetPedestals(&pedcammean); MPedestalCam pedcamrms; /*****************************/ /* SECOND LOOP: CALIBRATION */ /*****************************/ calloop.SetRelTimeCalibration(caltimeflag); calloop.SetExtractor(extractor); calloop.SetInput(&caliter); calloop.SetBadPixels(pedloop1.GetBadPixels()); if(calflag==2) calloop.SetUseBlindPixel(); if(calflag>0) if (!calloop.Process(pedcammean)) return; MCalibrationChargeCam &chargecam = calloop.GetCalibrationCam(); chargecam.Print(); /************************************************************************/ /* THIRD (SMALL) LOOP TO GET INITIAl PEDESTALS */ /************************************************************************/ MParList plist3; MTaskList tlist3; plist3.AddToList(&tlist3); plist3.AddToList(&geomcam); plist3.AddToList(&pedcammean); //tasks MReadMarsFile read3("Events"); static_cast(read3).AddFiles(datiter); read3.DisableAutoScheme(); MPedCalcFromLoGain pedlo1; pedlo1.SetPedestalUpdate(kFALSE); const Float_t win = extractor->GetNumHiGainSamples(); pedlo1.SetExtractWindow(15, (UShort_t)TMath::Nint(win)); pedlo1.SetNumEventsDump(500); pedlo1.SetMaxSignalVar(40); pedlo1.SetPedestalsOut(&pedcammean); tlist3.AddToList(&read3); tlist3.AddToList(&geomapl); tlist3.AddToList(&pedlo1); // Create and setup the eventloop MEvtLoop evtloop3; evtloop3.SetParList(&plist3); if (!evtloop3.Eventloop(500)) return; tlist3.PrintStatistics(); /************************************************************************/ /* FOURTH LOOP TO READ THE RANDOM FOREST FILE(S) */ /************************************************************************/ MParList plistrf; MTaskList tlistrf; MRanForest ranforest; MRanTree rantree; if(rffilename.Length()) { plistrf.AddToList(&tlistrf); plistrf.AddToList(&ranforest); plistrf.AddToList(&rantree); MReadTree readrf("Tree",rffilename); readrf.DisableAutoScheme(); MRanForestFill rffill; rffill.SetNumTrees(100); tlistrf.AddToList(&readrf); tlistrf.AddToList(&rffill); MEvtLoop evtlooprf; evtlooprf.SetParList(&plistrf); if (!evtlooprf.Eventloop()) return; tlistrf.PrintStatistics(); } /************************************************************************/ /* FIFTH LOOP: PEDESTAL+DATA CALIBRATION INTO PHOTONS */ /************************************************************************/ MParList plist5; MTaskList tlist5; plist5.AddToList(&tlist5); //tasks // MReadMarsFile read5("Events"); MReadReports read5; read5.AddTree("Drive"); read5.AddTree("Events","MTime.",kTRUE); read5.AddFiles(datiter); read5.AddToBranchList("MReportDrive.*"); read5.AddToBranchList("MRawEvtData.*"); read5.AddToBranchList("MRawEvtHeader.*"); // read5.DisableAutoScheme(); // static_cast(read5).AddFiles(datiter); tlist5.AddToList(&read5); tlist5.AddToList(&geomapl); tlist5.AddToList(&pedlo1); pedlo1.SetPedestalUpdate(kTRUE); MCalibrateData::CalibrationMode_t calMode=MCalibrateData::kDefault; if(calflag==0) calMode=MCalibrateData::kNone; if(calflag==-1) calMode=MCalibrateData::kDummy; MCalibrateData photcalc; photcalc.SetCalibrationMode(calMode); photcalc.EnablePedestalType(MCalibrateData::kRun); photcalc.EnablePedestalType(MCalibrateData::kEvent); photcalc.SetPedestalCamMean(&pedcammean); if (extractor->InheritsFrom("MExtractTimeAndCharge")) { /************************************************************************/ /* SIXTH (SMALL) LOOP TO GET INITIAl PEDESTALS RMS */ /************************************************************************/ MParList plist4; MTaskList tlist4; plist4.AddToList(&tlist4); plist4.AddToList(&geomcam); plist4.AddToList(&pedcamrms); //tasks MReadMarsFile read4("Events"); static_cast(read4).AddFiles(datiter); read4.DisableAutoScheme(); MPedCalcFromLoGain pedlo2; pedlo2.SetPedestalUpdate(kFALSE); pedlo2.SetExtractor((MExtractTimeAndCharge*)extractor); pedlo2.SetNumEventsDump(500); pedlo2.SetMaxSignalVar(40); pedlo2.SetPedestalsIn(&pedcammean); pedlo2.SetPedestalsOut(&pedcamrms); tlist4.AddToList(&read4); tlist4.AddToList(&geomapl); tlist4.AddToList(&pedlo2); // Create and setup the eventloop MEvtLoop evtloop4; evtloop4.SetParList(&plist4); if (!evtloop4.Eventloop(500)) return; tlist4.PrintStatistics(); tlist5.AddToList(&pedlo2); pedlo2.SetPedestalUpdate(kTRUE); photcalc.SetPedestalCamRms(&pedcamrms); } // containers MCerPhotEvt nphot; MPedPhotCam nphotrms; MArrivalTime arrtime; MHillas hillas; MNewImagePar newimagepar; MSrcPosCam source; MRawRunHeader runhead; MArrivalTimeCam timecam; MReportDrive reportdrive; // islands MIslands isl; MIslands isl2; MIslands isl3; MTopology topology; // disp // (read in optimum Disp parameter values) MDispParameters dispin; TArrayD dispPar; if(dispfilename.Length()) { TFile inparam(dispfilename); dispin.Read("MDispParameters"); cout << "Optimum parameter values taken for calculating Disp : " << endl; dispPar = dispin.GetParameters(); for (Int_t i=0; iSetIslandsName("MIslands"); if (islflag == 2) disphillas->SetIslandsName("MIslands2"); if (islflag == 3) disphillas->SetIslandsName("MIslands3"); } tlist5.AddToList(extractor); tlist5.AddToList(&timecalc); tlist5.AddToList(&photcalc); if(calflag==11 || calflag==21) tlist5.AddToList(&interpolatebadpixels); tlist5.AddToList(&clean); tlist5.AddToList(&island); if (islflag == 2) { tlist5.AddToList(&islclean); tlist5.AddToList(&island2); } if (islflag == 3) { tlist5.AddToList(&islclean); tlist5.AddToList(&island3); } tlist5.AddToList(&hcalc); tlist5.AddToList(&csrc1); if(filter.Length()) tlist5.AddToList(&applycut); if(topflag>0) tlist5.AddToList(&topcalc); if(rffilename.Length()) tlist5.AddToList(&hadrcalc); tlist5.AddToList(&pointingposcalc); if(dispfilename.Length()) { tlist5.AddToList(&dispcalc); tlist5.AddToList(&filldisp); } tlist5.AddToList(&write); if(display) { disphillas->SetPSFile(); disphillas->SetPSFileName(psfilename); if(display==2) disphillas->SetPause(kFALSE); tlist5.AddToList(disphillas); } // Create and setup the eventloop MEvtLoop datloop; datloop.SetParList(&plist5); cout << "*************************************************************" << endl; cout << "*** COMPUTING DATA USING EXTRACTED SIGNAL (IN PHOTONS) ***" << endl; cout << "*************************************************************" << endl; if (!datloop.Eventloop(nmaxevents)) return; tlist5.PrintStatistics(); delete extractor; } //------------------------------------------------------------------------------- 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 for calibration if(strcmp(word.Data(),"PCRUNS")==0) { if(pedcaliter.GetNumRuns()) cout << "readDataCards Warning: adding pedestal runs for calibration to the existing list" << endl; ifun >> word; pedcaliter.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; } // exclusion cut if(strcmp(word.Data(),"FILTER")==0) { if(filter.Length()) cout << "readDataCards Warning: overriding existing cut" << endl; char ch; while((ch=ifun.get())!='\n') filter.Append(ch); } // display flag if(strcmp(word.Data(),"DISPLAY")==0) ifun >> display; // ps file name if(strcmp(word.Data(),"PSFILENAME")==0) ifun >> psfilename; // excluded pixels file name if(strcmp(word.Data(),"BADPIXELFILE")==0) ifun >> bpfilename; // calibration flag if(strcmp(word.Data(),"CALFLAG")==0) ifun >> calflag; // calibration flag if(strcmp(word.Data(),"CALTIME")==0) ifun >> caltimeflag; // cleaning level if(strcmp(word.Data(),"CLEANLEVEL")==0) { ifun >> lcore; ifun >> ltail; if(ifun.get()!='\n'){ ifun.unget(); ifun >> lrings; if(ifun.get()!='\n'){ ifun.unget(); ifun >> cleanflag; } } } // random forest file name if(strcmp(word.Data(),"RANFOREST")==0) { if(rffilename.Length()) cout << "readDataCards Warning: overriding random forest file name" << endl; ifun >> rffilename; } // disp file name if(strcmp(word.Data(),"DISP")==0) { if(dispfilename.Length()) cout << "readDataCards Warning: overriding disp file name" << endl; ifun >> dispfilename; } // cleaning level if(strcmp(word.Data(),"EXTRACTOR")==0) ifun >> sext >> hifirst >> hilast >> lofirst >> lolast >> wsize; if(strcmp(word.Data(),"TOPFLAG")==0) ifun >> topflag; if(strcmp(word.Data(),"ISLFLAG")==0) { ifun >> islflag; // if (islflag == 1 || islflag == 2) ifun >> kalgorithm; } // island cleaning if (islflag == 2){ if(strcmp(word.Data(),"ISLANDCLEAN")==0) { ifun >> kmethod; ifun >> lnew; } } } pedcaliter.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) for calibration: " << endl; while(!(pfile=pedcaliter.Next()).IsNull()) cout << pfile << endl; cout << "Calibration file (s): " << endl; while(!(pfile=caliter.Next()).IsNull()) cout << pfile << endl; cout << "Warning: Pedestals for data will be computed from data themselves" << endl; cout << "Data file (s): " << endl; while(!(pfile=datiter.Next()).IsNull()) cout << pfile << endl; cout << "Maximum number of events: " << nmaxevents << endl; if(filter.Length()) cout << "Applying selection cut: " << filter << endl; cout << "Output file name: " << outname << endl; if(display) cout << "Generating PS file: " << psfilename << endl; if(bpfilename.Length()) cout << "Bad pixels will be read from " << bpfilename << endl; cout << "Signal Extractor: " << chext[sext] << " with bounds ("< 0) cout << "Topology parameters will be computed" << endl; if (islflag == 2) { cout << "Island Cleaning: "<< kmethod <<" method "<< lnew << " new threshold" << endl; } cout << "***********" << endl << endl; if(!pedcaliter.GetNumEntries() && calflag>0) { cout << "No pedestal file for calibration specified" << endl; return kFALSE; } if(!caliter.GetNumEntries() && calflag>0) { cout << "No calibration file name specified" << endl; return kFALSE; } if(!datiter.GetNumEntries()) { cout << "No data file name specified" << endl; return kFALSE; } if(!outname.Length()) { cout << "No output file name specified" << endl; return kFALSE; } return kTRUE; }