void ReadSetup(TString fname, MAstroCamera &cam) { MMcConfigRunHeader *config=0; MGeomCam *geom=0; TFile file(fname); TTree *tree = (TTree*)file.Get("RunHeaders"); tree->SetBranchAddress("MMcConfigRunHeader", &config); if (tree->GetBranch("MGeomCam")) tree->SetBranchAddress("MGeomCam", &geom); tree->GetEntry(0); cam.SetMirrors(*config->GetMirrors()); cam.SetGeom(*geom); } //-------------------------------------------------------------------------- // // findTelAxisFromStars // // macro to // - determine the positions of stars in the camera from the DC currents // - and to find the expected position of the source in the camera // void findTelAxisFromStars(const TString filename="20040422_23213_D_Mrk421_E.root", const TString directory="/.magic/magicserv01/MAGIC/rootdata/2004_04_22/", const UInt_t numEvents = 0) { gLog.SetNoColors(); //------------------------------------------------------------------- // book some histograms which are to be filled in the event loop // TH2D *aberr = new TH2D("OpticalAberr", "Opt_Aberr", 100, 0.0, 650.0, 100, 0.0, 650.0); aberr->SetXTitle("Ideal distance from center [mm]"); aberr->SetYTitle("Expected distance from center [mm]"); TProfile *reldiff = new TProfile("RelativeDiff", "Rel_Diff", 100, 0.0, 650.0, -0.2, 0.2, "S"); reldiff->SetXTitle("Ideal distance from center [mm]"); reldiff->SetYTitle("(R - R0)/R0"); //------------------------------------------------------------------- MParList plist; MTaskList tlist; plist.AddToList(&tlist); MGeomCamMagic geomcam; MCameraDC dccam; MStarCam starcam; //$$$$$$$$$$$$$$$$ ww MObservatory obs; MStarCam sourcecam; sourcecam.SetName("MSourceCam"); MHTelAxisFromStars htelaxis; plist.AddToList(&obs); plist.AddToList(&htelaxis); //$$$$$$$$$$$$$$$$ ww plist.AddToList(&geomcam); plist.AddToList(&dccam); plist.AddToList(&starcam); plist.AddToList(&sourcecam); // Reads the trees of the root file and the analysed branches MReadReports read; read.AddTree("Currents"); read.AddTree("Drive"); // If you do not include Drive info, the MFindStars class // will not use the star catalog method read.AddFile(directory+filename); // after the reading of the trees!!! read.AddToBranchList("MReportCurrents.*"); read.AddToBranchList("MReportDrive.*"); MGeomApply geomapl; // TString continuoslightfile = // "/data/MAGIC/Period016/rootdata/2004_04_16/20040416_22368_P_Off3c279-2CL100_E.root"; TString continuoslightfile = "/.magic/magicserv01/MAGIC/rootdata/2004_04_16/20040416_22368_P_Off3c279-2CL100_E.root"; Float_t mindc = 0.9; //[uA] MCalibrateDC dccal; dccal.SetFileName(continuoslightfile); dccal.SetMinDCAllowed(mindc); const Int_t numblind = 5; const Short_t x[numblind] = { 47, 124, 470, 475, 571}; const TArrayS blindpixels(numblind,(Short_t*)x); Float_t ringinterest = 100; //[mm] Float_t tailcut = 2.5; UInt_t integratedevents = 1; // We need the MAGIC mirror geometry from a MC header: //TString geometryfile = "/mcdata/standard/camera/NSB_013/Gamma/Gamma_zbin9_90_7_1480to1489_w0.root"; TString geometryfile = "/.magic/magicserv01/MAGIC/mcdata/standard/camera/NSB_013/Gamma/Gamma_zbin9_90_7_1480to1489_w0.root"; // We need the Bright Star Catalog: //TString catalogfile = "/home/rwagner/bsc5.dat"; TString catalogfile = "mtemp/mmpi/macros/bsc5.dat"; //$$$$$$$$$$$$$$$$ ww MSourceDirections sdirs; sdirs.SetGeometryFile(geometryfile); const Double_t ra = MAstro::Hms2Rad(11, 4, 26); const Double_t dec = MAstro::Dms2Rad(38, 12, 36); sdirs.SetRaDec(ra,dec); sdirs.AddDirection(ra,dec,1,"Mkn 421"); const Double_t ra = MAstro::Hms2Rad(11, 0, 50); const Double_t dec = MAstro::Dms2Rad(39, 12, 44); sdirs.AddDirection(ra,dec,1,"My_UMa 49"); const Double_t ra = MAstro::Hms2Rad(11, 4, 31); const Double_t dec = MAstro::Dms2Rad(38, 14, 29); sdirs.AddDirection(ra,dec,1,"My_UMa 51"); // dummy source dirtections const Double_t ra = MAstro::Hms2Rad(11, 4, 31); const Double_t dec = MAstro::Dms2Rad(38, 20, 0); sdirs.AddDirection(ra,dec,1,"Dummy1"); const Double_t ra = MAstro::Hms2Rad(11, 4, 31); const Double_t dec = MAstro::Dms2Rad(38, 20, 0); sdirs.AddDirection(ra,dec,1,"Dummy2"); const Double_t ra = MAstro::Hms2Rad(11, 4, 31); const Double_t dec = MAstro::Dms2Rad(38, 30, 0); sdirs.AddDirection(ra,dec,1,"Dummy3"); const Double_t ra = MAstro::Hms2Rad(11, 4, 31); const Double_t dec = MAstro::Dms2Rad(38, 40, 0); sdirs.AddDirection(ra,dec,1,"Dummy4"); const Double_t ra = MAstro::Hms2Rad(11, 4, 31); const Double_t dec = MAstro::Dms2Rad(38, 50, 0); sdirs.AddDirection(ra,dec,1,"Dummy5"); const Double_t ra = MAstro::Hms2Rad(11, 4, 31); const Double_t dec = MAstro::Dms2Rad(39, 0, 0); sdirs.AddDirection(ra,dec,1,"Dummy6"); const Double_t ra = MAstro::Hms2Rad(11, 4, 31); const Double_t dec = MAstro::Dms2Rad(39, 10, 0); sdirs.AddDirection(ra,dec,1,"Dummy7"); const Double_t ra = MAstro::Hms2Rad(11, 4, 31); const Double_t dec = MAstro::Dms2Rad(39, 20, 0); sdirs.AddDirection(ra,dec,1,"Dummy8"); const Double_t ra = MAstro::Hms2Rad(11, 4, 31); const Double_t dec = MAstro::Dms2Rad(39, 30, 0); sdirs.AddDirection(ra,dec,1,"Dummy9"); const Double_t ra = MAstro::Hms2Rad(11, 4, 31); const Double_t dec = MAstro::Dms2Rad(39, 40, 0); sdirs.AddDirection(ra,dec,1,"Dummy10"); const Double_t ra = MAstro::Hms2Rad(11, 4, 31); const Double_t dec = MAstro::Dms2Rad(39, 50, 0); sdirs.AddDirection(ra,dec,1,"Dummy11"); const Double_t ra = MAstro::Hms2Rad(11, 4, 31); const Double_t dec = MAstro::Dms2Rad(40, 0, 0); sdirs.AddDirection(ra,dec,1,"Dummy12"); const Double_t ra = MAstro::Hms2Rad(11, 4, 31); const Double_t dec = MAstro::Dms2Rad(40, 10, 0); sdirs.AddDirection(ra,dec,1,"Dummy13"); const Double_t ra = MAstro::Hms2Rad(11, 4, 31); const Double_t dec = MAstro::Dms2Rad(40, 20, 0); sdirs.AddDirection(ra,dec,1,"Dummy14"); sdirs.SetRadiusFOV(3); //$$$$$$$$$$$$$$$$ ww MFindStars findstars; findstars.SetBlindPixels(blindpixels); findstars.SetRingInterest(ringinterest); findstars.SetDCTailCut(tailcut); findstars.SetNumIntegratedEvents(integratedevents); findstars.SetMinuitPrintOutLevel(-1); findstars.SetGeometryFile(geometryfile); findstars.SetBSCFile(catalogfile); const Double_t ra = MAstro::Hms2Rad(11, 4, 26); const Double_t dec = MAstro::Dms2Rad(38, 12, 36); findstars.SetRaDec(ra,dec); findstars.SetLimMag(8); findstars.SetRadiusFOV(1.5); //$$$$$$$$$$$$$$$$ ww Int_t InputType = 0; //findstars.SetUseCorrelatedGauss(kFALSE); MTelAxisFromStars telaxis; //telaxis.FixRotationAngleAt(-1.0); //telaxis.FixScaleFactorAt(-1.0); telaxis.SetInputType(InputType); MFillH fillhisto("MHTelAxisFromStars[MHTelAxisFromStars]",""); htelaxis.SetInputType(InputType); //$$$$$$$$$$$$$$$$ ww tlist.AddToList(&geomapl); tlist.AddToList(&read); tlist.AddToList(&dccal); tlist.AddToList(&findstars, "Currents"); //$$$$$$$$$$$$$$$$ ww tlist.AddToList(&sdirs); tlist.AddToList(&telaxis); tlist.AddToList(&fillhisto); //$$$$$$$$$$$$$$$$ ww // The following lines you only need if in addition you want to display // independent MAstroCamera output // // TString fname = "/mcdata/standard/camera/NSB_013/Gamma/Gamma_zbin9_90_7_1480to1489_w0.root"; // MObservatory magic1; // const Double_t ra = MAstro::Hms2Rad(11, 4, 26); //Mkn421 // const Double_t dec = MAstro::Dms2Rad(38, 12, 36); // // MAstroCamera stars; // ReadSetup(fname, stars); // stars.SetLimMag(9); // stars.SetRadiusFOV(3); // stars.SetRaDec(ra, dec); // stars.ReadBSC("/home/rwagner/bsc5.dat"); // stars.SetObservatory(magic1); MEvtLoop evtloop; evtloop.SetParList(&plist); if (!evtloop.PreProcess()) return; MHCamera display(geomcam); display.SetPrettyPalette(); display.Draw(); gPad->cd(1); starcam.Draw(); UInt_t numevents=0; while (tlist.Process()) { //gLog << "---------------------------------------------" << endl; //gLog << "Macro : MStarPos content = " << endl; //starcam.Print("namepossizchierr"); //gLog << "---------------------------------------------" << endl; //gLog << "Macro : MSourcePos content = " << endl; //sourcecam.Print("namepossizchierr"); //gLog << "---------------------------------------------" << endl; numevents++; if (numevents%integratedevents==0) { display.SetCamContent(findstars.GetDisplay()); gPad->Modified(); gPad->Update(); // This line prints the results: // starcam.Print(); // This is how to access the TList of stars: // TList* starlist = starcam.GetList(); // This is how to iterate over stars found: // TIter Next(starlist); // MStarPos* star; // UInt_t starnum = 0; // cout << filename << " "; // cout << "Iterating over list" << endl; // while ( (star=(MStarPos*)Next()) ) // { // cout << "star[" << starnum << "] "; // cout << star->GetMeanX() << " " // << star->GetMeanY() << " "; // starnum++; // } // cout << endl; } // fill the histograms // loop over stars in "MSourceCam" TList* sourcelist = sourcecam.GetList(); TIter Nextsource(sourcelist); MStarPos* source; while ( (source=(MStarPos*)Nextsource()) ) { Double_t R = sqrt( source->GetXExp()*source->GetXExp() + source->GetYExp()*source->GetYExp() ); Double_t R0= sqrt( source->GetXIdeal()*source->GetXIdeal() + source->GetYIdeal()*source->GetYIdeal() ); //cout << "R0, R = " << R0 << ", " << R << endl; aberr->Fill(R0, R, 1.0); reldiff->Fill(R0, (R-R0)/R0), 1.0); } //-------------------------------------- MTime time; time.Set(2004, 4, 22, 21, 51, 15); //superimpose star picture // stars.SetTime(time); // TObject *o = stars.Clone(); // o->SetBit(kCanDelete); // o->Draw(); // wait after each event //if (!HandleInput()) // break; } evtloop.PostProcess(); tlist.PrintStatistics(); plist.Print(); //$$$$$$$$$$$$$$$$ ww //gLog << "Event loop finished; call DrawClone of MHTelAxisFromStars" << endl; //TObject *srccam = plist.FindObject("MSourceCam"); //gLog << "srccam = " << srccam << endl; TObject *obj = plist.FindObject("MHTelAxisFromStars"); if (obj) { //obj->Print(); //obj->Dump(); obj->DrawClone(); } else gLog << "address of MHTelAxisFromStars container is zero" << endl; //$$$$$$$$$$$$$$$$ ww // plot the histograms TCanvas *canv = new TCanvas("Optical_Aberration", "Optical aberration", 0, 0, 600, 300); canv->Divide(2,1); canv->SetBorderMode(0); canv->cd(1); aberr->Draw(); canv->cd(2); reldiff->Draw(); } //---------------------------------------------------------------------- // // Bool_t HandleInput() { TTimer timer("gSystem->ProcessEvents();", 50, kFALSE); while (1) { // // While reading the input process gui events asynchronously // timer.TurnOn(); TString input = Getline("Type 'q' to exit, to go on: "); timer.TurnOff(); if (input=="q\n") return kFALSE; if (input=="\n") return kTRUE; }; return kFALSE; } //=========================================================================