source: trunk/MagicSoft/Mars/mtemp/mifae/macros/distancebetweenstars.C@ 5138

 /* ======================================================================== *\
!
! *
! * 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, 05/2004 <mailto:jlopez@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

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, <return> to go on: ");
        timer.TurnOff();

        if (input=="q\n")
            return kFALSE;

        if (input=="\n")
            return kTRUE;
    };

    return kFALSE;
}

Double_t fitfunc(Double_t *x, Double_t *par);

void distancebetweenstars(const TString filename="dc_2004_03_17_01_16_51_20440_Mrk421.root", const TString directory="/nfs/magic/CaCodata/rootdata/Mrk421/Period015/2004_03_17/", const UInt_t numEvents = 100000000)
{

    // general settings
    gROOT->Reset();
    gStyle->SetCanvasColor(0);
    gStyle->SetCanvasBorderMode(0);
    gStyle->SetPadBorderMode(0);
    gStyle->SetFrameBorderMode(0);
    gStyle->SetOptTitle(0);
    gStyle->SetTitleOffset(1.7,"y");
    gStyle->SetPadLeftMargin(0.15);
    gStyle->SetOptStat(111110);
    gStyle->SetOptFit(1);
    gStyle->SetStatColor(0);
    gStyle->SetStatBorderSize(1);
    gStyle->SetStatW(0.2);
    gStyle->SetStatH(0.1);
    gStyle->SetStatX(0.9);
    gStyle->SetStatY(0.9);

  Int_t nbins  = 100;
  Float_t mindist = 0.0;
  Float_t maxdist = 400.0;
  
  TH1F* histStarsDintances = new TH1F("StarsDintances","Distance between stars",nbins,mindist,maxdist);
  histStarsDintances->SetXTitle("Distance [mm]");
  histStarsDintances->SetYTitle("Counts [#]");
  TH1F* histStarsDintances1 = new TH1F("StarsDintances1","Distance between stars",nbins,mindist,maxdist);
  TH1F* histStarsDintances2 = new TH1F("StarsDintances2","Distance between stars [1]-[3]",nbins,mindist,maxdist);
  TH1F* histStarsDintances3 = new TH1F("StarsDintances3","Distance between stars [2]-[3]",nbins,mindist,maxdist);
  //
  // 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);

  MGeomCamMagic geomcam;
  MCameraDC     dccam;
  MStarLocalCam starcam;

  plist.AddToList(&geomcam);
  plist.AddToList(&dccam);
  plist.AddToList(&starcam);

  //
  // Now setup the tasks and tasklist:
  // ---------------------------------
  //

  // Reads the trees of the root file and the analysed branches
  MReadReports read;
  read.AddTree("Currents"); 
  read.AddFile(directory+filename);     // after the reading of the trees!!!
  read.AddToBranchList("MReportCurrents.*");

  MGeomApply geomapl;
  TString continuoslightfile = 
    //    "/home/Javi/mnt_magic_data/CaCo/rootdata/Miscellaneous/Period016/2004_04_16/dc_2004_04_16_04_46_18_22368_Off3c279-2CL100.root";
        "/nfs/magic/CaCodata/rootdata/Miscellaneous/Period016/2004_04_16/dc_2004_04_16_04_46_18_22368_Off3c279-2CL100.root";

  Float_t mindc = 0.7;
  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 = 4.0;
  UInt_t integratedevents = 1;

  MFindStars findstars;
  findstars.SetBlindPixels(blindpixels);
  findstars.SetRingInterest(ringinterest);
  findstars.SetDCTailCut(tailcut);
  findstars.SetNumIntegratedEvents(integratedevents);
  findstars.SetMinuitPrintOutLevel(-1);

  tlist.AddToList(&geomapl);
  tlist.AddToList(&read);
  tlist.AddToList(&dccal);
  tlist.AddToList(&findstars, "Currents");

  //
  // Create and setup the eventloop
  //
  MEvtLoop evtloop;
  evtloop.SetParList(&plist);
  
//   MProgressBar bar;
//   evtloop.SetProgressBar(&bar);
  
  //
  // Execute your analysis
  //
  
//   if (!evtloop.Eventloop(numEvents))
//     return;

  if (!evtloop.PreProcess())
    return;

  Float_t maxchindof = 4.;
  
  while (tlist.Process())
    {
  
      Int_t numStars = starcam.GetNumStars();
      if ( numStars == 3)
        {
          
          for (Int_t first=0; first<numStars; first++)
            {
              if (starcam[first].GetChiSquareNdof()>0. && starcam[first].GetChiSquareNdof()<maxchindof)
                {
                  for (Int_t second=first+1; second<numStars; second++)
                    {
                      if (starcam[second].GetChiSquareNdof()>0. && starcam[second].GetChiSquareNdof()<maxchindof)
                        {
                          Float_t dist = TMath::Sqrt((starcam[first].GetMeanX()-starcam[second].GetMeanX())*
                                                     (starcam[first].GetMeanX()-starcam[second].GetMeanX()) +
                                                     (starcam[first].GetMeanY()-starcam[second].GetMeanY())*
                                                     (starcam[first].GetMeanY()-starcam[second].GetMeanY()));

                            histStarsDintances->Fill(dist);
                          if (first == 0 && second == 1)
                            histStarsDintances1->Fill(dist);
                          else if (first == 0 && second == 2)
                            histStarsDintances2->Fill(dist);
                          else if (first == 1 && second == 2)
                            histStarsDintances3->Fill(dist);

                        }
                    }
                }
            }
        }
    }
  
  

  evtloop.PostProcess();

  tlist.PrintStatistics();

  //Draw results

  //   histStarsDintances->Draw();

// Creates a Root function based on function fitf above
   TF1 *func = new TF1("fitfunc",fitfunc,mindist,maxdist,7);

// Sets initial values and parameter names
   func->SetParNames("ConvF","Max0","Sig0","Max1","Sig1","Max2","Sig2");
   func->SetParameters(300.,500.,5.,500.,5.,500.,5.);

// Fit histogram in range defined by function
   histStarsDintances->Fit("fitfunc","R");

//   histStarsDintances1->Draw();
//   histStarsDintances2->Fit("gaus","0");
//   histStarsDintances2->Draw("same");
//   histStarsDintances2->Fit("gaus","0");
//   histStarsDintances3->Draw("same");
//   histStarsDintances3->Fit("gaus","0");

  if (!HandleInput())

  delete histStarsDintances;
  delete histStarsDintances1;
  delete histStarsDintances2;
  delete histStarsDintances3;
  
}

Double_t fitfunc(Double_t *x, Double_t *par)
{
  Double_t dist[3] = 
    {
      0.64,0.750,1.203
    };
  
  
   Double_t fitval = 
     par[1]*TMath::Exp(-0.5*(x[0]-par[0]*dist[0])*(x[0]-par[0]*dist[0])/(par[2]*par[2])) +
     par[3]*TMath::Exp(-0.5*(x[0]-par[0]*dist[1])*(x[0]-par[0]*dist[1])/(par[4]*par[4])) +
     par[5]*TMath::Exp(-0.5*(x[0]-par[0]*dist[2])*(x[0]-par[0]*dist[2])/(par[6]*par[6]));

   //   cout << "x "  << x[0] << " fitval " << fitval << endl;
   

   return fitval;
}