Index: /trunk/MagicSoft/Mars/mtemp/mifae/macros/makeHillas.C =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/macros/makeHillas.C (revision 4310) +++ /trunk/MagicSoft/Mars/mtemp/mifae/macros/makeHillas.C (revision 4310) @@ -0,0 +1,561 @@ +/*********************************/ +/* 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.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.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.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); + + // 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; + + // set bad pixels + MBlindPixelCalc blind; + MBlindPixelCalc blind2; + const Short_t x[16] = {330,395,329,396,389, + 323,388,322,384,385, + 386,387,321,320,319, + 394}; + const TArrayS bp(16,(Short_t*)x); + blind.SetPixelIndices(bp); + blind2.SetPixelIndices(bp); + + 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"); + //tlist4.AddToList(&blind,"Events"); + 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(&srcposcalc,"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: ("<