void wobblemap(Float_t maxhadronness = 0.15, Float_t minsize=2000., Float_t maxsize = 1.e10)
{

  //// CHANGE NAME OF FILE TO NOT OVERWRITE FORMER FILES!!! ////////
  TString outname = "Hists_";
  if (maxsize > 1.e6)
      outname += "above";

  outname += Form("%4d", (Int_t)minsize);

  if (maxsize < 1.e6)
    {
      outname += "to";
      outname += Form("%4d", (Int_t)maxsize);
    }

  outname += "_had";

  TString dummy = Form("%3.2f", (Double_t)maxhadronness); // cast is necessary!!
  outname += dummy;
  outname += ".root";

  cout << outname << endl;

  TFile* hfile = new TFile(outname,"recreate");

  MGeomCamMagic mgeom;
  Double_t mm2deg = mgeom.GetConvMm2Deg();

  Double_t CrabDec = 22.014444444; // deg
  Double_t CrabRA  = 5.5755555555; // hours

//   Double_t CrabDec = 19.8197221; // deg
//   Double_t CrabRA  = 8.895194; // hours  // FALSE source for the 3EG0853 data!

  Float_t onregion_radius = 0.25; // degree
  Float_t onr2 = onregion_radius * onregion_radius; // deg2
  Float_t datatime = 0.;

  // CHANGE MC PARAMETERS AND FILE BELOW ACCORDING TO STUDIED DATA SETS!!!!
  // Period 025:
//    TString spotsize = "2.0cm";
//    Float_t minMCtheta = 5.73;   // lower edge of zbin 1
//    Float_t maxMCtheta = 23.79;  // upper edge of zbin 8
//    Float_t n_gamma_produced = 1.e5 * 8. * 6. / 2.; 
  // #Events per file * #zbins * #files per zbin / 2 (test sample only)

  // YOU ALSO HAVE TO CHANGE THE SCALING OF W, L AND THE EXPRESSION FOR DISP!

  // Period 021
  TString spotsize = "1.4cm";
  Float_t minMCtheta = 15.20;   // lower edge of zbin 4
  Float_t maxMCtheta = 28.96;   // upper edge of zbin 12
  Float_t n_gamma_produced = 1.e5 * 9. * 6. / 2.; 
  // #Events per file * #zbins * #files per zbin / 2 (test sample only)



  TChain *MCchain = new TChain("Events");
  //  MCchain->Add("/data2/magic/data/rootdata/Crab/Period025/2005_01_10/RandomForest_1.0_spot/star_S1000_mcgammas.root");

  if (spotsize == "1.4cm")
    MCchain->Add("/data2/magic/data/rootdata/Crab/Period021/2004_09_15/RandomForest_1.4_spot/star_S1000_mcgammas.root");
  else if (spotsize == "2.0cm")
    MCchain->Add("/data2/magic/data/rootdata/Crab/Period025/2005_01_10/RandomForest_2.0_spot/star_S1000_mcgammas.root");

  //////////////////////////////////////////////////////////////////////////////////////////////////

  TChain *Wchain = new TChain("Events");

  Wchain->Add("/data2/magic/data/rootdata/Crab/Period021/2004_09_15/RandomForest_1.4_spot/star_S1000_20040915_CrabNebulaW1.root");
  datatime += 1062.3;

  Wchain->Add("/data2/magic/data/rootdata/Crab/Period021/2004_09_24/RandomForest_1.4_spot/star_S1000_20040924_CrabNebulaW1.root");
  datatime += 1316.2;

  Wchain->Add("/data2/magic/data/rootdata/Crab/Period021/2004_09_15/RandomForest_1.4_spot/star_S1000_20040915_CrabNebulaW2.root");
  datatime += 1241.1;

  Wchain->Add("/data2/magic/data/rootdata/Crab/Period021/2004_09_18/RandomForest_1.4_spot/star_S1000_20040918_CrabNebulaW2.root");
  datatime += 1388.1;

  Wchain->Add("/data2/magic/data/rootdata/Crab/Period021/2004_09_24/RandomForest_1.4_spot/star_S1000_20040924_CrabNebulaW2.root");
  datatime += 767.0;


//   Wchain->Add("/data2/magic/data/rootdata/Crab/Period025/2005_01_10/RandomForest_2.0_spot/star_S1000_20050110_CrabNebulaW1.root");
//   datatime += 9469.7;
  
//   Wchain->Add("/data2/magic/data/rootdata/Crab/Period025/2005_01_11/RandomForest_2.0_spot/star_S1000_20050111_CrabNebulaW1.root");
//   datatime += 4400.1;

//   Wchain->Add("/data2/magic/data/rootdata/Crab/Period025/2005_01_11/RandomForest_2.0_spot/star_S1000_20050111_CrabNebulaW2.root");
//   datatime += 4673.1;


  // November data: WE SKIP THEM, they are very few!
  //   Wchain->Add("/data2/magic/data/rootdata/Crab/Period023/2004_11_23/RandomForest_1.0_spot/star_S1000_20041123_CrabNebulaW1.root");
  //   datatime += 394.2;
  //   Wchain->Add("/data2/magic/data/rootdata/Crab/Period023/2004_11_23/RandomForest_1.0_spot/star_S1000_20041123_CrabNebulaW2.root");
  //   datatime += 0.0;  // No data with theta below 28.96 deg!

   
//   Wchain->Add("/data2/magic/data/rootdata/3EG0853/Period024/2004_12_10/RandomForest_2.0_spot/star_S1000_20041210_3EG0853.root");
//   datatime += 4487.9;


 
  TH2F *CamMap = new TH2F("CamMap", "Events incoming direction map ",
			    50, -2.5, 2.5,50, -2.5, 2.5);


  TH2F *SkyMap = new TH2F("SkyMap", "Events incoming direction map ",
			  70, CrabRA-0.175, CrabRA+0.175, 
			  70, CrabDec-2.8315, CrabDec+2.8315);

  TH2F *ONandOFF = new TH2F("ONandOFF", "ON and OFF regions",
			    140, CrabRA-0.175, CrabRA+0.175, 
			    140, CrabDec-2.8315, CrabDec+2.8315);

  ////////////// BACKGROUND MODEL: //////////////////////////////////////////////////

  // Radial distribution, dN/dr with trespect to camera center:
  // SIZE > 2000
  TF1* myfunc = new TF1("myfunc","x*([0]+[1]*log(x)+[2]*log(x)*log(x))", 0., 2.);
  myfunc->SetParameters(2.06765e+02, -2.36253e+02, -2.13230e+01);

  // SIZE 1000-2000
  // TF1* myfunc = new TF1("myfunc","x*([0]+[1]*log(x)+[2]*log(x)*log(x)+[3]*log(x)**4)", 0., 2.);
  // myfunc->SetParameters(6.45011e+02, -1.01606e+03, 1.63394e+02, -3.60749e+01);

  TH1F *rndmrad = new TH1F("rndmrad","",200, 0., 2.);
  rndmrad->FillRandom("myfunc",1e6);
  rndmrad->Draw("e");

  TH2F *BgMap = (TH2F*) SkyMap->Clone("BgMap");
  TRandom dice;

  ///////////////////////////////////////////////////////////////////////////////////


  TH1F *Alpha = new TH1F("Alpha", "", 60, -90., 90.);
  TH1F *AbsAlpha = new TH1F("AbsAlpha", "", 45, 0., 90.);

  TH1F *WidthOn = new TH1F("WidthOn", "WIDTH", 120, 0., 0.60);
  TH1F *WidthOff  = (TH1F*) WidthOn->Clone("WidthOff");
  WidthOn->Sumw2();
  WidthOff->Sumw2();

  TH1F *LengthOn = new TH1F("LengthOn", "LENGTH", 50, 0., 1.);
  TH1F *LengthOff  = (TH1F*) LengthOn->Clone("LengthOff");
  LengthOn->Sumw2();
  LengthOff->Sumw2();

  TH1F *SWidthOn = new TH1F("SWidthOn", "SWIDTH", 75, 0.3, 5.3);
  TH1F *SWidthOff  = (TH1F*) SWidthOn->Clone("SWidthOff");
  SWidthOn->Sumw2();
  SWidthOff->Sumw2();

  TH1F *SLengthOn = new TH1F("SLengthOn", "SLENGTH", 45, 0.2, 3.2);
  TH1F *SLengthOff  = (TH1F*) SLengthOn->Clone("SLengthOff");
  SLengthOn->Sumw2();
  SLengthOff->Sumw2();

  TH1F *ConcOn = new TH1F("ConcOn", "CONC", 30, 0., 0.6);
  TH1F *ConcOff  = (TH1F*) ConcOn->Clone("ConcOff");
  ConcOn->Sumw2();
  ConcOff->Sumw2();

  TH1F *M3LongOn = new TH1F("M3LongOn", "M3LONG", 80, -0.8, 0.8);
  TH1F *M3LongOff  = (TH1F*) M3LongOn->Clone("M3LongOff");
  M3LongOn->Sumw2();
  M3LongOff->Sumw2();

  TH1F *DistOn = new TH1F("DistOn", "DIST", 40, 0., 2.);
  TH1F *DistOff  = (TH1F*) DistOn->Clone("DistOff");
  DistOn->Sumw2();
  DistOff->Sumw2();

  TH1F *SizeOn = new TH1F("SizeOn", "LOG10SIZE", 60, 2.5, 5.5);
  TH1F *SizeOff  = (TH1F*) SizeOn->Clone("SizeOff");
  SizeOn->Sumw2();
  SizeOff->Sumw2();


  TH1F *HadronnessOn = new TH1F("HadronnessOn", "Hadronness", 101, 0., 1.01);
  TH1F *HadronnessOff = (TH1F*) HadronnessOn->Clone("HadronnessOff");
  HadronnessOn->Sumw2();
  HadronnessOff->Sumw2();

  TH1F* r_source = new TH1F("r_source", "r_source", 50, 0., 2.5);


  // Now the MC histograms:

  TH1F *AlphaMC = (TH1F*) Alpha->Clone("AlphaMC");
  TH1F *AbsAlphaMC = (TH1F*) AbsAlpha->Clone("AbsAlphaMC");
  TH1F *WidthMC = (TH1F*) WidthOn->Clone("WidthMC");
  TH1F *LengthMC  = (TH1F*) LengthOn->Clone("LengthMC");
  TH1F *SWidthMC  = (TH1F*) SWidthOn->Clone("SWidthMC");
  TH1F *SLengthMC  = (TH1F*) SLengthOn->Clone("SLengthMC");
  
  TH1F *ConcMC  = (TH1F*) ConcOn->Clone("ConcMC");
  TH1F *M3LongMC = (TH1F*) M3LongOn->Clone("M3LongMC");
  TH1F *M3LongMC_all = (TH1F*) M3LongOn->Clone("M3LongMC_all");
  TH1F *DistMC  = (TH1F*) DistOn->Clone("DistMC");
  TH1F *SizeMC  = (TH1F*) SizeOn->Clone("SizeMC");


  TH1F *HadronnessMC = (TH1F*) HadronnessOn->Clone("HadronnessMC");

  TH1F *EnergyMC = new TH1F("EnergyMC", "MC, log_{10}E_{\gamma} (GeV)", 1000, 0., 2000.);
  EnergyMC->Sumw2();

  ////////////////////////////////////////////////////////////////

  MObservatory mobserv;
  MStarCamTrans* starcamtrans = new MStarCamTrans(mgeom, mobserv);


  MHillas *mhil = new MHillas();
  MHillasSrc *msrc = new MHillasSrc();
  MHillasExt *mhilext = new MHillasExt();
  MNewImagePar *mnew = new MNewImagePar();
  MPointingPos *mpoint = new MPointingPos();
  MHadronness  *mhad = new MHadronness();
  MPointingPos *mpoint = new MPointingPos();
  MTime        *mtime  = new MTime();

  Wchain->SetBranchStatus("*", 0);
  Wchain->SetBranchStatus("MHillas.fSize", 1);
  Wchain->SetBranchStatus("MHillas.fWidth", 1);
  Wchain->SetBranchStatus("MHillas.fLength", 1);
  Wchain->SetBranchStatus("MHillas.fMeanX", 1);
  Wchain->SetBranchStatus("MHillas.fMeanY", 1);
  Wchain->SetBranchStatus("MHillas.fCosDelta", 1);
  Wchain->SetBranchStatus("MHillas.fSinDelta", 1);
  Wchain->SetBranchStatus("MHillasSrc.fCosDeltaAlpha", 1);
  Wchain->SetBranchStatus("MHillasExt.fM3Long", 1);
  Wchain->SetBranchStatus("MNewImagePar.fNumCorePixels", 1);
  Wchain->SetBranchStatus("MNewImagePar.fConc", 1);
  Wchain->SetBranchStatus("MPointingPos.fZd",1);
  Wchain->SetBranchStatus("MHadronness.fHadronness",1);
  Wchain->SetBranchStatus("MPointingPos.fDec",1);
  Wchain->SetBranchStatus("MPointingPos.fRa",1);
  Wchain->SetBranchStatus("MTime.*",1);
  Wchain->SetBranchStatus("MHillasSrc.fDist",1);

  Wchain->SetBranchAddress("MHillas.",      &mhil);
  Wchain->SetBranchAddress("MHillasSrc.",   &msrc);
  Wchain->SetBranchAddress("MHillasExt.",   &mhilext);
  Wchain->SetBranchAddress("MNewImagePar.", &mnew);
  Wchain->SetBranchAddress("MPointingPos.", &mpoint);
  Wchain->SetBranchAddress("MHadronness.",  &mhad);
  Wchain->SetBranchAddress("MTime.",        &mtime);


  for (Int_t i = 0; i < Wchain->GetEntries(); i++) // event count starts in 0
    {
      if ((i+1)%10000 == 0)
	cout << i+1 << " of " << Wchain->GetEntries() << endl;

      Wchain->GetEvent(i);

      Float_t size = mhil->GetSize();

      if (size < minsize)
	continue;
      else if (size > maxsize)
	continue;

      Float_t log10size = log10(size);

      Float_t ZA = mpoint.GetZd();
      // Exclude events with no drive info. For Crab ZAmin = 6.8
      if (ZA < minMCtheta)
	continue;
      else if (ZA > maxMCtheta)
	continue;   // MC up to ~29 deg loaded, but data may have lower upper end!


      Float_t width  = mhil->GetWidth()* mm2deg;;
      Float_t length = mhil->GetLength()* mm2deg;;


      Float_t disp;

      // spot 1 cm:
      //      Float_t disp   = (1.-(width/length))*(2.68576-0.830937*log10size+0.1227*log10size*log10size);


      if (spotsize == "1.4cm")
	disp   = (1.-(width/length))*(3.64292-1.28713*log10size+0.179722*log10size*log10size);

      else if (spotsize == "2.0cm")
	disp   = (1.-(width/length))*(4.11949-1.39848*log10size+0.183514*log10size*log10size);


      Float_t xmean = mhil->GetMeanX()* mm2deg;;
      Float_t ymean = mhil->GetMeanY()* mm2deg;;

      Float_t cosdelta = mhil->GetCosDelta();
      Float_t sindelta = mhil->GetSinDelta();

      Float_t asylon = mhilext->GetM3Long();


      Double_t sourcex, sourcey;

      // Choose right orientation:

      if (asylon < 0.)
      	{
	  sourcex = xmean+cosdelta*disp;
	  sourcey = ymean+sindelta*disp;
 	}
       else
 	{
 	  sourcex = xmean-cosdelta*disp;
 	  sourcey = ymean-sindelta*disp;
 	}

      Double_t SourceDec = 0;
      Double_t SourceRA  = 0.;
      Double_t SourceHA  = 0.;

      Double_t CenterDec = mpoint->GetDec();
      Double_t CenterRA  = mpoint->GetRa();

      Double_t gmst = mtime->GetGmst(); // Greenwich mean sidereal time [rad]

      gmst *= 24. / (2.*TMath::Pi()); // Conversion to hours

      Double_t Wlongitude = -1. * mobserv.GetLongitudeDeg()/15.; // h
      // Warning: GetLongitude return is positive if location is E of Greenwich!!
      Double_t lst =  gmst - Wlongitude; // Local sidereal time [h]

      Double_t CenterHA = lst - CenterRA;

      if (CenterHA < 0.)
	CenterHA += 24.;
      else if (CenterHA > 24.)
	CenterHA -= 24.;
      // Hour angle of center of camera [h]


      // Find Sky coordinates of reconstructed direction:

      starcamtrans->Cel0CamToCel(CenterDec, CenterHA, sourcex/mm2deg, 
				 sourcey/mm2deg, SourceDec, SourceHA);
      SourceRA = lst - SourceHA; // Estimated source RA [h]
      if (SourceRA < 0.)
	SourceRA += 24.;
      else if (SourceRA > 24.)
	SourceRA -= 24.;


      // Generate background

      Float_t r = rndmrad->GetRandom();
      Float_t phi = 2.*TMath::Pi()*dice.Rndm();
      Float_t xbg = r * cos(phi);
      Float_t ybg = r * sin(phi);
      Double_t BGDec = 0.;
      Double_t BGHA = 0.;
      Double_t BGRA = 0.;
      // Sky coordinates:
      starcamtrans->Cel0CamToCel(CenterDec, CenterHA, xbg/mm2deg, ybg/mm2deg, BGDec, BGHA);
      BGRA = lst - BGHA; // Estimated source RA [h]
      if (BGRA < 0.)
	BGRA += 24.;
      else if (BGRA > 24.)
	BGRA -= 24.;


      // Now calculate Alpha with respect to true Crab position

      Double_t CrabHA = lst - CrabRA;
      if (CrabHA > 24.)
	CrabHA -= 24.;
      else if (CrabHA < 0.)
	CrabHA += 24.;
      
      Double_t xcrab, ycrab; // mm

      starcamtrans->Cel0CelToCam(CenterDec, CenterHA, CrabDec, CrabHA, 
				 xcrab, ycrab);

      xcrab *= mm2deg;
      ycrab *= mm2deg; // Convert to deg

      // Calculate Alpha relative to Crab:

      Float_t scalar = ((xmean-xcrab)*cosdelta + (ymean-ycrab)*sindelta) /
	sqrt( pow(xmean-xcrab, 2.) + pow(ymean-ycrab, 2.) );

      Float_t alpha = 180./TMath::Pi()*acos(scalar);

      if (alpha > 90.)
	alpha -= 180.;


      // Hadronness cut

      Float_t hadronness = mhad->GetHadronness();
      if (hadronness > maxhadronness)
	continue;

      // WARNING!: If scaling is changed, do it also in the MC part !!

//       Float_t swidth  = width/(-0.028235+(0.03231*log10size));
//       Float_t slength = length/(-0.13527+(0.09876*log10size));

// Spot 1 cm:
//      Float_t swidth  = width/(-0.0332763+(0.0320227*log10size));
//      Float_t slength = length/(-0.174697+(0.107548*log10size));

      Float_t swidth, slength;

      if (spotsize == "1.4cm")
	{
	  swidth  = width/(-0.0308984+(0.0327119*log10size));
	  slength = length/(-0.181605+(0.109609*log10size));
	}
      else if (spotsize == "2.0cm")
	{
	  Float_t swidth  = width/(-0.0279008+(0.0344538*log10size));
	  Float_t slength = length/(-0.186394+(0.111055*log10size));
	}

      Float_t conc = mnew->GetConc();
      Float_t dist = msrc->GetDist()*mm2deg;


      Float_t crab_from_center = sqrt(xcrab*xcrab+ycrab*ycrab);
      Float_t source_from_center = sqrt(sourcex*sourcex+sourcey*sourcey);
      // Skip events reconstructed further than 2 deg from center:

      if (source_from_center > 2.)
	continue;

      // Fill histograms

      CamMap->Fill(sourcex, sourcey);
      SkyMap->Fill(SourceRA, SourceDec);

      BgMap->Fill(BGRA, BGDec);


      // Fill histogram of distance from reconstructed source position to camera
      // center, but exclude the quadrant where the candidate source is.... (dirty trick)

      if ( (xcrab*sourcex+ycrab*sourcey)/
	   source_from_center / crab_from_center < 0.707) // cos 45 deg
	r_source->Fill(source_from_center);


      // ON region
      
      Float_t ondist2 = (sourcex-xcrab)*(sourcex-xcrab) + 
	(sourcey-ycrab)*(sourcey-ycrab);

      Float_t offdist2_a = (sourcex+xcrab)*(sourcex+xcrab) + 
	(sourcey+ycrab)*(sourcey+ycrab);
      Float_t offdist2_b = (sourcex-ycrab)*(sourcex-ycrab) + 
	(sourcey+xcrab)*(sourcey+xcrab);
      Float_t offdist2_c = (sourcex+ycrab)*(sourcex+ycrab) + 
	(sourcey-xcrab)*(sourcey-xcrab);


      if (ondist2 < onr2)
	{
	  WidthOn->Fill(width);
	  LengthOn->Fill(length);
	  SWidthOn->Fill(swidth);
	  SLengthOn->Fill(slength);
	  ConcOn->Fill(conc);
	  DistOn->Fill(dist);


	  // 3rd moment with sign referred to source position:

	  //      Float_t m3long = asylon * TMath::Sign(1., msrc->GetCosDeltaAlpha()) * mm2deg;
	  Float_t m3long = asylon * TMath::Sign(1., xmean-xcrab) * mm2deg;

	  M3LongOn->Fill(m3long);
	  SizeOn->Fill(log10size);

	  HadronnessOn->Fill(hadronness);
	  ONandOFF->Fill(SourceRA, SourceDec);
	}

      if (offdist2_a < onr2 || offdist2_b < onr2 || offdist2_c < onr2)
	{
	  WidthOff->Fill(width, 1./3.);
	  LengthOff->Fill(length, 1./3.);
	  SWidthOff->Fill(swidth, 1./3.);
	  SLengthOff->Fill(slength, 1./3.);
	  ConcOff->Fill(conc, 1./3.);
	  SizeOff->Fill(log10size, 1./3.);

	  HadronnessOff->Fill(hadronness, 1./3.);
	  ONandOFF->Fill(SourceRA, SourceDec);

	  if (offdist2_a < onr2)
	    {
	      DistOff->Fill(sqrt((xmean+xcrab)*(xmean+xcrab)+(ymean+ycrab)*(ymean+ycrab)), 1./3.);

	      m3long = asylon * TMath::Sign(1.,xmean+xcrab) * mm2deg;
	      M3LongOff->Fill(m3long, 1./3.);
	    }
	  else if (offdist2_b < onr2)
	    {
	      DistOff->Fill(sqrt((xmean-ycrab)*(xmean-ycrab)+(ymean+xcrab)*(ymean+xcrab)), 1./3.);

	      m3long = asylon * TMath::Sign(1.,xmean-ycrab) * mm2deg;
	      M3LongOff->Fill(m3long, 1./3.);
	    }

	  else if (offdist2_c < onr2)
	    {
	      DistOff->Fill(sqrt((xmean+ycrab)*(xmean+ycrab)+(ymean-xcrab)*(ymean-xcrab)), 1./3.);

	      m3long = asylon * TMath::Sign(1.,xmean+ycrab) * mm2deg;
	      M3LongOff->Fill(m3long, 1./3.);
	    }
	}


      // Very rough cut on 3rd moment related to Crab position:
      // (only for alpha plot)

      if ( (asylon < 0. && xcrab < xmean) ||
	   (asylon > 0. && xcrab > xmean) )
	continue;

      Alpha->Fill(alpha);
      AbsAlpha->Fill(abs(alpha));

    }

  // Convert y-axis units to Hz

  WidthOn->Scale(1./datatime);
  LengthOn->Scale(1./datatime);
  SWidthOn->Scale(1./datatime);
  ConcOn->Scale(1./datatime);
  DistOn->Scale(1./datatime);
  M3LongOn->Scale(1./datatime);
  SizeOn->Scale(1./datatime);
  SLengthOn->Scale(1./datatime);
  HadronnessOn->Scale(1./datatime);

  WidthOff->Scale(1./datatime);
  LengthOff->Scale(1./datatime);
  SWidthOff->Scale(1./datatime);
  ConcOff->Scale(1./datatime);
  DistOff->Scale(1./datatime);
  M3LongOff->Scale(1./datatime);
  SizeOff->Scale(1./datatime);
  SLengthOff->Scale(1./datatime);
  HadronnessOff->Scale(1./datatime);

  TH1F* WidthDiff = (TH1F*) WidthOn->Clone("WidthDiff");
  WidthDiff->Add(WidthOff, -1.);

  TH1F* LengthDiff = (TH1F*) LengthOn->Clone("LengthDiff");
  LengthDiff->Add(LengthOff, -1.);

  TH1F* SWidthDiff = (TH1F*) SWidthOn->Clone("SWidthDiff");
  SWidthDiff->Add(SWidthOff, -1.);

  TH1F* SLengthDiff = (TH1F*) SLengthOn->Clone("SLengthDiff");
  SLengthDiff->Add(SLengthOff, -1.);

  TH1F* ConcDiff = (TH1F*) ConcOn->Clone("ConcDiff");
  ConcDiff->Add(ConcOff, -1.);

  TH1F* DistDiff = (TH1F*) DistOn->Clone("DistDiff");
  DistDiff->Add(DistOff, -1.);

  TH1F* M3LongDiff = (TH1F*) M3LongOn->Clone("M3LongDiff");
  M3LongDiff->Add(M3LongOff, -1.);

  TH1F* SizeDiff = (TH1F*) SizeOn->Clone("SizeDiff");
  SizeDiff->Add(SizeOff, -1.);

  TH1F* HadronnessDiff = (TH1F*) HadronnessOn->Clone("HadronnessDiff");
  HadronnessDiff->Add(HadronnessOff, -1.);


  TH2F* SkyMapRaw = SkyMap->Clone("SkyMapRaw");
  TH2F* SkyMapSubtracted= SkyMap->Clone("SkyMapSubtracted");

  SkyMapSubtracted->Add(BgMap, -1.*(SkyMap->Integral() - 
				    datatime*WidthDiff->Integral()) /
			BgMap->Integral());

  SkyMap->Scale(1./datatime);
  BgMap->Scale(1./datatime);

  SkyMap->Add(BgMap, -1.*(SkyMap->Integral()-WidthDiff->Integral())/BgMap->Integral());


  //  gStyle->SetPalette(8,0);   // Greyscale
  gStyle->SetPalette(1,0);   // Pretty palette
  //  gStyle->SetPalette(51,0);  // Deep sea palette

  CamMap->SetStats(0);
  SkyMapRaw->SetStats(0);
  SkyMapSubtracted->SetStats(0);
  SkyMap->SetStats(0);

  //  CamMap->Draw("zcol");
  //  SkyMap->Draw("zcol");

  //  Alpha->Draw();
  //  AbsAlpha->Draw();

//   WidthOn->Draw("e");
//   WidthOff->Draw("same");

//   LengthOn->Draw("e");
//   LengthOff->Draw("histo,same");


  ///////////////////// MONTE CARLO /////////////////////////

  Double_t xcrab = 120.; //mm
  Double_t ycrab = 0.;

  xcrab *= mm2deg;
  ycrab *= mm2deg; // Convert to deg


  mhil = new MHillas();
  mhilsrc = new MHillasSrc();
  mhilext = new MHillasExt();
  mnew = new MNewImagePar();
  mpoint = new MPointingPos();
  mhad = new MHadronness();
  MMcEvt* mcevt = new MMcEvt();

  MCchain->SetBranchStatus("*", 0);
  MCchain->SetBranchStatus("MHillas.fSize", 1);
  MCchain->SetBranchStatus("MHillas.fWidth", 1);
  MCchain->SetBranchStatus("MHillas.fLength", 1);
  MCchain->SetBranchStatus("MHillas.fMeanX", 1);
  MCchain->SetBranchStatus("MHillas.fMeanY", 1);
  MCchain->SetBranchStatus("MHillas.fCosDelta", 1);
  MCchain->SetBranchStatus("MHillas.fSinDelta", 1);
  MCchain->SetBranchStatus("MHillasSrc.fCosDeltaAlpha", 1);
  MCchain->SetBranchStatus("MHillasSrc.fDist", 1);
  MCchain->SetBranchStatus("MHillasExt.fM3Long", 1);
  MCchain->SetBranchStatus("MNewImagePar.fNumCorePixels", 1);
  MCchain->SetBranchStatus("MNewImagePar.fConc", 1);
  MCchain->SetBranchStatus("MHadronness.fHadronness",1);
  MCchain->SetBranchStatus("MMcEvt.fTelescopeTheta",1);
  MCchain->SetBranchStatus("MMcEvt.fEnergy",1);

  MCchain->SetBranchAddress("MHillas.",      &mhil);
  MCchain->SetBranchAddress("MHillasExt.",   &mhilext);
  MCchain->SetBranchAddress("MNewImagePar.", &mnew);
  MCchain->SetBranchAddress("MHadronness.",  &mhad);
  MCchain->SetBranchAddress("MMcEvt.",       &mcevt);
  MCchain->SetBranchAddress("MHillasSrc.",   &msrc);


  // Original # gammas in the MC sample between 1 and 30 TeV (estimate!):
  Float_t n_gamma_produced_1to30TeV =  
    n_gamma_produced*(pow(1000.,-1.6) - pow(30000.,-1.6)) / 
    (pow(10.,-1.6) - pow(30000.,-1.6));

  // Correction due to the modification of the spectrum: -2.62 instead of -2.6
  n_gamma_produced_1to30TeV *=  (pow(1000.,-1.62) - pow(30000.,-1.62)) /
    (pow(1000.,-1.6) - pow(30000.,-1.6));

  Float_t evts2rate =   1.74e-11 / n_gamma_produced_1to30TeV; // cm-2 s-1
  evts2rate *=  TMath::Pi()*pow(300.e2, 2.); // * MC production  area in cm2 => s-1

  cout << "Processing Monte Carlo gammas...." << endl;

  Int_t MC_gammas_in_off = 0;
  Int_t MC_gammas_in_on  = 0;

  for (Int_t i = 0; i < MCchain->GetEntries(); i++) // event count starts in 0
    {
      if ((i+1)%10000 == 0)
	cout << i+1 << " of " << MCchain->GetEntries() << endl;

      MCchain->GetEvent(i);

      Float_t size = mhil->GetSize();

      if (size < minsize)
	continue;
      else if (size > maxsize)
	continue;

      Float_t log10size = log10(size);

      Float_t ZA = mcevt->GetTelescopeTheta()*180./TMath::Pi();

      // Exclude MC from zbin0. For Crab ZAmin = 6.8 deg. We do not make exactly that cut
      // to simplify the flux calculation.

      if (ZA < minMCtheta)  // Corresponds to the upper theta limit of zbin0
	continue;
      else if (ZA > maxMCtheta)
	continue;

      Float_t width  = mhil->GetWidth()* mm2deg;;
      Float_t length = mhil->GetLength()* mm2deg;;

      //      Float_t disp   = (1.-(width/length))*(0.9776+0.101062*log10(size));

      // spot 1 cm:
      //      Float_t disp   = (1.-(width/length))*(2.68576-0.830937*log10size+0.1227*log10size*log10size);

      Float_t disp;

      if (spotsize == "1.4cm")
	disp   = (1.-(width/length))*(3.64292-1.28713*log10size+0.179722*log10size*log10size);
      else if (spotsize == "2.0cm")
	disp   = (1.-(width/length))*(4.11949-1.39848*log10size+0.183514*log10size*log10size);

      Float_t xmean = mhil->GetMeanX()* mm2deg;;
      Float_t ymean = mhil->GetMeanY()* mm2deg;;

      Float_t cosdelta = mhil->GetCosDelta();
      Float_t sindelta = mhil->GetSinDelta();

      Float_t asylon = mhilext->GetM3Long();

      Double_t sourcex, sourcey;

      // Choose right orientation:

      if (asylon < 0.)
      	{
	  sourcex = xmean+cosdelta*disp;
	  sourcey = ymean+sindelta*disp;
 	}
       else
 	{
 	  sourcex = xmean-cosdelta*disp;
 	  sourcey = ymean-sindelta*disp;
 	}


      // Calculate Alpha relative to Crab:

      Float_t scalar = ((xmean-xcrab)*cosdelta + (ymean-ycrab)*sindelta) /
	sqrt( pow(xmean-xcrab, 2.) + pow(ymean-ycrab, 2.) );

      Float_t alpha = 180./TMath::Pi()*acos(scalar);

      if (alpha > 90.)
	alpha -= 180.;


      // Hadronness cut

      Float_t hadronness = mhad->GetHadronness();
      if (hadronness > maxhadronness)
	continue;

//       Float_t swidth  = width/(-0.028235+(0.03231*log10size));
//       Float_t slength = length/(-0.13527+(0.09876*log10size));

// Spot 1 cm:
//      Float_t swidth  = width/(-0.0332763+(0.0320227*log10size));
//      Float_t slength = length/(-0.174697+(0.107548*log10size));

      Float_t swidth, slength;

      if (spotsize == "1.4cm")
	{
	  swidth  = width/(-0.0308984+(0.0327119*log10size));
	  slength = length/(-0.181605+(0.109609*log10size));
	}
      else if (spotsize == "2.0cm")
	{
	  swidth  = width/(-0.0279008+(0.0344538*log10size));
	  slength = length/(-0.186394+(0.111055*log10size));
	}

      Float_t conc = mnew->GetConc();
      Float_t dist = msrc->GetDist()*mm2deg;

      // ON region
      
      Float_t ondist2 = (sourcex-xcrab)*(sourcex-xcrab) + 
	(sourcey-ycrab)*(sourcey-ycrab);

      Float_t offdist2_a = (sourcex+xcrab)*(sourcex+xcrab) + 
	(sourcey+ycrab)*(sourcey+ycrab);
      Float_t offdist2_b = (sourcex-ycrab)*(sourcex-ycrab) + 
	(sourcey+xcrab)*(sourcey+xcrab);
      Float_t offdist2_c = (sourcex+ycrab)*(sourcex+ycrab) + 
	(sourcey-xcrab)*(sourcey-xcrab);


      Float_t energy = mcevt->GetEnergy(); // GeV

      // Calculate weight to account for true Crab spectrum:
      // We take the (curved) shape of the spectrum BELOW 1 TeV from
      // Astropart.Phys. 19 (2003) 339-350 (Fabrizio, Konopelko):
      // dN/dE = C*(E/1TeV)^[-2.47-0.11*log(E/1TeV)]
      // normalizing it above 1 TeV with the value from ApJ614

      // Integrated Crab flux 1 - 30 TeV from Aharonian et al. 2004, ApJ 614: 
      // 1.74e-11 cm-2 s-1 (shape: power law with index -2.62)

      Float_t mcweight;
      if (energy > 1000.)
	mcweight = pow(energy/1000., 2.6-2.62); // correct spectrum
      else
	{
	  mcweight = pow(energy/1000., (2.6-2.47-0.11*log(energy/1000.)));
	}

      mcweight *= evts2rate;  // weight to get results in Hz

      // 3rd moment with sign referred to source position, convert to degree:

      //      Float_t m3long = asylon * TMath::Sign(1., msrc->GetCosDeltaAlpha()) * mm2deg;

      Float_t m3long = asylon * TMath::Sign(1., xmean-xcrab) * mm2deg;

      M3LongMC_all->Fill(m3long, mcweight);

      if (ondist2 < onr2)
	{
	  WidthMC->Fill(width, mcweight);
	  LengthMC->Fill(length, mcweight);
	  SWidthMC->Fill(swidth, mcweight);
	  SLengthMC->Fill(slength, mcweight);

	  ConcMC->Fill(conc, mcweight);
	  DistMC->Fill(dist, mcweight);
	  M3LongMC->Fill(m3long, mcweight);
	  SizeMC->Fill(log10size, mcweight);

	  HadronnessMC->Fill(hadronness, mcweight);
	  EnergyMC->Fill(energy, mcweight);

	  MC_gammas_in_on++;
	}

      if (offdist2_a < onr2 || offdist2_b < onr2 || offdist2_c < onr2)
	MC_gammas_in_off++;

      // Very rough cut on 3rd moment related to Crab position:
      // (only for alpha plot)

      if ( (asylon < 0. && xcrab < xmean) ||
	   (asylon > 0. && xcrab > xmean) )
	continue;

      AlphaMC->Fill(alpha);
      AbsAlphaMC->Fill(abs(alpha));

    }

  cout << "MC gammas in ON region: "  << MC_gammas_in_on << endl;
  cout << "MC gammas in OFF region: " << MC_gammas_in_off << " / 3" << endl;

  WidthMC->SetLineColor(4);
  LengthMC->SetLineColor(4);
  SWidthMC->SetLineColor(4);
  SLengthMC->SetLineColor(4);
  HadronnessMC->SetLineColor(4);
  ConcMC->SetLineColor(4);
  DistMC->SetLineColor(4);
  M3LongMC->SetLineColor(4);
  M3LongMC_all->SetLineColor(2);
  SizeMC->SetLineColor(4);

  /////////////////////////    DRAW    /////////////////////////////////


  WidthDiff->Draw();
  WidthMC->Draw("chisto,same");

  cout << "Data rate: " << WidthDiff->Integral() << " Hz" << endl;
  cout << "MC   rate: " << WidthMC->Integral() << " Hz" << endl;
  cout << "MC   events: " << WidthMC->GetEntries() << endl;

  hfile->Write();

  return;
}

