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; }