Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 6395) +++ /trunk/MagicSoft/Mars/Changelog (revision 6396) @@ -21,4 +21,9 @@ -*-*- END OF LINE -*-*- + + 2005/02/12 Abelardo Moralejo + + * mtemp/mpadova/macros/wobblemap.C + - Added. Macro to do some studies on wobble mode data. 2005/02/12 Markus Gaug Index: /trunk/MagicSoft/Mars/mtemp/mpadova/macros/wobblemap.C =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpadova/macros/wobblemap.C (revision 6396) +++ /trunk/MagicSoft/Mars/mtemp/mpadova/macros/wobblemap.C (revision 6396) @@ -0,0 +1,982 @@ +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; +}