///////////////////////////////////////

void Smooth(Char_t* filename="Hists_above2000_had0.15.root")
{
  TH2F* SkyMap = new TH2F();

  TFile *f = new TFile(filename);
  SkyMap->Read("SkyMapSubtracted");


  Int_t binfactor = 5;
  Int_t range = 11;     // odd number both!

  Int_t nx = binfactor*SkyMap->GetNbinsX();
  Int_t ny = binfactor*SkyMap->GetNbinsY();

  TH2F* SkyMap2 = new TH2F("SkyMap2", "Events incoming direction map ",
			    nx, 
			    SkyMap->GetXaxis()->GetXmin(),
			    SkyMap->GetXaxis()->GetXmax(),
			    ny, 
			    SkyMap->GetYaxis()->GetXmin(),
			    SkyMap->GetYaxis()->GetXmax());


  // Parameters from fit to real Crab data

  TF2* angres = new TF2("g", 
			"exp(-(x/0.01065)**2-(y/0.135684)**2)", 
			-SkyMap->GetXaxis()->GetBinWidth(1)*range/2., 
			SkyMap->GetXaxis()->GetBinWidth(1)*range/2.,
			-SkyMap->GetYaxis()->GetBinWidth(1)*range/2., 
			SkyMap->GetYaxis()->GetBinWidth(1)*range/2. );

  Int_t nbinsx = (angres->GetXmax() - angres->GetXmin()) / 
    (SkyMap2->GetXaxis()->GetBinWidth(1));

  Int_t nbinsy = (angres->GetYmax() - angres->GetYmin()) / 
    (SkyMap2->GetYaxis()->GetBinWidth(1));

  angres->SetNpx(nbinsx);
  angres->SetNpy(nbinsy);

  angres->Draw("lego1");
 

  TH2F* hist = (TH2F*)angres->GetHistogram()->Clone("hist");

  Float_t norm = 1./hist->Integral();

  hist->Scale(norm);

  cout << hist->Integral() << endl;

  Float_t maxx = SkyMap2->GetXaxis()->GetXmax();
  Float_t minx = SkyMap2->GetXaxis()->GetXmin();
  Float_t maxy = SkyMap2->GetYaxis()->GetXmax();
  Float_t miny = SkyMap2->GetYaxis()->GetXmin();

  for (Int_t ix = 1; ix <= SkyMap->GetNbinsX(); ix++)
    {
      for (Int_t iy = 1; iy <= SkyMap->GetNbinsY(); iy++)
	{
	  Float_t events = SkyMap->GetBinContent(ix, iy);

	  for (Int_t ix2 = 1; ix2 <= nbinsx; ix2++)
	    for (Int_t iy2 = 1; iy2 <= nbinsy; iy2++)
	      {
		Float_t events2 = events * hist->GetBinContent(ix2, iy2);

		Float_t newx = SkyMap->GetXaxis()->GetBinCenter(ix)
		  + hist->GetXaxis()->GetBinCenter(ix2);
		Float_t newy = SkyMap->GetYaxis()->GetBinCenter(iy) 
		  + hist->GetYaxis()->GetBinCenter(iy2);

		SkyMap2->Fill(newx, newy, events2);
	    }
	}
      if (ix%10 == 0)
	cout << ix << endl;
    }

  gStyle->SetPalette(1,0);
  SkyMap2->SetStats(0);
  SkyMap2->Draw("zcol");

  return;
}
