/*********************************/ /* Compute the hillas parameters */ /*********************************/ Bool_t readDatacards(TString& filename); // initial and final time slices to be used in signal extraction const Byte_t hifirst = 0; const Byte_t hilast = 13; const Byte_t lofirst = 3; const Byte_t lolast = 12; // declaration of variables read from datacards TString outname; TString idirname; MRunIter pedcaliter; 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 islflag = 0; Float_t lnew = 40; Int_t kmethod = 1; Int_t nfiles = 0; const TString defaultcard="input.datacard"; //Slow control varialbles Short_t slowflag=0; Bool_t isSlowControlAvailable = kFALSE; TString hvConfigFile="/mnt/users/jlopez/Mars/Files4Mars/Config/HVSettings_FF35q.conf"; TString continuosLightFile="/local_disk/rootdata/Miscellaneous/Period016/2004_04_16/dc_2004_04_16_04_46_18_22368_Off3c279-2CL100.root"; /*************************************************************/ int makeHillas(TString datacard = defaultcard) { if(!datacard.Length()) datacard = defaultcard; if(!readDatacards(datacard)) { cout << "Error reading datacards. Stoping" << endl; return -1; } // Set the general tasks/containers MExtractFixedWindow extractor; extractor.SetRange(hifirst,hilast,lofirst,lolast); MCalibrationQECam qecam; MBadPixelsCam badcam; MGeomCamMagic geomcam; MGeomApply geomapl; cout << endl; cout << "/****************************************************/" << endl; cout << "/* FIRST LOOP: PEDESTAL FOR CALIBRATION COMPUTATION */" << endl; cout << "/****************************************************/" << endl; cout << endl; // If you want to exclude pixels from the beginning, read // an ascii-file with the corr. pixel numbers (see MBadPixelsCam) ifstream fin("badpixels.dat"); badcam.AsciiRead(fin); fin.close(); MJPedestal pedloop0; pedloop0.SetOutputPath("./"); pedloop0.SetInput(&pedcaliter); pedloop0.SetExtractor(&extractor); pedloop0.SetBadPixels(badcam); if (!pedloop0.Process()) return; cout << endl; cout << "/*****************************/" << endl; cout << "/* SECOND LOOP: CALIBRATION */" << endl; cout << "/*****************************/" << endl; cout << endl; MJCalibration calloop; calloop.SetOutputPath("./"); calloop.SetExtractor(&extractor); calloop.SetInput(&caliter); calloop.SetQECam(qecam); calloop.SetBadPixels(pedloop0.GetBadPixels()); if (!calloop.Process(pedloop0.GetPedestalCam())) return; cout << endl; cout << "/*************************************************/" << endl; cout << "/* THIRD LOOP: PEDESTAL CALIBRATION INTO PHOTONS */" << endl; cout << "/*************************************************/" << endl; cout << endl; // First Compute the pedestals MJPedestal pedloop; pedloop.SetOutputPath("./"); pedloop.SetInput(&pediter); if (!pedloop.Process()) return; MPedestalCam pedcam = pedloop.GetPedestalCam(); // Now convert this pedestals in photons MParList plist3; MTaskList tlist3; plist3.AddToList(&tlist3); // containers MCerPhotEvt nphot; MPedPhotCam nphotrms; MExtractedSignalCam sigcam; plist3.AddToList(&geomcam); plist3.AddToList(&pedcam); plist3.AddToList(&calloop.GetCalibrationCam()); plist3.AddToList(&qecam); plist3.AddToList(&badcam); plist3.AddToList(&sigcam); plist3.AddToList(&nphot); plist3.AddToList(&nphotrms); MCalibrate::CalibrationMode_t calMode=MCalibrate::kDefault; if(calflag==0) calMode=MCalibrate::kNone; if(calflag==-1) calMode=MCalibrate::kDummy; //tasks MReadMarsFile read3("Events"); static_cast(read3).AddFiles(pediter); read3.DisableAutoScheme(); // MReadReports read3; // read3.AddTree("Events","MTime.",kFALSE); // static_cast(read3).AddFiles(pediter); // MExtractSignal extsig; // extsig.SetRange(hifirst,hilast,lofirst,lolast); MCalibrate photcalc(calMode); MPedPhotCalc photrmscalc; tlist3.AddToList(&read3); tlist3.AddToList(&geomapl); tlist3.AddToList(&extractor); tlist3.AddToList(&photcalc); tlist3.AddToList(&photrmscalc); // Create and setup the eventloop MEvtLoop evtloop3; evtloop3.SetParList(&plist3); if (!evtloop3.Eventloop()) return; tlist3.PrintStatistics(); cout << "/**********************************************/" << endl; cout << "/* FOURTH LOOP: DATA CALIBRATION INTO PHOTONS */" << endl; cout << "/**********************************************/" << endl; MParList plist4; MTaskList tlist4; plist4.AddToList(&tlist4); // containers MHillas hillas; MSrcPosCam source; MRawRunHeader runhead; MArrivalTimeCam timecam; MIslands isl; isl.SetName("MIslands1"); MIslands isl2; isl2.SetName("MIslands2"); MReportDrive drive; MStarLocalCam starcam; if (islflag == 1 || islflag == 2) { plist4.AddToList(&timecam); plist4.AddToList(&isl); } if (islflag == 2) { plist4.AddToList(&isl2); } if (isSlowControlAvailable) { plist4.AddToList(&starcam); } plist4.AddToList(&geomcam); plist4.AddToList(&pedcam); plist4.AddToList(&calloop.GetCalibrationCam()); plist4.AddToList(&qecam); plist4.AddToList(&badcam); plist4.AddToList(&nphot); plist4.AddToList(&nphotrms); plist4.AddToList(&source); plist4.AddToList(&hillas); plist4.AddToList(&runhead); //tasks // MReadMarsFile read4("Events"); // static_cast(read4).AddFiles(datiter); // read4.DisableAutoScheme(); MReadReports read4; read4.AddTree("Events","MTime.",kTRUE); read4.AddTree("Currents"); read4.AddTree("Camera"); read4.AddTree("Drive"); static_cast(read4).AddFiles(datiter); read4.AddToBranchList("MReportCurrents.*"); // Task that use Slow control information MFHVNotNominal fhvnotnominal; fhvnotnominal.SetHVNominalValues(hvConfigFile); fhvnotnominal.SetMaxNumPixelsDeviated(40); MContinue hvnotnominal(&fhvnotnominal); MCalibrateDC dccalib; dccalib.SetFileName(continuosLightFile); const Int_t numblind = 5; const Short_t w[numblind] = { 47, 124, 470, 475, 571}; const TArrayS blindpixels(numblind,(Short_t*)w); Float_t ringinterest = 100; //[mm] Float_t tailcut = 3.5; UInt_t integratedevents = 4; MFindStars findstars; findstars.SetBlindPixels(blindpixels); findstars.SetRingInterest(ringinterest); findstars.SetDCTailCut(tailcut); findstars.SetNumIntegratedEvents(integratedevents); findstars.SetMinuitPrintOutLevel(-1); MSrcPosFromStars srcposfromstar; MBadPixelsTreat interpolatebadpixels; interpolatebadpixels.SetUseInterpolation(); MImgCleanStd clean(lcore,ltail); MArrivalTimeCalc2 timecalc; MIslandCalc island; island.SetOutputName("MIslands1"); MIslandClean islclean(lnew); islclean.SetInputName("MIslands1"); islclean.SetMethod(kmethod); MIslandCalc island2; island2.SetOutputName("MIslands2"); MHillasCalc hcalc; MHillasSrcCalc csrc1; MWriteRootFile write(outname,"RECREATE"); write.AddContainer("MTime" , "Parameters"); 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"); if (islflag == 1 || islflag == 2) write.AddContainer("MIslands1" , "Parameters"); if (islflag == 2) write.AddContainer("MIslands2" , "Parameters"); if (isSlowControlAvailable) { write.AddContainer("MStarLocalCam" , "Parameters"); write.AddContainer("MReportDrive" , "Parameters"); } tlist4.AddToList(&read4); tlist4.AddToList(&geomapl); if (isSlowControlAvailable) { tlist4.AddToList(&hvnotnominal,"Events"); tlist4.AddToList(&dccalib,"Currents"); tlist4.AddToList(&findstars,"Currents"); tlist4.AddToList(&srcposfromstar,"Events"); } tlist4.AddToList(&extractor,"Events"); tlist4.AddToList(&photcalc,"Events"); if(calflag==11) tlist4.AddToList(&interpolatebadpixels); tlist4.AddToList(&clean,"Events"); if (islflag == 1 || islflag == 2) { tlist4.AddToList(&timecalc,"Events"); tlist4.AddToList(&island,"Events"); } if (islflag == 2) { tlist4.AddToList(&islclean,"Events"); tlist4.AddToList(&island2,"Events"); } //tlist4.AddToList(&blind2,"Events"); tlist4.AddToList(&hcalc,"Events"); tlist4.AddToList(&csrc1,"Events"); tlist4.AddToList(&write,"Events"); // Create and setup the eventloop MEvtLoop datloop; datloop.SetParList(&plist4); cout << endl; cout << "******************************************************" << endl; cout << "* COMPUTING DATA USING EXTRACTED SIGNAL (IN PHOTONS) *" << endl; cout << "******************************************************" << 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 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()); } // pedestal runs for data if(strcmp(word.Data(),"PRUNS")==0) { if(pediter.GetNumRuns()) cout << "readDataCards Warning: adding pedestal runs for data to the existing list" << endl; ifun >> word; pediter.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; } if(strcmp(word.Data(),"ISLFLAG")==0) { ifun >> islflag; } // island cleaning if (islflag == 2){ if(strcmp(word.Data(),"ISLANDCLEAN")==0) { ifun >> kmethod; ifun >> lnew; } } // slow control data available if(strcmp(word.Data(),"SLOWCRT")==0) { ifun >> slowflag; if (slowflag == 1) isSlowControlAvailable = kTRUE; } // hv configuration file if(strcmp(word.Data(),"HVCONFFILE")==0) { ifun >> hvConfigFile; } // countinous light file to calibrate the dc data if(strcmp(word.Data(),"CLFILE")==0) { ifun >> continuosLightFile; } } pedcaliter.Reset(); 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) for calibration: " << endl; while(kTRUE) { pfile=pedcaliter.Next(); if (pfile.IsNull()) break; cout << pfile << endl; } cout << "Calibration file (s): " << endl; while(kTRUE) { pfile=caliter.Next(); if (pfile.IsNull()) break; cout << pfile << endl; } cout << "Pedestal file (s) for data: " << endl; while(kTRUE) { pfile=pediter.Next(); if (pfile.IsNull()) break; cout << pfile << endl; } cout << "Data file (s): " << endl; while(kTRUE) { pfile=datiter.Next(); if (pfile.IsNull()) break; cout << pfile << endl; } cout << "Maximum number of events: " << nmaxevents << endl; cout << "Output file name: " << outname << endl; cout << "Calibration flag: " << calflag << endl; cout << "Cleaning level: ("<