source: branches/MarsMoreSimulationTruth/mtemp/mpisa/macros/AnalCurrents.C@ 20115

/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Javier López, 04/2004 <mailto:jlopez@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */
/*
#include <iostream>
#include <stdlib.h>
#include <string.h>
#include <TString.h>
#include <TArrayS.h>

#include "MParList.h"
#include "MTaskList.h"
#include "MGeomCamMagic.h"
#include "MCameraDC.h"

#include "MReadReports.h"
#include "MGeomApply.h"
#include "MEvtLoop.h"
*/
const Int_t gsNpix = 577;
MGeomCamMagic geomcam;

using namespace std;

void AnalCurrents(const TString filename)
{
    //UInt_t numEvents = 1000000;
    UInt_t numEvents = 1;
    
    //
    // Create a empty Parameter List and an empty Task List
    // The tasklist is identified in the eventloop by its name
    //
    MParList  plist;  
    MTaskList tlist;
    plist.AddToList(&tlist);
    
    
    MCameraDC     dccam;
    
    plist.AddToList(&geomcam);
    plist.AddToList(&dccam);
    
    //
    // Now setup the tasks and tasklist:
    // ---------------------------------
    //
    // Reads the trees of the root file and the analysed branches
    MReadReports read;
    read.AddTree("Currents"); 
    read.AddFile(filename);     // after the reading of the trees!!!
    read.AddToBranchList("MReportCurrents.*");
    
    MGeomApply geomapl;
    
    tlist.AddToList(&geomapl);
    tlist.AddToList(&read);

    //
    // Create and setup the eventloop
    //
    MEvtLoop evtloop;
    evtloop.SetParList(&plist);
     
    MCerPhotEvt dcevt;
    
    const Int_t NSteps = 5;  // Number of steps to average DC over events
    Int_t count = 0;         // counter for averaging NSteps events
    Int_t NumAverEvt = 1;    // counter for averaged events stored into dcevt
    Double_t val[gsNpix];    // array for managing the averaging of gsNpix
                             // DC values over NSteps events

    // Initialize to 0 array val
    for(Int_t i=0; i<gsNpix; i++)
        val[i] = 0.0;
    //
    // Execute your analysis
    //
    //if (!evtloop.Eventloop(numEvents))
    //  return kFALSE;
    if (!evtloop.PreProcess())
        return;
    
    while (tlist.Process())
    {
        ++count;
        for(Int_t i=0; i<gsNpix; i++)
        {
            val[i] += dccam.GetCurrent(i);

            if( (count%NSteps) == 0 )
            {
                dcevt.AddPixel(i,val[i]/(Double_t)NSteps,0);
                dcevt.GetPixelContent(val[i],i,geomcam,0);
                val[i] = 0.0;
            }
        }

        // Reset count and call findstar function on this macroevent
        if( (count%NSteps) == 0 )
        {
            count = 0;
            ++NumAverEvt;
            cout << "Calling test_findstar function" << endl;
            test_findstar(&dcevt);
        }
    }

    // Last loop overall pixels to fill with the average 
    // of last count events
    if(count!=0)
    {
        for(Int_t i=0; i<gsNpix; i++)
            dcevt.AddPixel(i,val[i]/(Double_t)count,0);

        cout << "Calling test_findstar function" << endl;
        test_findstar(&dcevt);
    }

    cout << "Number of averaged events stored into container = " << NumAverEvt << endl;
    tlist.PrintStatistics();
    
}


// This macro creates a fake MCerPhotEvt container (with some clusters
// that simulate the starfield seen by the DC currents) and then applies
// the image cleaning and the Hillas algorithm to recognize and classify
// the clusters.
//
void test_findstar(MCerPhotEvt *wDCEvt)
{
    //MGeomCamMagic *geom; 
    MCerPhotEvt *DCEvt = new MCerPhotEvt();
    Int_t cluster[gsNpix];
    Float_t startp[gsNpix];
    
    memset(cluster,-1,gsNpix);
    memset(startp,-1,gsNpix);
    
    TVectorD *startpixs = new TVectorD();
    Double_t tempDC;
    Int_t pixid;
    // fill a fake MCerPhotEvt 
    // loop over all pixels
    for (pixid=0; pixid<gsNpix ; pixid++)
    {
        switch (pixid)
        {
            case 9:
                DCEvt->AddPixel(pixid,1.1,0.0);
                break;
            case 2:
            case 8:
            case 21:
            case 22:
            case 23:
            case 10:
                DCEvt->AddPixel(pixid,1.1,0.0);
                break;
            case 39:
                DCEvt->AddPixel(pixid,1.2,0.0);
                break;
            case 64:
                DCEvt->AddPixel(pixid,1.71,0.0);
                break;
            case 65:
                DCEvt->AddPixel(pixid,1.70,0.0);
                break;
            case 40:
                DCEvt->AddPixel(pixid,1.72,0.0);
                break;
            default:
                DCEvt->AddPixel(pixid);
                break;
        }
    } // end loop over pixels

    if(!(FindStartPixels(wDCEvt, startpixs)))
        cout << "ARGH!@!! In function FindStartPixel(): no starting pixel found" << endl;
    
}

Int_t FindStartPixels(MCerPhotEvt *evt, TVectorD *startpixs)
{
    Double_t currDC;
    Int_t i = 0;
    Double_t DCthr = 10.0;
    
    MCerPhotPix *dcpix;
    MGeomPix *pix;
    
    // look for all the pixels with a DC higher than the DC of neighbour pixels 
    // loop over all pixels
    for (Int_t pixid=0; pixid<gsNpix; pixid++)
    {
        Double_t tempDC;

        pix = (MGeomPix *)geomcam[pixid]; // MGeomCamMagic pixel
        
        // Get the DC value of the current pixel
        evt->GetPixelContent(currDC,pixid,geomcam,0);
        
        Float_t maxDC = 0.0;
        // look for the max DC in the neighbors pixels
        for (Int_t j=0; j < pix->GetNumNeighbors(); j++)
        {
            evt->GetPixelContent(tempDC,pix->GetNeighbor(j),geomcam,0);
            if(tempDC > maxDC) 
                maxDC = tempDC;
        }
        
        // a starting pixel was found: it is added to the array
        if ((currDC>maxDC) && (currDC>DCthr))
        {
            startpixs->ResizeTo(++i);
            startpixs(i-1)=pixid;
            cout << "Starting pixel PixID=" << startpixs(i-1) << " curr=" << currDC << endl;
        }
    }
    
    if(startpixs->GetNrows()>0)
    {
        cout << "Number of starting pixels = " << startpixs->GetNrows() << endl;
        return 1;
    }
    else
        return 0;
}
/*
Int_t FindCluster(Int_t startpix, Int_t cluster, MGeomCam *geom)
{
  return 1;

}
*/