Index: trunk/MagicSoft/Mars/mtemp/Changelog
===================================================================
--- trunk/MagicSoft/Mars/mtemp/Changelog	(revision 4434)
+++ trunk/MagicSoft/Mars/mtemp/Changelog	(revision 4435)
@@ -18,4 +18,11 @@
 
                                                  -*-*- END OF LINE -*-*-
+
+ 2004/07/30: Nicola Galante
+
+   * mtemp/mpisa/macros/AnalCurrents.C
+      - Added new macros AnalCurrents.C which reads DC files
+        and looks for DC clusters. Still preliminary! 
+        It currently just looks for starting pixels of clusters.
 
  2004/07/22: Sabrina Stark
Index: trunk/MagicSoft/Mars/mtemp/mpisa/macros/AnalCurrents.C
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mpisa/macros/AnalCurrents.C	(revision 4435)
+++ trunk/MagicSoft/Mars/mtemp/mpisa/macros/AnalCurrents.C	(revision 4435)
@@ -0,0 +1,224 @@
+/* ======================================================================== *\
+!
+! *
+! * 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;
+  Int_t count = 0;
+  Double_t val;
+  //
+  // Execute your analysis
+  //
+  //if (!evtloop.Eventloop(numEvents))
+  //  return kFALSE;
+  if (!evtloop.PreProcess())
+    return;
+    
+  while (tlist.Process())
+    {
+      //cout << "Evt #" << ++count << endl;
+      for(Int_t i=0; i<gsNpix; i++){
+	dcevt.AddPixel(i,dccam.GetCurrent(i)/geomcam.GetPixRatio(i),0);
+	dcevt.GetPixelContent(val,i,geomcam,0);
+	//cout << i << "-" << dccam.GetCurrent(i) << "-" << val << endl;
+      }
+      //cout << endl;
+    }
+  tlist.PrintStatistics();
+
+  cout << "Calling test_findstar function" << endl;
+  test_findstar(&dcevt);
+}
+
+
+// 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;
+
+}
+*/
+