Ignore:
Timestamp:
01/20/05 10:09:24 (20 years ago)
Author:
domingo
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mtemp/mifae/library/MHDisp.cc

    r5671 r5904  
    2929//                                                                         //
    3030//   container holding the histograms for Disp                             //
    31 //   also computes the Chi^2 of the Disp optimization                      //
     31//   also computes the minimization parameter of the Disp optimization     //
    3232//                                                                         //
    3333/////////////////////////////////////////////////////////////////////////////
    … …  
    5454
    5555#include "MHMatrix.h"
     56#include "MParameters.h"
    5657
    5758#include "MParList.h"
    … …  
    6364
    6465using namespace std;
     66
     67enum dispVar_t {kX,kY,kMeanX,kMeanY,kDelta,kSize,kM3Long,kAsym,
     68                kEnergy,kImpact,kLongitmax,kZd,kAz,kTotVar};    // enum variables for the
     69                                                                // matrix columns mapping
    6570
    6671// --------------------------------------------------------------------------
    … …  
    7580    fTitle = title ? title : "Histograms for Disp";
    7681
    77     fSelectedPos = 1;  // default MHDisp flag (selects Correct Disp srcpos)
     82    fSelectedPos = 1;  // default MHDisp flag (selects Correct Disp source position solution)
    7883
    7984    fMatrix = NULL;
     85
     86    // initialize mapping array dimension to the number of columns of fMatrix
     87    fMap.Set(kTotVar);
    8088
    8189    //--------------------------------------------------
    … …  
    110118    fSkymapXY->SetYTitle("Y Disp source position [deg]");
    111119
    112     fHistChi2   = new TH1F("fHistChi2"  ,
     120    fHistMinPar   = new TH1F("fHistMinPar"  ,
    113121         "Distance^2 between Disp and real srcpos", 100, 0., 2.);
    114     fHistChi2->SetDirectory(NULL);
    115     fHistChi2->UseCurrentStyle();
    116     fHistChi2->SetXTitle("Chi^2 = distance^2 Disp to real srcpos [deg^2]");
    117     fHistChi2->SetYTitle("# events");
     122    fHistMinPar->SetDirectory(NULL);
     123    fHistMinPar->UseCurrentStyle();
     124    fHistMinPar->SetXTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]");
     125    fHistMinPar->SetYTitle("# events");
    118126   
    119127    fHistDuDv   = new TH2F("fHistDuDv",
    … …  
    125133    fHistDuDv->SetYTitle("Du = longitudinal distance [deg]");
    126134
    127     fHistChi2Energy   = new TH2F("fHistChi2Energy",
    128          "Chi^2 vs. Energy", 50, 1., 3., 100, 0., 2.);
    129     fHistChi2Energy->SetDirectory(NULL);
    130     fHistChi2Energy->UseCurrentStyle();
    131     fHistChi2Energy->SetXTitle("log10(Energy) [GeV]");
    132     fHistChi2Energy->SetYTitle("Chi^2 = distance^2 Disp to real srcpos [deg^2]");
    133 
    134     fHistChi2Size   = new TH2F("fHistChi2Size",
    135          "Chi^2 vs. Size", 50, 2., 4., 100, 0., 2.);
    136     fHistChi2Size->SetDirectory(NULL);
    137     fHistChi2Size->UseCurrentStyle();
    138     fHistChi2Size->SetXTitle("log10(Size) [#phot]");
    139     fHistChi2Size->SetYTitle("Chi^2 = distance^2 Disp to real srcpos [deg^2]");
     135    fHistMinParEnergy   = new TH2F("fHistMinParEnergy",
     136         "Minimization parameter vs. Energy", 50, 1., 3., 100, 0., 2.);
     137    fHistMinParEnergy->SetDirectory(NULL);
     138    fHistMinParEnergy->UseCurrentStyle();
     139    fHistMinParEnergy->SetXTitle("log10(Energy) [GeV]");
     140    fHistMinParEnergy->SetYTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]");
     141
     142    fHistMinParSize   = new TH2F("fHistMinParSize",
     143         "Minimization parameter vs. Size", 50, 2., 4., 100, 0., 2.);
     144    fHistMinParSize->SetDirectory(NULL);
     145    fHistMinParSize->UseCurrentStyle();
     146    fHistMinParSize->SetXTitle("log10(Size) [#phot]");
     147    fHistMinParSize->SetYTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]");
    140148
    141149    fHistDuEnergy   = new TH2F("fHistDuEnergy",
    … …  
    170178         "Fraction events srcpos well assign vs. Size", 50, 2., 4., 0., 1.);
    171179    fEvCorrAssign->SetDirectory(NULL);
    172     fEvCorrAssign->UseCurrentStyle();
     180    fEvCorrAssign->SetStats(0);
    173181    fEvCorrAssign->SetXTitle("log10(Size) [#phot]");
    174182    fEvCorrAssign->SetYTitle("Fraction of events with srcpos well assign");
    … …  
    186194  delete fHistcosZA;
    187195  delete fSkymapXY;
    188   delete fHistChi2;
     196  delete fHistMinPar;
    189197  delete fHistDuDv;
    190   delete fHistChi2Energy;
    191   delete fHistChi2Size;
     198  delete fHistMinParEnergy;
     199  delete fHistMinParSize;
    192200  delete fHistDuEnergy;
    193201  delete fHistDuSize;
    … …  
    204212Bool_t MHDisp::SetupFill(const MParList *pList)
    205213{
    206     // reset all histograms and Chi^2 computing variables
     214    // reset all histograms and Minimization parameter computing variables
    207215    // before each new eventloop
    208216    fNumEv = 0;
    209     fSumChi2  = 0.;
    210     fChi2  = 0.;
     217    fSumMinPar  = 0.;
     218    fMinPar = (MParameterD*)const_cast<MParList*>(pList)->FindCreateObj("MParameterD", "MinimizationParameter");
     219    if (!fMinPar)
     220    {
     221      *fLog << err << "MParameterD (MinimizationParameter) not found and could not be created... aborting."
     222            << endl;
     223      return kFALSE;
     224    }
     225    fMinPar->SetVal(0);
    211226
    212227    fHistEnergy->Reset();
    … …  
    214229    fHistcosZA->Reset();
    215230    fSkymapXY->Reset();
    216     fHistChi2->Reset();
     231    fHistMinPar->Reset();
    217232    fHistDuDv->Reset();
    218     fHistChi2Energy->Reset();
    219     fHistChi2Size->Reset();
     233    fHistMinParEnergy->Reset();
     234    fHistMinParSize->Reset();
    220235    fHistDuEnergy->Reset();
    221236    fHistDuSize->Reset();
    … …  
    250265    // if fMatrix exists, skip preprocessing and just read events from matrix;
    251266    // if doesn't exist, read variables values from containers, so look for them
     267    cout << "MHDisp::SetupFill:  fMatrix = " << fMatrix << endl;
     268
    252269    if (fMatrix)
    253270      return kTRUE;
    … …  
    291308    if (!fPointing)
    292309    {
    293         *fLog << err << "MPointingPos not found... "
    294               << "MMcEvt is going to be used to get Theta and Phi."
    295               << endl;
    296         //        return kFALSE;
     310        *fLog << err << "MPointingPos not found... aborting." << endl;
     311        return kFALSE;
    297312    }
    298313
    … …  
    306321//
    307322// Set which selection algorithm for the Disp estimated source position
    308 // we want to follow when filling the histograms:
     323// solutions we want to follow when filling the histograms:
    309324//  * iflag == 1 --> Correct source position, the closest to the real srcpos
    310325//                   (only applicable to MC or well known source real data)
    … …  
    323338//
    324339// Sets the Matrix and the array of mapped values for each Matrix column
    325 // (see MDispCalc)
    326340//
    327341void MHDisp::SetMatrixMap(MHMatrix *matrix, TArrayI &map)
    … …  
    334348// --------------------------------------------------------------------------
    335349//
     350// Returns the Matrix and the mapped value for each Matrix column
     351//
     352void MHDisp::GetMatrixMap(MHMatrix* &matrix, TArrayI &map)
     353{
     354    map = fMap;
     355    matrix = fMatrix;
     356}
     357
     358
     359// --------------------------------------------------------------------------
     360//
     361// You can use this function if you want to use a MHMatrix instead of the
     362// given containers for the Disp optimization. This function adds all
     363// necessary columns to the given matrix. Afterwards, you should fill
     364// the matrix with the corresponding data (eg from a file by using
     365// MHMatrix::Fill). Then, if you loop through the matrix (eg using
     366// MMatrixLoop), MHDisp::Fill will take the values from the matrix
     367// instead of the containers.
     368//
     369void MHDisp::InitMapping(MHMatrix *mat)
     370{
     371    if (fMatrix)
     372      return;
     373
     374    fMatrix = mat;
     375
     376    fMap[kX]          = fMatrix->AddColumn("MSrcPosCam.fX");
     377    fMap[kY]          = fMatrix->AddColumn("MSrcPosCam.fY");
     378    fMap[kMeanX]      = fMatrix->AddColumn("MHillas.fMeanX");
     379    fMap[kMeanY]      = fMatrix->AddColumn("MHillas.fMeanY");
     380    fMap[kDelta]      = fMatrix->AddColumn("MHillas.fDelta");
     381
     382    fMap[kSize]       = fMatrix->AddColumn("MHillas.fSize");
     383   
     384    fMap[kM3Long]     = fMatrix->AddColumn("MHillasExt.fM3Long");
     385    fMap[kAsym]       = fMatrix->AddColumn("MHillasExt.fAsym");
     386   
     387    fMap[kEnergy]     = fMatrix->AddColumn("MMcEvt.fEnergy");
     388    fMap[kImpact]     = fMatrix->AddColumn("MMcEvt.fImpact");
     389    fMap[kLongitmax]  = fMatrix->AddColumn("MMcEvt.fLongitmax");
     390   
     391    fMap[kZd]         = fMatrix->AddColumn("MPointingPos.fZd");
     392    fMap[kAz]         = fMatrix->AddColumn("MPointingPos.fAz");
     393}
     394
     395
     396// --------------------------------------------------------------------------
     397//
    336398// Returns a mapped value from the Matrix
    337399//
    … …  
    349411// --------------------------------------------------------------------------
    350412//
    351 // Fill the histograms and sum up the Chi^2 outcome of each processed event
     413// Fill the histograms and sum up the Minimization paramter outcome
     414// of each processed event
    352415//
    353416Bool_t MHDisp::Fill(const MParContainer *par, const Stat_t w)
    … …  
    356419    Double_t impact = 0.;
    357420    Double_t xmax = 0.;
    358     Double_t theta = 0.;
    359     Double_t phi = 0.;
    360421
    361422    if ( fMatrix || (!fMatrix && fMcEvt) )
    362423    { 
    363       energy   = fMatrix ? GetVal(13)  : fMcEvt->GetEnergy();
    364       impact   = fMatrix ? GetVal(14)  : fMcEvt->GetImpact();
    365       xmax     = fMatrix ? GetVal(15)  : fMcEvt->GetLongitmax();
    366      
    367       if (fPointing)
    368       {
    369         theta  = fMatrix ? GetVal(16)  : fPointing->GetZd();
    370         phi    = fMatrix ? GetVal(17)  : fPointing->GetAz();
    371       }
    372       else
    373       {
    374         theta  = fMatrix ? GetVal(16)  : fMcEvt->GetTelescopeTheta()*kRad2Deg;
    375         phi    = fMatrix ? GetVal(17)  : fMcEvt->GetTelescopePhi()*kRad2Deg;
    376       }
     424      energy   = fMatrix ? GetVal(kEnergy)     : fMcEvt->GetEnergy();
     425      impact   = fMatrix ? GetVal(kImpact)     : fMcEvt->GetImpact();
     426      xmax     = fMatrix ? GetVal(kLongitmax)  : fMcEvt->GetLongitmax();     
    377427    }
    378428    else
    379429      *fLog << "No MMcEvt container available (no Energy,ImpactPar or Xmax)" << endl;
    380430
    381     Double_t xsrcpos0   = fMatrix ? GetVal(0)   : fSrcPos->GetX();
    382     Double_t ysrcpos0   = fMatrix ? GetVal(1)   : fSrcPos->GetY();
    383     Double_t xmean0     = fMatrix ? GetVal(2)   : fHil->GetMeanX(); 
    384     Double_t ymean0     = fMatrix ? GetVal(3)   : fHil->GetMeanY();
    385     Double_t delta      = fMatrix ? GetVal(4)   : fHil->GetDelta();
    386    
    387     Double_t size       = fMatrix ? GetVal(5)   : fHil->GetSize();
    388     //    Double_t width0     = fMatrix ? GetVal(6)   : fHil->GetWidth();
    389     //    Double_t length0    = fMatrix ? GetVal(7)   : fHil->GetLength();
    390    
    391     //    Double_t conc       = fMatrix ? GetVal(8)   : fNewPar->GetConc();
    392     //    Double_t leakage1   = fMatrix ? GetVal(9)   : fNewPar->GetLeakage1();
    393     //    Double_t leakage2   = fMatrix ? GetVal(10)  : fNewPar->GetLeakage2();
    394    
    395     Double_t m3long     = fMatrix ? GetVal(11)  : fHilExt->GetM3Long();
    396     Double_t asym       = fMatrix ? GetVal(12)  : fHilExt->GetAsym();
     431    Double_t theta      = fMatrix ? GetVal(kZd)  : fPointing->GetZd();
     432    //    Double_t phi        = fMatrix ? GetVal(kAz)  : fPointing->GetAz();
     433
     434    Double_t xsrcpos0   = fMatrix ? GetVal(kX)       : fSrcPos->GetX();
     435    Double_t ysrcpos0   = fMatrix ? GetVal(kY)       : fSrcPos->GetY();
     436    Double_t xmean0     = fMatrix ? GetVal(kMeanX)   : fHil->GetMeanX(); 
     437    Double_t ymean0     = fMatrix ? GetVal(kMeanY)   : fHil->GetMeanY();
     438    Double_t delta      = fMatrix ? GetVal(kDelta)   : fHil->GetDelta();
     439   
     440    Double_t size       = fMatrix ? GetVal(kSize)    : fHil->GetSize();
     441   
     442    Double_t m3long     = fMatrix ? GetVal(kM3Long)  : fHilExt->GetM3Long();
     443    Double_t asym       = fMatrix ? GetVal(kAsym)    : fHilExt->GetAsym();
    397444
    398445    //------------------------------------------
    399446    // convert to deg
    400     //    Double_t width    = width0  * fMm2Deg;
    401     //    Double_t length   = length0 * fMm2Deg;
    402447    Double_t xsrcpos  = xsrcpos0 * fMm2Deg;
    403448    Double_t ysrcpos  = ysrcpos0 * fMm2Deg;
    … …  
    453498    // Distance between estimated Disp and real source position
    454499    Double_t d2 = (xdisp-xsrcpos)*(xdisp-xsrcpos) +  (ydisp-ysrcpos)*(ydisp-ysrcpos);
    455     fHistChi2->Fill(d2);
     500    fHistMinPar->Fill(d2);
    456501   
    457502    // Longitudinal and transversal distances between Disp and real source positon
    … …  
    465510    fHistcosZA->Fill(cos(theta/kRad2Deg));
    466511   
    467     // to check the size dependence of the optimization
    468     fHistChi2Energy->Fill(log10(energy),d2);
    469     fHistChi2Size->Fill(log10(size),d2);
     512    // to check the size and energy dependence of the optimization
     513    fHistMinParEnergy->Fill(log10(energy),d2);
     514    fHistMinParSize->Fill(log10(size),d2);
    470515    fHistDuEnergy->Fill(log10(energy),Du);
    471516    fHistDuSize->Fill(log10(size),Du);
    … …  
    473518    fHistDvSize->Fill(log10(size),Dv);
    474519   
    475     // variables for Chi^2 calculation (= d^2)
     520    // variables for the Minimization parameter calculation (= d^2)
    476521    fNumEv += 1;
    477     fSumChi2 += d2; 
     522    fSumMinPar += d2; 
    478523   
    479524    return kTRUE;
    … …  
    483528// --------------------------------------------------------------------------
    484529//
    485 // Calculates the final Chi^2 of the Disp optimization
     530// Calculates the final Minimization parameter of the Disp optimization
    486531//
    487532Bool_t MHDisp::Finalize()
    488533{
    489     fChi2 = fSumChi2/fNumEv;
     534    fMinPar->SetVal(fSumMinPar/fNumEv);
    490535    *fLog << endl;
    491     *fLog << "MHDisp::Finalize: SumChi2, NumEv = " << fSumChi2 << ", " << fNumEv << endl;
    492     *fLog << "MHDisp::Finalize: Final Chi2 = " << fChi2 << endl;
    493    
    494     fChi2 = fHistChi2->GetMean();
    495     *fLog << "MHDisp::Finalize: Final Chi2 = " << fChi2 << endl;
     536    *fLog << "MHDisp::Finalize: SumMinPar, NumEv = " << fSumMinPar << ", " << fNumEv << endl;
     537    *fLog << "MHDisp::Finalize: Final MinPar = " << fMinPar->GetVal() << endl;
     538   
     539    fMinPar->SetVal(fHistMinPar->GetMean());
     540    *fLog << "MHDisp::Finalize: Final MinPar = " << fMinPar->GetVal() << endl;
    496541
    497542    SetReadyToSave();
    … …  
    542587    pad->cd(4);
    543588    gPad->SetBorderMode(0);
    544     fHistChi2->SetTitleOffset(1.2,"Y");
    545     fHistChi2->DrawClone(opt);
    546     fHistChi2->SetBit(kCanDelete);
    547     gPad->Modified();
    548 
    549     TProfile *profChi2Energy = fHistChi2Energy->ProfileX("Chi^2 vs. Energy",0,99999,"s");
    550     profChi2Energy->SetXTitle("log10(Energy) [GeV]");
    551     profChi2Energy->SetYTitle("Chi^2 = distance^2 Disp to real srcpos [deg^2]");
     589    fHistMinPar->SetTitleOffset(1.2,"Y");
     590    fHistMinPar->DrawClone(opt);
     591    fHistMinPar->SetBit(kCanDelete);
     592    gPad->Modified();
     593
     594    TProfile *profMinParEnergy = fHistMinParEnergy->ProfileX("Minimization parameter vs. Energy",0,99999,"s");
     595    profMinParEnergy->SetXTitle("log10(Energy) [GeV]");
     596    profMinParEnergy->SetYTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]");
    552597    pad->cd(5);
    553598    gPad->SetBorderMode(0);
    554     profChi2Energy->SetTitleOffset(1.2,"Y");
    555     profChi2Energy->DrawClone(opt);
    556     profChi2Energy->SetBit(kCanDelete);
    557     gPad->Modified();
    558 
    559     TProfile *profChi2Size = fHistChi2Size->ProfileX("Chi^2 vs. Size",0,99999,"s");
    560     profChi2Size->SetXTitle("log10(Size) [#phot]");
    561     profChi2Size->SetYTitle("Chi^2 = distance^2 Disp to real srcpos [deg^2]");
     599    profMinParEnergy->SetTitleOffset(1.2,"Y");
     600    profMinParEnergy->SetStats(0);
     601    profMinParEnergy->DrawClone(opt);
     602    profMinParEnergy->SetBit(kCanDelete);
     603    gPad->Modified();
     604
     605    TProfile *profMinParSize = fHistMinParSize->ProfileX("Minimization parameter vs. Size",0,99999,"s");
     606    profMinParSize->SetXTitle("log10(Size) [#phot]");
     607    profMinParSize->SetYTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]");
    562608    pad->cd(6);
    563609    gPad->SetBorderMode(0);
    564     profChi2Size->SetTitleOffset(1.2,"Y");
    565     profChi2Size->DrawClone(opt);
    566     profChi2Size->SetBit(kCanDelete);
     610    profMinParSize->SetTitleOffset(1.2,"Y");
     611    profMinParSize->SetStats(0);
     612    profMinParSize->DrawClone(opt);
     613    profMinParSize->SetBit(kCanDelete);
    567614    gPad->Modified();
    568615
    … …  
    604651    gPad->SetBorderMode(0);
    605652    profDuEnergy->SetTitleOffset(1.2,"Y");
     653    profDuEnergy->SetStats(0);
    606654    profDuEnergy->DrawClone(opt);
    607655    profDuEnergy->SetBit(kCanDelete);
    … …  
    614662    gPad->SetBorderMode(0);
    615663    profDuSize->SetTitleOffset(1.2,"Y");
     664    profDuSize->SetStats(0);
    616665    profDuSize->DrawClone(opt);
    617666    profDuSize->SetBit(kCanDelete);
    … …  
    634683    gPad->SetBorderMode(0);
    635684    profDvEnergy->SetTitleOffset(1.2,"Y");
     685    profDvEnergy->SetStats(0);
    636686    profDvEnergy->DrawClone(opt);
    637687    profDvEnergy->SetBit(kCanDelete);
    … …  
    644694    gPad->SetBorderMode(0);
    645695    profDvSize->SetTitleOffset(1.2,"Y");
     696    profDvSize->SetStats(0);
    646697    profDvSize->DrawClone(opt);
    647698    profDvSize->SetBit(kCanDelete);
Note: See TracChangeset for help on using the changeset viewer.