diff --git a/analysis/AnaLinkDef.h b/analysis/AnaLinkDef.h index a726522b8..459f434b9 100644 --- a/analysis/AnaLinkDef.h +++ b/analysis/AnaLinkDef.h @@ -47,5 +47,6 @@ #pragma link C++ class R3BFibervsTofDOnlineSpectra+; #pragma link C++ class R3BGeneralOnlineSpectra+; #pragma link C++ class R3BIncomingIDOnlineSpectra+; +#pragma link C++ class R3BFiberTrackingOnlineSpectra+; #endif diff --git a/analysis/online/R3BFiberTrackingOnlineSpectra.cxx b/analysis/online/R3BFiberTrackingOnlineSpectra.cxx new file mode 100644 index 000000000..c6b70cac9 --- /dev/null +++ b/analysis/online/R3BFiberTrackingOnlineSpectra.cxx @@ -0,0 +1,2267 @@ +/****************************************************************************** + * Copyright (C) 2019 GSI Helmholtzzentrum für Schwerionenforschung GmbH * + * Copyright (C) 2019 Members of R3B Collaboration * + * * + * This software is distributed under the terms of the * + * GNU General Public Licence (GPL) version 3, * + * copied verbatim in the file "LICENSE". * + * * + * In applying this license GSI does not waive the privileges and immunities * + * granted to it by virtue of its status as an Intergovernmental Organization * + * or submit itself to any jurisdiction. * + ******************************************************************************/ + +// ------------------------------------------------------------ +// ----- R3BFiberTrackingOnlineSpectra ----- +// ----- Created 18/05/22 by J.L. Rodriguez-Sanchez & Enis Lorenz ----- +// ----- Fill tracking online histograms ----- +// ------------------------------------------------------------ + +/* + * This task should fill histograms for tracking after GLAD + */ + +#include "R3BFiberTrackingOnlineSpectra.h" +#include "R3BEventHeader.h" +#include "R3BLogger.h" +#include "R3BFiberMAPMTHitData.h" +#include "R3BFiberMappingPar.h" +#include "R3BTofdHitData.h" +#include "R3BLosCalData.h" +//#include "R3BMwpcHitData.h" +#include "R3BTGeoPar.h" +//#include "R3BTwimHitData.h" +#include "R3BFrsData.h" + +#include "FairLogger.h" +#include "FairRootManager.h" +#include "FairRunAna.h" +#include "FairRunOnline.h" +#include "FairRuntimeDb.h" + +#include "TArrow.h" +#include "TStyle.h" +#include "TColor.h" +#include "TCanvas.h" +#include "TClonesArray.h" +#include "TFolder.h" +#include "TH1F.h" +#include "TH2F.h" +#include "TH3D.h" +#include "TNtuple.h" +#include "THttpServer.h" +#include "TLatex.h" +#include "TLine.h" +#include "TMath.h" +#include "TRandom.h" +#include "TVector3.h" + +#define IS_NAN(x) TMath::IsNaN(x) + +R3BFiberTrackingOnlineSpectra::R3BFiberTrackingOnlineSpectra() + : R3BFiberTrackingOnlineSpectra("FiberTrackingOnlineSpectra", 1) +{ +} + +R3BFiberTrackingOnlineSpectra::R3BFiberTrackingOnlineSpectra(const TString& name, Int_t iVerbose) + : FairTask(name, iVerbose) + , fNEvents(0) + , fTrigger(-1) + , fTPatMask(0xffff)//16bits, all TPats + , fFi31MapPar(NULL) + , fFi32MapPar(NULL) + , fFi33MapPar(NULL) + , fFi30HitDataCA(NULL) + , fFi31HitDataCA(NULL) + , fFi32HitDataCA(NULL) + , fFi33HitDataCA(NULL) + , fHitTofdItems(NULL) + , fNbFibers31(512) + , fNbFibers32(512) + , fNbFibers33(512) + , fNbTofdPlanes(1) + , fNbTofdPaddlesPerPlane(44) + , fDist_acelerator_glad(3718.0) + , fPosTarget(2656.) + , fWidthTarget(40.) + , fZ_max(40.) + , fZ_min(0.) + , fFRSGATE(false) + , fAQmin(0.) + , fAQmax(100.) + , fQmin(0.) + , fQmax(100.) + , fTpat1(-1) + , fTpat2(-1) +{ +} + +R3BFiberTrackingOnlineSpectra::~R3BFiberTrackingOnlineSpectra() +{ + R3BLOG(debug1, "Destructor"); + if (fFi31HitDataCA) + { + delete fFi31HitDataCA; + } + if (fFi32HitDataCA) + { + delete fFi32HitDataCA; + } + if (fFi33HitDataCA) + { + delete fFi33HitDataCA; + } + if(fFi30HitDataCA) + { + delete fFi30HitDataCA; + } + if (fHitFrs) + delete fHitFrs; +} + +// ----- Public method SetParContainers -------------------------------- +void R3BFiberTrackingOnlineSpectra::SetParContainers() +{ + // Parameter Container + FairRuntimeDb* rtdb = FairRuntimeDb::instance(); + R3BLOG_IF(fatal, NULL == rtdb, "FairRuntimeDb not found"); + + // fMw0GeoPar = (R3BTGeoPar*)rtdb->getContainer("FiberGeoPar"); + // R3BLOG_IF(error, !FiberGeoPar, "Could not get access to FiberGeoPar container."); + + // fTargetGeoPar = (R3BTGeoPar*)rtdb->getContainer("TargetGeoPar"); + // R3BLOG_IF(error, !fTargetGeoPar, "Could not get access to TargetGeoPar container."); + + fMw1GeoPar = (R3BTGeoPar*)rtdb->getContainer("Mwpc1GeoPar"); + R3BLOG_IF(error, !fMw1GeoPar, "Could not get access to Mwpc1GeoPar container."); + return; + + fFi31MapPar = (R3BFiberMappingPar*)FairRuntimeDb::instance()->getContainer("Fi31MappingPar"); + R3BLOG_IF(error, !fFi31MapPar, "Couldn't get " << "Fi31" << "MappingPar"); + if (fFi31MapPar) + { + fNbFibers31 = fFi31MapPar->GetNbChannels(); + R3BLOG(info, "Fi31 " << "MappingPar found with " << fNbFibers31 << " fibers"); + } + fFi32MapPar = (R3BFiberMappingPar*)FairRuntimeDb::instance()->getContainer("Fi32MappingPar"); + R3BLOG_IF(error, !fFi32MapPar, "Couldn't get " << "Fi32" << "MappingPar"); + if (fFi32MapPar) + { + fNbFibers32 = fFi32MapPar->GetNbChannels(); + R3BLOG(info, "Fi32 " << "MappingPar found with " << fNbFibers32 << " fibers"); + } + fFi33MapPar = (R3BFiberMappingPar*)FairRuntimeDb::instance()->getContainer("Fi33MappingPar"); + R3BLOG_IF(error, !fFi33MapPar, "Couldn't get " << "Fi33" << "MappingPar"); + if (fFi33MapPar) + { + fNbFibers33 = fFi33MapPar->GetNbChannels(); + R3BLOG(info, "Fi33 " << "MappingPar found with " << fNbFibers33 << " fibers"); + } + +} + +// ----- Public method ReInit ---------------------------------------------- +InitStatus R3BFiberTrackingOnlineSpectra::ReInit() +{ + SetParContainers(); + return kSUCCESS; +} + +InitStatus R3BFiberTrackingOnlineSpectra::Init() +{ + R3BLOG(info, "For Fibers 31, 32, 33"); + FairRootManager* mgr = FairRootManager::Instance(); + R3BLOG_IF(fatal, NULL == mgr, "FairRootManager not found"); + + header = (R3BEventHeader*)mgr->GetObject("EventHeader."); + if (!header) + { + R3BLOG(warning, "EventHeader. not found"); + header = (R3BEventHeader*)mgr->GetObject("R3BEventHeader"); + } + else + R3BLOG(info, " EventHeader. found"); + + // uncomment lines below when ucesb avaliable + // FairRunOnline* run = FairRunOnline::Instance(); + // run->GetHttpServer()->Register("", this); + + fFi31HitDataCA = (TClonesArray*)mgr->GetObject("Fi31Hit"); + if (!fFi31HitDataCA) + { + R3BLOG(warning, "Fi31HitData not found"); + // return kwarning; + } + + fFi32HitDataCA = (TClonesArray*)mgr->GetObject("Fi32Hit"); + if (!fFi32HitDataCA) + { + R3BLOG(warning, "Fi32HitData not found"); + // return kwarning; + } + + fFi33HitDataCA = (TClonesArray*)mgr->GetObject("Fi33Hit"); + if (!fFi33HitDataCA) + { + R3BLOG(warning, "Fi33HitData not found"); + // return kwarning; + } + fFi30HitDataCA = (TClonesArray*)mgr->GetObject("Fi30Hit"); + if(!fFi30HitDataCA) + { + R3BLOG(warning, "Fi30HitData not found"); + } + + fHitTofdItems = (TClonesArray*)mgr->GetObject("TofdHit"); + R3BLOG_IF(warning, NULL == fHitTofdItems, "TofdHit not found"); + + fHitFrs = (TClonesArray*)mgr->GetObject("FrsData"); + R3BLOG_IF(warning, !fHitFrs, "Branch FrsData not found"); + + fCalLos = (TClonesArray*)mgr->GetObject("LosCal"); + R3BLOG_IF(warning, !fCalLos, "Branch LosCal not found"); + + // Create histograms for detectors + TString Name1; + TString Name2; + + cFibInfo = new TCanvas("FibInfo", "FibInfo", 10, 10, 800, 700); + + fh_ToT_Fib = new TH2F( + "Fib32Hit_ToT_iFib", "Fib32Hit ToT vs Fib", 512, 1., 512 + 1, 200, 0., 100.); + fh_ToT_Fib->GetXaxis()->SetTitle("Fiber number"); + fh_ToT_Fib->GetYaxis()->SetTitle("ToT / ns"); + + gPad->SetLogz(); + fh_ToT_Fib->Draw("colz"); + + // Hit data, fiber tracking plane X-Z + cTrackingXZ = new TCanvas("Tracking_after_GLAD_XZ", "Fragment Tracking (Fibers) plane XZ Fibers + ToFD", 10, 10, 800, 700); + + Name1 = "fh2_fibtracking_planeXZ"; + Name2 = "Fragment Tracking (Fibers) plane XZ Fibers + ToFD"; + Int_t histoYlim = 150; + fh2_fibtracking_planeXZ = new TH2F(Name1, Name2, 250, 0., 1250, 1400, -1. * 70., 70.); + fh2_fibtracking_planeXZ->GetXaxis()->SetTitle("Beam direction-Z [cm]"); + fh2_fibtracking_planeXZ->GetYaxis()->SetTitle("X Position [cm]"); + fh2_fibtracking_planeXZ->GetYaxis()->SetTitleOffset(1.1); + fh2_fibtracking_planeXZ->GetXaxis()->CenterTitle(true); + fh2_fibtracking_planeXZ->GetYaxis()->CenterTitle(true); + fh2_fibtracking_planeXZ->GetXaxis()->SetLabelSize(0.045); + fh2_fibtracking_planeXZ->GetXaxis()->SetTitleSize(0.045); + fh2_fibtracking_planeXZ->GetYaxis()->SetLabelSize(0.045); + fh2_fibtracking_planeXZ->GetYaxis()->SetTitleSize(0.045); + gPad->SetLogz(); + fh2_fibtracking_planeXZ->Draw("colz"); + + // Target indicated as a hole + TLine* l1 = new TLine(594.5, -0.6, 594.5, 50.6);//Fi33 + l1->SetLineColor(2); + l1->SetLineWidth(2); + l1->Draw(); + TLine* l3 = new TLine(611.4, 0.6, 611.4, -50.6);//Fi31 + l3->SetLineColor(2); + l3->SetLineWidth(2); + l3->Draw(); + TLine* l2 = new TLine(434.1, -25.6, 434.1, 25.6);//Fi32 + l2->SetLineColor(2); + l2->SetLineWidth(2); + l2->Draw(); + + TLine* l4 = new TLine(0., 0., 1250., 0.);//14° line + l4->SetLineColor(1); + l4->SetLineWidth(2); + l4->Draw(); + + TLine* l5 = new TLine(1138.5, -60., 1138.5, 60);//ToFD + l5->SetLineColor(1); + l5->SetLineWidth(2); + l5->Draw(); + + TArrow* arrow = new TArrow(100., 58., 300., 58., 0.02, ">"); + arrow->SetLineColor(3); + arrow->SetFillStyle(1001); + arrow->SetLineWidth(3); + arrow->Draw(); + + TLatex latex; + latex.SetTextSize(0.045); + latex.SetTextAlign(13); + latex.SetTextColor(3); + latex.DrawLatex(100., 65., "Beam"); + latex.SetTextColor(2); + latex.DrawLatex(434.1, -26.5, "Fi32"); + latex.SetTextColor(2); + latex.DrawLatex(631.4, -45., "Fi31"); + latex.SetTextColor(2); + latex.DrawLatex(604.5, 50., "Fi33"); + latex.SetTextColor(1); + latex.DrawLatex(1150., -50., "ToFD"); + latex.SetTextColor(1); + latex.DrawLatex(50., -2., "14 deg"); + latex.SetTextColor(1); + latex.DrawLatex(50., -50., "Umlenkpunkt"); + + + + cslopeVsX = new TCanvas("slopeX_vs_positionX_Fibers", "slope_XZ vs position X of Fibers", 10, 10, 800, 700); + fh2_fib_slope_vs_x = + new TH2F("slopeXvsPosX", "slope vs position on Fibers", 600, -60., 60., 1000, -50., 50.); + fh2_fib_slope_vs_x->GetXaxis()->SetTitle("RPC <-- X Position [cm] --> NeuLand"); + fh2_fib_slope_vs_x->GetYaxis()->SetTitle("NL <- slope rel. to 14 deg [mrad] -> RPC"); + fh2_fib_slope_vs_x->GetYaxis()->SetTitleOffset(1.1); + fh2_fib_slope_vs_x->GetXaxis()->CenterTitle(true); + fh2_fib_slope_vs_x->GetYaxis()->CenterTitle(true); + fh2_fib_slope_vs_x->GetXaxis()->SetLabelSize(0.045); + fh2_fib_slope_vs_x->GetXaxis()->SetTitleSize(0.04); + fh2_fib_slope_vs_x->GetYaxis()->SetLabelSize(0.045); + fh2_fib_slope_vs_x->GetYaxis()->SetTitleSize(0.04); + + fh2_fib_slope_vs_x4 = + new TH2F("slopeXvsPosX4", "slope vs position on Fibers Charge 4", 600, -60., 60., 1000, -50., 50.); + fh2_fib_slope_vs_x4->GetXaxis()->SetTitle("RPC <-- X Position [cm] --> NeuLand"); + fh2_fib_slope_vs_x4->GetYaxis()->SetTitle("NL <- slope rel. to 14 deg [mrad] -> RPC"); + fh2_fib_slope_vs_x4->GetYaxis()->SetTitleOffset(1.1); + fh2_fib_slope_vs_x4->GetXaxis()->CenterTitle(true); + fh2_fib_slope_vs_x4->GetYaxis()->CenterTitle(true); + fh2_fib_slope_vs_x4->GetXaxis()->SetLabelSize(0.045); + fh2_fib_slope_vs_x4->GetXaxis()->SetTitleSize(0.04); + fh2_fib_slope_vs_x4->GetYaxis()->SetLabelSize(0.045); + fh2_fib_slope_vs_x4->GetYaxis()->SetTitleSize(0.04); + + fh2_fib_slope_vs_x5 = + new TH2F("slopeXvsPosX5", "slope vs position on Fibers Charge 5", 600, -60., 60., 1000, -50., 50.); + fh2_fib_slope_vs_x5->GetXaxis()->SetTitle("RPC <-- X Position [cm] --> NeuLand"); + fh2_fib_slope_vs_x5->GetYaxis()->SetTitle("NL <- slope rel. to 14 deg [mrad] -> RPC"); + fh2_fib_slope_vs_x5->GetYaxis()->SetTitleOffset(1.1); + fh2_fib_slope_vs_x5->GetXaxis()->CenterTitle(true); + fh2_fib_slope_vs_x5->GetYaxis()->CenterTitle(true); + fh2_fib_slope_vs_x5->GetXaxis()->SetLabelSize(0.045); + fh2_fib_slope_vs_x5->GetXaxis()->SetTitleSize(0.04); + fh2_fib_slope_vs_x5->GetYaxis()->SetLabelSize(0.045); + fh2_fib_slope_vs_x5->GetYaxis()->SetTitleSize(0.04); + + fh2_fib_slope_vs_x6 = + new TH2F("slopeXvsPosX6", "slope vs position on Fibers Charge 6", 600, -60., 60., 1000, -50., 50.); + fh2_fib_slope_vs_x6->GetXaxis()->SetTitle("RPC <-- X Position [cm] --> NeuLand"); + fh2_fib_slope_vs_x6->GetYaxis()->SetTitle("NL <- slope rel. to 14 deg [mrad] -> RPC"); + fh2_fib_slope_vs_x6->GetYaxis()->SetTitleOffset(1.1); + fh2_fib_slope_vs_x6->GetXaxis()->CenterTitle(true); + fh2_fib_slope_vs_x6->GetYaxis()->CenterTitle(true); + fh2_fib_slope_vs_x6->GetXaxis()->SetLabelSize(0.045); + fh2_fib_slope_vs_x6->GetXaxis()->SetTitleSize(0.04); + fh2_fib_slope_vs_x6->GetYaxis()->SetLabelSize(0.045); + fh2_fib_slope_vs_x6->GetYaxis()->SetTitleSize(0.04); + + cslopeVsX->Divide(2,2); + gStyle->SetPalette(kRainBow); + cslopeVsX->cd(1); + gPad->SetLogz(); + fh2_fib_slope_vs_x->Draw("colz,cp55"); + cslopeVsX->cd(2); + gPad->SetLogz(); + fh2_fib_slope_vs_x4->Draw("colz,cp55"); + cslopeVsX->cd(3); + gPad->SetLogz(); + fh2_fib_slope_vs_x5->Draw("colz,cp55"); + cslopeVsX->cd(4); + gPad->SetLogz(); + fh2_fib_slope_vs_x6->Draw("colz,cp55"); + + + + // AngleX and positionX on Fibers relative to center / 14° line - positive angles mean higher bending, vice versa + cAngVsX = new TCanvas("AngleX_vs_positionX_Fibers", "Angle_XZ vs position X of Fibers", 10, 10, 800, 700); + fh2_fib_ang_vs_x = + new TH2F("AngXvsPosX", "Angle vs position on Fibers", 600, -60., 60., 1000, -50., 50.); + fh2_fib_ang_vs_x->GetXaxis()->SetTitle("RPC <-- X Position [cm] --> NeuLand"); + fh2_fib_ang_vs_x->GetYaxis()->SetTitle("NL <- Angle rel. to 14 deg [mrad] -> RPC"); + fh2_fib_ang_vs_x->GetYaxis()->SetTitleOffset(1.1); + fh2_fib_ang_vs_x->GetXaxis()->CenterTitle(true); + fh2_fib_ang_vs_x->GetYaxis()->CenterTitle(true); + fh2_fib_ang_vs_x->GetXaxis()->SetLabelSize(0.045); + fh2_fib_ang_vs_x->GetXaxis()->SetTitleSize(0.04); + fh2_fib_ang_vs_x->GetYaxis()->SetLabelSize(0.045); + fh2_fib_ang_vs_x->GetYaxis()->SetTitleSize(0.04); + + fh2_fib_ang_vs_x4 = + new TH2F("AngXvsPosX4", "Angle vs position on Fibers Charge 4", 600, -60., 60., 1000, -50., 50.); + fh2_fib_ang_vs_x4->GetXaxis()->SetTitle("RPC <-- X Position [cm] --> NeuLand"); + fh2_fib_ang_vs_x4->GetYaxis()->SetTitle("NL <- Angle rel. to 14 deg [mrad] -> RPC"); + fh2_fib_ang_vs_x4->GetYaxis()->SetTitleOffset(1.1); + fh2_fib_ang_vs_x4->GetXaxis()->CenterTitle(true); + fh2_fib_ang_vs_x4->GetYaxis()->CenterTitle(true); + fh2_fib_ang_vs_x4->GetXaxis()->SetLabelSize(0.045); + fh2_fib_ang_vs_x4->GetXaxis()->SetTitleSize(0.04); + fh2_fib_ang_vs_x4->GetYaxis()->SetLabelSize(0.045); + fh2_fib_ang_vs_x4->GetYaxis()->SetTitleSize(0.04); + + fh2_fib_ang_vs_x5 = + new TH2F("AngXvsPosX5", "Angle vs position on Fibers Charge 5", 600, -60., 60., 1000, -50., 50.); + fh2_fib_ang_vs_x5->GetXaxis()->SetTitle("RPC <-- X Position [cm] --> NeuLand"); + fh2_fib_ang_vs_x5->GetYaxis()->SetTitle("NL <- Angle rel. to 14 deg [mrad] -> RPC"); + fh2_fib_ang_vs_x5->GetYaxis()->SetTitleOffset(1.1); + fh2_fib_ang_vs_x5->GetXaxis()->CenterTitle(true); + fh2_fib_ang_vs_x5->GetYaxis()->CenterTitle(true); + fh2_fib_ang_vs_x5->GetXaxis()->SetLabelSize(0.045); + fh2_fib_ang_vs_x5->GetXaxis()->SetTitleSize(0.04); + fh2_fib_ang_vs_x5->GetYaxis()->SetLabelSize(0.045); + fh2_fib_ang_vs_x5->GetYaxis()->SetTitleSize(0.04); + + fh2_fib_ang_vs_x6 = + new TH2F("AngXvsPosX6", "Angle vs position on Fibers Charge 6", 600, -60., 60., 1000, -50., 50.); + fh2_fib_ang_vs_x6->GetXaxis()->SetTitle("RPC <-- X Position [cm] --> NeuLand"); + fh2_fib_ang_vs_x6->GetYaxis()->SetTitle("NL <- Angle rel. to 14 deg [mrad] -> RPC"); + fh2_fib_ang_vs_x6->GetYaxis()->SetTitleOffset(1.1); + fh2_fib_ang_vs_x6->GetXaxis()->CenterTitle(true); + fh2_fib_ang_vs_x6->GetYaxis()->CenterTitle(true); + fh2_fib_ang_vs_x6->GetXaxis()->SetLabelSize(0.045); + fh2_fib_ang_vs_x6->GetXaxis()->SetTitleSize(0.04); + fh2_fib_ang_vs_x6->GetYaxis()->SetLabelSize(0.045); + fh2_fib_ang_vs_x6->GetYaxis()->SetTitleSize(0.04); + + cAngVsX->Divide(2,2); + cAngVsX->cd(1); + gPad->SetLogz(); + fh2_fib_ang_vs_x->Draw("colz,cp55"); + cAngVsX->cd(2); + gPad->SetLogz(); + fh2_fib_ang_vs_x4->Draw("colz,cp55"); + cAngVsX->cd(3); + gPad->SetLogz(); + fh2_fib_ang_vs_x5->Draw("colz,cp55"); + cAngVsX->cd(4); + gPad->SetLogz(); + fh2_fib_ang_vs_x6->Draw("colz,cp55"); + + + //----------------------------fh_ang_vs_pos_rot + cAngRot = new TCanvas("AngX_vs_positionX_Fibers rotated", "Angle_XZ vs position X of Fibers", 10, 10, 800, 700); + + fh_ang_vs_pos_rot = + new TH2F("AngXvsPosXRot", "Angle vs Position on Fibers rotated", 600, -60., 60., 1000, -50., 50.); + fh_ang_vs_pos_rot->GetXaxis()->SetTitle("Arb. X"); + fh_ang_vs_pos_rot->GetYaxis()->SetTitle("Arb. Y"); + fh_ang_vs_pos_rot->GetYaxis()->SetTitleOffset(1.1); + fh_ang_vs_pos_rot->GetXaxis()->CenterTitle(true); + fh_ang_vs_pos_rot->GetYaxis()->CenterTitle(true); + fh_ang_vs_pos_rot->GetXaxis()->SetLabelSize(0.045); + fh_ang_vs_pos_rot->GetXaxis()->SetTitleSize(0.04); + fh_ang_vs_pos_rot->GetYaxis()->SetLabelSize(0.045); + fh_ang_vs_pos_rot->GetYaxis()->SetTitleSize(0.04); + + fh_pid = + new TH2F("PID_QvsProj", "Charge vs Trajectory by angle vs position", 240, -60., 60., 100, 0., 10.); + fh_pid->GetXaxis()->SetTitle("Arb. X / projection of angular correlation against position"); + fh_pid->GetYaxis()->SetTitle("Charge (ToFD)"); + fh_pid->GetYaxis()->SetTitleOffset(1.1); + fh_pid->GetXaxis()->CenterTitle(true); + fh_pid->GetYaxis()->CenterTitle(true); + fh_pid->GetXaxis()->SetLabelSize(0.045); + fh_pid->GetXaxis()->SetTitleSize(0.04); + fh_pid->GetYaxis()->SetLabelSize(0.045); + fh_pid->GetYaxis()->SetTitleSize(0.04); + + fh_pid_slope = + new TH2F("PID_QvsProjSlope", "Charge vs Linear Fitted Trajectory", 600, -60., 60., 750, 0., 30.); + fh_pid_slope->GetXaxis()->SetTitle("Arb. X / projection of slope against intercept"); + fh_pid_slope->GetYaxis()->SetTitle("Charge (ToFD)"); + fh_pid_slope->GetYaxis()->SetTitleOffset(1.1); + fh_pid_slope->GetXaxis()->CenterTitle(true); + fh_pid_slope->GetYaxis()->CenterTitle(true); + fh_pid_slope->GetXaxis()->SetLabelSize(0.045); + fh_pid_slope->GetXaxis()->SetTitleSize(0.04); + fh_pid_slope->GetYaxis()->SetLabelSize(0.045); + fh_pid_slope->GetYaxis()->SetTitleSize(0.04); + + fh_pid_3d = new TH3D("pid3d","pid3d",400, -40., 40., 200, -10., 10., 80, 0., 8.); + fh_pid_3d->GetXaxis()->SetTitle("RPC <-- X Position [cm] --> NeuLand"); + fh_pid_3d->GetYaxis()->SetTitle("NL <- Angle rel. to 14 deg [mrad] -> RPC"); + fh_pid_3d->GetZaxis()->SetTitle("Charge (ToFD)"); + fh_pid_3d->GetYaxis()->SetTitleOffset(1.1); + fh_pid_3d->GetXaxis()->CenterTitle(true); + fh_pid_3d->GetYaxis()->CenterTitle(true); + fh_pid_3d->GetZaxis()->CenterTitle(true); + fh_pid_3d->GetXaxis()->SetLabelSize(0.045); + fh_pid_3d->GetXaxis()->SetTitleSize(0.04); + fh_pid_3d->GetYaxis()->SetLabelSize(0.045); + fh_pid_3d->GetYaxis()->SetTitleSize(0.04); + fh_pid_3d->GetZaxis()->SetLabelSize(0.045); + fh_pid_3d->GetZaxis()->SetTitleSize(0.04); + + + + cAngRot->Divide(2,2); + gStyle->SetPalette(kRainBow); + cAngRot->cd(2); + gPad->SetLogz(); + fh_ang_vs_pos_rot->Draw("colz,cp55"); + gStyle->SetPalette(kRainBow); + cAngRot->cd(1); + gPad->SetLogz(); + fh2_fib_ang_vs_x->Draw("colz"); + cAngRot->cd(3); + gPad->SetLogz(); + fh_pid->Draw("colz"); + cAngRot->cd(4); + gPad->SetLogy(); + //fh_pid_3d->Draw("col"); + + cPID_3D = new TCanvas("PID", "PID", 10, 10, 800, 700); + //ntuple = new TNtuple("ntuple","data from ascii file","x:y:z:count"); + gStyle->SetPalette(kRainBow); + //ntuple->Draw("x:y:z:cont","","COLZ"); + fh_pid_3d->Draw("colt"); + + + cSlopeRot = new TCanvas("Slope_vs_positionX_Fibers rotated", "Angle_XZ vs position X of Fibers", 10, 10, 800, 700); + + fh_slope_vs_pos_rot = + new TH2F("SlopevsPosXRot", "Slope vs Position on Fibers rotated", 600, -60., 60., 1000, -50., 50.); + fh_slope_vs_pos_rot->GetXaxis()->SetTitle("Arb. X"); + fh_slope_vs_pos_rot->GetYaxis()->SetTitle("Arb. Y"); + fh_slope_vs_pos_rot->GetYaxis()->SetTitleOffset(1.1); + fh_slope_vs_pos_rot->GetXaxis()->CenterTitle(true); + fh_slope_vs_pos_rot->GetYaxis()->CenterTitle(true); + fh_slope_vs_pos_rot->GetXaxis()->SetLabelSize(0.045); + fh_slope_vs_pos_rot->GetXaxis()->SetTitleSize(0.04); + fh_slope_vs_pos_rot->GetYaxis()->SetLabelSize(0.045); + fh_slope_vs_pos_rot->GetYaxis()->SetTitleSize(0.04); + + + cSlopeRot->Divide(2,2); + gStyle->SetPalette(kRainBow); + cSlopeRot->cd(2); + gPad->SetLogz(); + fh_slope_vs_pos_rot->Draw("colz,cp55"); + gStyle->SetPalette(kRainBow); + cSlopeRot->cd(1); + gPad->SetLogz(); + fh2_fib_slope_vs_x->Draw("colz"); + cSlopeRot->cd(3); + gPad->SetLogz(); + fh_pid_slope->Draw("colz"); + /* + // AngleY and positionY on the target position + auto cAPY = new TCanvas("AngleY_vs_positionY_target", "Angle_YZ vs position Y on target", 10, 10, 800, 700); + fh2_angvsposy = + new TH2F("AngYvsPosY", "Angle vs position on target", 500, -fWidthTarget, fWidthTarget, 500, -10., 10.); + fh2_angvsposy->GetXaxis()->SetTitle("Y [mm]"); + fh2_angvsposy->GetYaxis()->SetTitle("Angle plane_YZ [mrad]"); + fh2_angvsposy->GetYaxis()->SetTitleOffset(1.1); + fh2_angvsposy->GetXaxis()->CenterTitle(true); + fh2_angvsposy->GetYaxis()->CenterTitle(true); + fh2_angvsposy->GetXaxis()->SetLabelSize(0.045); + fh2_angvsposy->GetXaxis()->SetTitleSize(0.045); + fh2_angvsposy->GetYaxis()->SetLabelSize(0.045); + fh2_angvsposy->GetYaxis()->SetTitleSize(0.045); + cAPY->cd(); + fh2_angvsposy->Draw("colz"); +*/ + // Hit data, Z versus beta + + cToFD_Z = new TCanvas("Charge_ToFD", "Q in ToFD", 10, 10, 800, 700); + fh_ToFD_Charge = + new TH1F("Charge Number", "Counts", 1200, 0., 12.); + fh_ToFD_Charge->GetXaxis()->SetTitle("Charge-Z"); + fh_ToFD_Charge->GetYaxis()->SetTitle("Counts"); + fh_ToFD_Charge->GetYaxis()->SetTitleOffset(1.1); + fh_ToFD_Charge->GetXaxis()->CenterTitle(true); + fh_ToFD_Charge->GetYaxis()->CenterTitle(true); + fh_ToFD_Charge->GetXaxis()->SetLabelSize(0.045); + fh_ToFD_Charge->GetXaxis()->SetTitleSize(0.045); + fh_ToFD_Charge->GetYaxis()->SetLabelSize(0.045); + fh_ToFD_Charge->GetYaxis()->SetTitleSize(0.045); + cToFD_Z->cd(); + gPad->SetLogy(); + fh_ToFD_Charge->Draw("colz"); + + AngCorrFibToFD = new TCanvas("Angular_Correlation_Fib_ToFD","Angle fibers vs angle ToFD",10,10,800,700); + fh2_ang_Fib_vs_ToFD = new TH2F("Angles Fiber","Angles ToFD",400,-20.,20.,400,-20.,20.); + fh2_ang_Fib_vs_ToFD->GetXaxis()->SetTitle("Angle Fiber32 -> Fibers31/33 [mrad]"); + fh2_ang_Fib_vs_ToFD->GetYaxis()->SetTitle("Angle Fiber mean -> ToFD [mrad]"); + fh2_ang_Fib_vs_ToFD->GetYaxis()->SetTitleOffset(1.1); + fh2_ang_Fib_vs_ToFD->GetXaxis()->CenterTitle(true); + fh2_ang_Fib_vs_ToFD->GetYaxis()->CenterTitle(true); + fh2_ang_Fib_vs_ToFD->GetXaxis()->SetLabelSize(0.045); + fh2_ang_Fib_vs_ToFD->GetXaxis()->SetTitleSize(0.045); + fh2_ang_Fib_vs_ToFD->GetYaxis()->SetLabelSize(0.045); + fh2_ang_Fib_vs_ToFD->GetYaxis()->SetTitleSize(0.045); + //AngCorrFibToFD->cd(); + //gPad->SetLogz(); + //fh2_ang_Fib_vs_ToFD->Draw("colz"); + + fh2_ang_Fib_vs_ToFD_nocuts = new TH2F("Angles Fiber no cuts","Angles ToFD",400,-20.,20.,400,-20.,20.); + fh2_ang_Fib_vs_ToFD_nocuts->GetXaxis()->SetTitle("Angle Fiber32 -> Fibers31/33 [mrad]"); + fh2_ang_Fib_vs_ToFD_nocuts->GetYaxis()->SetTitle("Angle Fiber mean -> ToFD [mrad]"); + fh2_ang_Fib_vs_ToFD_nocuts->GetYaxis()->SetTitleOffset(1.1); + fh2_ang_Fib_vs_ToFD_nocuts->GetXaxis()->CenterTitle(true); + fh2_ang_Fib_vs_ToFD_nocuts->GetYaxis()->CenterTitle(true); + fh2_ang_Fib_vs_ToFD_nocuts->GetXaxis()->SetLabelSize(0.045); + fh2_ang_Fib_vs_ToFD_nocuts->GetXaxis()->SetTitleSize(0.045); + fh2_ang_Fib_vs_ToFD_nocuts->GetYaxis()->SetLabelSize(0.045); + fh2_ang_Fib_vs_ToFD_nocuts->GetYaxis()->SetTitleSize(0.045); + AngCorrFibToFD->Divide(2,1); + AngCorrFibToFD->cd(1); + gPad->SetLogz(); + fh2_ang_Fib_vs_ToFD_nocuts->Draw("colz"); + AngCorrFibToFD->cd(2); + gPad->SetLogz(); + fh2_ang_Fib_vs_ToFD->Draw("colz"); + + + + cFibToT_Q = new TCanvas("FibToT_Charge_ToFD", "FibToT vs Q ToFD", 10, 10, 800, 700); + + fh2_FibToT_vs_QToFD = new TH2F("FibToTvsQToFD", "FibToT vs QToFD", 400, 0., 40.,200, 0.,20.); + fh2_FibToT_vs_QToFD->GetXaxis()->SetTitle("Fib32 ToT [ns]"); + fh2_FibToT_vs_QToFD->GetYaxis()->SetTitle("Charge (ToFD)"); + fh2_FibToT_vs_QToFD->GetYaxis()->SetTitleOffset(1.1); + fh2_FibToT_vs_QToFD->GetXaxis()->CenterTitle(true); + fh2_FibToT_vs_QToFD->GetYaxis()->CenterTitle(true); + fh2_FibToT_vs_QToFD->GetXaxis()->SetLabelSize(0.045); + fh2_FibToT_vs_QToFD->GetXaxis()->SetTitleSize(0.045); + fh2_FibToT_vs_QToFD->GetYaxis()->SetLabelSize(0.045); + fh2_FibToT_vs_QToFD->GetYaxis()->SetTitleSize(0.045); + cFibToT_Q->cd(); + gPad->SetLogz(); + fh2_FibToT_vs_QToFD->Draw("colz"); + + cFibToTCorr = new TCanvas("FibToT_Corr", "Fib ToT Corr", 10, 10, 800, 700); + + fh2_ToT_Fib32_vs_Fib31 = new TH2F("ToTFib32vsFib31", "ToT Fib32 vs Fib31", 400, 0., 40.,400, 0., 40.); + fh2_ToT_Fib32_vs_Fib31->GetXaxis()->SetTitle("Fib32 ToT [ns]"); + fh2_ToT_Fib32_vs_Fib31->GetYaxis()->SetTitle("Fib31 ToT [ns]"); + fh2_ToT_Fib32_vs_Fib31->GetYaxis()->SetTitleOffset(1.1); + fh2_ToT_Fib32_vs_Fib31->GetXaxis()->CenterTitle(true); + fh2_ToT_Fib32_vs_Fib31->GetYaxis()->CenterTitle(true); + fh2_ToT_Fib32_vs_Fib31->GetXaxis()->SetLabelSize(0.045); + fh2_ToT_Fib32_vs_Fib31->GetXaxis()->SetTitleSize(0.045); + fh2_ToT_Fib32_vs_Fib31->GetYaxis()->SetLabelSize(0.045); + fh2_ToT_Fib32_vs_Fib31->GetYaxis()->SetTitleSize(0.045); + + fh2_ToT_Fib32_vs_Fib30 = new TH2F("ToTFib32vsFib30", "ToT Fib32 vs Fib30", 400, 0., 40.,400, 0., 40.); + fh2_ToT_Fib32_vs_Fib30->GetXaxis()->SetTitle("Fib32 ToT [ns]"); + fh2_ToT_Fib32_vs_Fib30->GetYaxis()->SetTitle("Fib30 ToT [ns]"); + fh2_ToT_Fib32_vs_Fib30->GetYaxis()->SetTitleOffset(1.1); + fh2_ToT_Fib32_vs_Fib30->GetXaxis()->CenterTitle(true); + fh2_ToT_Fib32_vs_Fib30->GetYaxis()->CenterTitle(true); + fh2_ToT_Fib32_vs_Fib30->GetXaxis()->SetLabelSize(0.045); + fh2_ToT_Fib32_vs_Fib30->GetXaxis()->SetTitleSize(0.045); + fh2_ToT_Fib32_vs_Fib30->GetYaxis()->SetLabelSize(0.045); + fh2_ToT_Fib32_vs_Fib30->GetYaxis()->SetTitleSize(0.045); + + fh2_ToT_Fib32_vs_Fib33 = new TH2F("ToTFib32vsFib33", "ToT Fib32 vs Fib33", 400, 0., 40.,400, 0., 40.); + fh2_ToT_Fib32_vs_Fib33->GetXaxis()->SetTitle("Fib32 ToT [ns]"); + fh2_ToT_Fib32_vs_Fib33->GetYaxis()->SetTitle("Fib33 ToT [ns]"); + fh2_ToT_Fib32_vs_Fib33->GetYaxis()->SetTitleOffset(1.1); + fh2_ToT_Fib32_vs_Fib33->GetXaxis()->CenterTitle(true); + fh2_ToT_Fib32_vs_Fib33->GetYaxis()->CenterTitle(true); + fh2_ToT_Fib32_vs_Fib33->GetXaxis()->SetLabelSize(0.045); + fh2_ToT_Fib32_vs_Fib33->GetXaxis()->SetTitleSize(0.045); + fh2_ToT_Fib32_vs_Fib33->GetYaxis()->SetLabelSize(0.045); + fh2_ToT_Fib32_vs_Fib33->GetYaxis()->SetTitleSize(0.045); + + fh_pos_TPat = new TH2F("X_vs_TPat", "X vs Tpat", 600, -60., 60.,17, 0., 17.); + fh_pos_TPat->GetXaxis()->SetTitle("X [cm]"); + fh_pos_TPat->GetYaxis()->SetTitle("TPat"); + fh_pos_TPat->GetYaxis()->SetTitleOffset(1.1); + fh_pos_TPat->GetXaxis()->CenterTitle(true); + fh_pos_TPat->GetYaxis()->CenterTitle(true); + fh_pos_TPat->GetXaxis()->SetLabelSize(0.045); + fh_pos_TPat->GetXaxis()->SetTitleSize(0.045); + fh_pos_TPat->GetYaxis()->SetLabelSize(0.045); + fh_pos_TPat->GetYaxis()->SetTitleSize(0.045); + + fh_x_vs_tof = new TH2F("X_vs_ToF", "X vs ToF", 600, -60., 60., 1000, -200., 0.); + fh_x_vs_tof->GetXaxis()->SetTitle("X [cm]"); + fh_x_vs_tof->GetYaxis()->SetTitle("ToF [ns]"); + fh_x_vs_tof->GetYaxis()->SetTitleOffset(1.1); + fh_x_vs_tof->GetXaxis()->CenterTitle(true); + fh_x_vs_tof->GetYaxis()->CenterTitle(true); + fh_x_vs_tof->GetXaxis()->SetLabelSize(0.045); + fh_x_vs_tof->GetXaxis()->SetTitleSize(0.045); + fh_x_vs_tof->GetYaxis()->SetLabelSize(0.045); + fh_x_vs_tof->GetYaxis()->SetTitleSize(0.045); + + fh_x_vs_q = new TH2F("X_vs_Q", "X vs Q", 600, -60., 60., 300, 4., 7.); + fh_x_vs_q->GetXaxis()->SetTitle("X [cm]"); + fh_x_vs_q->GetYaxis()->SetTitle("charge"); + fh_x_vs_q->GetYaxis()->SetTitleOffset(1.1); + fh_x_vs_q->GetXaxis()->CenterTitle(true); + fh_x_vs_q->GetYaxis()->CenterTitle(true); + fh_x_vs_q->GetXaxis()->SetLabelSize(0.045); + fh_x_vs_q->GetXaxis()->SetTitleSize(0.045); + fh_x_vs_q->GetYaxis()->SetLabelSize(0.045); + fh_x_vs_q->GetYaxis()->SetTitleSize(0.045); + + cFibToTCorr->Divide(2,2); + cFibToTCorr->cd(1); + gPad->SetLogz(); + fh2_ToT_Fib32_vs_Fib30->Draw("colz"); + cFibToTCorr->cd(2); + gPad->SetLogz(); + //fh2_ToT_Fib32_vs_Fib31->Draw("colz"); + fh_x_vs_q->Draw("colz"); + cFibToTCorr->cd(3); + gPad->SetLogz(); + //fh2_ToT_Fib32_vs_Fib33->Draw("colz"); + fh_x_vs_tof->Draw("colz"); + cFibToTCorr->cd(4); + gPad->SetLogz(); + fh_pos_TPat->Draw("colz"); + + cFibToFDnocuts = new TCanvas("NoCuts", "NoCuts", 10, 10, 800, 700); + fh2_fib_ang_vs_x_nocuts = + new TH2F("AngXvsPosXnocuts", "Angle vs position on Fibers no cuts", 600, -60., 60., 1000, -50., 50.); + fh2_fib_ang_vs_x_nocuts->GetXaxis()->SetTitle("RPC <-- X Position [cm] --> NeuLand"); + fh2_fib_ang_vs_x_nocuts->GetYaxis()->SetTitle("NL <- Angle rel. to 14 deg [mrad] -> RPC"); + fh2_fib_ang_vs_x_nocuts->GetYaxis()->SetTitleOffset(1.1); + fh2_fib_ang_vs_x_nocuts->GetXaxis()->CenterTitle(true); + fh2_fib_ang_vs_x_nocuts->GetYaxis()->CenterTitle(true); + fh2_fib_ang_vs_x_nocuts->GetXaxis()->SetLabelSize(0.045); + fh2_fib_ang_vs_x_nocuts->GetXaxis()->SetTitleSize(0.04); + fh2_fib_ang_vs_x_nocuts->GetYaxis()->SetLabelSize(0.045); + fh2_fib_ang_vs_x_nocuts->GetYaxis()->SetTitleSize(0.04); + + fh2_fib_slope_vs_x_nocuts = + new TH2F("slopeXvsPosXnocuts", "slope vs position on Fibers", 600, -60., 60., 1000, -50., 50.); + fh2_fib_slope_vs_x_nocuts->GetXaxis()->SetTitle("RPC <-- X Position [cm] --> NeuLand"); + fh2_fib_slope_vs_x_nocuts->GetYaxis()->SetTitle("NL <- slope rel. to 14 deg [mrad] -> RPC"); + fh2_fib_slope_vs_x_nocuts->GetYaxis()->SetTitleOffset(1.1); + fh2_fib_slope_vs_x_nocuts->GetXaxis()->CenterTitle(true); + fh2_fib_slope_vs_x_nocuts->GetYaxis()->CenterTitle(true); + fh2_fib_slope_vs_x_nocuts->GetXaxis()->SetLabelSize(0.045); + fh2_fib_slope_vs_x_nocuts->GetXaxis()->SetTitleSize(0.04); + fh2_fib_slope_vs_x_nocuts->GetYaxis()->SetLabelSize(0.045); + fh2_fib_slope_vs_x_nocuts->GetYaxis()->SetTitleSize(0.04); + + fh2_fibtracking_planeXZ_nocuts = new TH2F(Name1 +"nocuts", Name2, 250, 0., 1250, 1400, -1. * 70., 70.); + fh2_fibtracking_planeXZ_nocuts->GetXaxis()->SetTitle("Beam direction-Z [cm]"); + fh2_fibtracking_planeXZ_nocuts->GetYaxis()->SetTitle("X Position [cm]"); + fh2_fibtracking_planeXZ_nocuts->GetYaxis()->SetTitleOffset(1.1); + fh2_fibtracking_planeXZ_nocuts->GetXaxis()->CenterTitle(true); + fh2_fibtracking_planeXZ_nocuts->GetYaxis()->CenterTitle(true); + fh2_fibtracking_planeXZ_nocuts->GetXaxis()->SetLabelSize(0.045); + fh2_fibtracking_planeXZ_nocuts->GetXaxis()->SetTitleSize(0.045); + fh2_fibtracking_planeXZ_nocuts->GetYaxis()->SetLabelSize(0.045); + fh2_fibtracking_planeXZ_nocuts->GetYaxis()->SetTitleSize(0.045); + + fh_Rsqr = new TH1F("R_Sqr", "R Squared", 1100, 0., 1.1); + fh_Rsqr->GetXaxis()->SetTitle("R Squared"); + fh_Rsqr->GetYaxis()->SetTitle("Counts"); + fh_Rsqr->GetYaxis()->SetTitleOffset(1.1); + fh_Rsqr->GetXaxis()->CenterTitle(true); + fh_Rsqr->GetYaxis()->CenterTitle(true); + fh_Rsqr->GetXaxis()->SetLabelSize(0.045); + fh_Rsqr->GetXaxis()->SetTitleSize(0.045); + fh_Rsqr->GetYaxis()->SetLabelSize(0.045); + fh_Rsqr->GetYaxis()->SetTitleSize(0.045); + + cFibToFDnocuts->Divide(2,2); + cFibToFDnocuts->cd(1); + gPad->SetLogz(); + fh2_fib_ang_vs_x_nocuts->Draw("colz"); + + cFibToFDnocuts->cd(2); + gPad->SetLogz(); + fh2_fib_slope_vs_x_nocuts->Draw("colz"); + cFibToFDnocuts->cd(3); + gPad->SetLogz(); + fh2_fibtracking_planeXZ_nocuts->Draw("colz"); + cFibToFDnocuts->cd(4); + gPad->SetLogy(); + fh_Rsqr->Draw(); + + + cFib30Rel = new TCanvas("Fib30Rel", "Fib30Rel", 10, 10, 800, 700); + + fh_f32X_vs_f30Y = new TH2F("Fib32XvsFib30Y", "Fib32X vs Fib30Y", 600, -30., 30.,600, -30., 30.); + fh_f32X_vs_f30Y->GetXaxis()->SetTitle("Fib32 X [cm]"); + fh_f32X_vs_f30Y->GetYaxis()->SetTitle("Fib30 Y [cm]"); + fh_f32X_vs_f30Y->GetYaxis()->SetTitleOffset(1.1); + fh_f32X_vs_f30Y->GetXaxis()->CenterTitle(true); + fh_f32X_vs_f30Y->GetYaxis()->CenterTitle(true); + fh_f32X_vs_f30Y->GetXaxis()->SetLabelSize(0.045); + fh_f32X_vs_f30Y->GetXaxis()->SetTitleSize(0.045); + fh_f32X_vs_f30Y->GetYaxis()->SetLabelSize(0.045); + fh_f32X_vs_f30Y->GetYaxis()->SetTitleSize(0.045); + + fh_f32Y_vs_f30Y = new TH2F("Fib32YvsFib30Y", "Fib32Y vs Fib30Y", 600, -30., 30.,600, -30., 30.); + fh_f32Y_vs_f30Y->GetXaxis()->SetTitle("Fib32 Y [cm]"); + fh_f32Y_vs_f30Y->GetYaxis()->SetTitle("Fib30 Y [cm]"); + fh_f32Y_vs_f30Y->GetYaxis()->SetTitleOffset(1.1); + fh_f32Y_vs_f30Y->GetXaxis()->CenterTitle(true); + fh_f32Y_vs_f30Y->GetYaxis()->CenterTitle(true); + fh_f32Y_vs_f30Y->GetXaxis()->SetLabelSize(0.045); + fh_f32Y_vs_f30Y->GetXaxis()->SetTitleSize(0.045); + fh_f32Y_vs_f30Y->GetYaxis()->SetLabelSize(0.045); + fh_f32Y_vs_f30Y->GetYaxis()->SetTitleSize(0.045); + + fh_f32tot_vs_f30Y = new TH2F("Fib32totvsFib30Y", "Fib32ToT vs Fib30Y", 400, 0., 40.,600, -30., 30.); + fh_f32tot_vs_f30Y->GetXaxis()->SetTitle("Fib32 ToT [ns]"); + fh_f32tot_vs_f30Y->GetYaxis()->SetTitle("Fib30 Y [cm]"); + fh_f32tot_vs_f30Y->GetYaxis()->SetTitleOffset(1.1); + fh_f32tot_vs_f30Y->GetXaxis()->CenterTitle(true); + fh_f32tot_vs_f30Y->GetYaxis()->CenterTitle(true); + fh_f32tot_vs_f30Y->GetXaxis()->SetLabelSize(0.045); + fh_f32tot_vs_f30Y->GetXaxis()->SetTitleSize(0.045); + fh_f32tot_vs_f30Y->GetYaxis()->SetLabelSize(0.045); + fh_f32tot_vs_f30Y->GetYaxis()->SetTitleSize(0.045); + + fh_tofdq_vs_f30Y = new TH2F("ToFDQvsFib30Y", "ToFDQ vs Fib30Y", 100, 0., 10.,600, -30., 30.); + fh_tofdq_vs_f30Y->GetXaxis()->SetTitle("ToFD Charge"); + fh_tofdq_vs_f30Y->GetYaxis()->SetTitle("Fib30 Y [cm]"); + fh_tofdq_vs_f30Y->GetYaxis()->SetTitleOffset(1.1); + fh_tofdq_vs_f30Y->GetXaxis()->CenterTitle(true); + fh_tofdq_vs_f30Y->GetYaxis()->CenterTitle(true); + fh_tofdq_vs_f30Y->GetXaxis()->SetLabelSize(0.045); + fh_tofdq_vs_f30Y->GetXaxis()->SetTitleSize(0.045); + fh_tofdq_vs_f30Y->GetYaxis()->SetLabelSize(0.045); + fh_tofdq_vs_f30Y->GetYaxis()->SetTitleSize(0.045); + + cFib30Rel->Divide(2,2); + cFib30Rel->cd(1); + gPad->SetLogz(); + fh_f32X_vs_f30Y->Draw("colz"); + + cFib30Rel->cd(2); + gPad->SetLogz(); + fh_f32Y_vs_f30Y->Draw("colz"); + cFib30Rel->cd(3); + gPad->SetLogz(); + fh_f32tot_vs_f30Y->Draw("colz"); + cFib30Rel->cd(4); + gPad->SetLogz(); + fh_tofdq_vs_f30Y->Draw("colz"); + + cXcorr = new TCanvas("Xcorr", "Xcorr", 10, 10, 800, 700); + + fh_X_fib_vs_tofd_no_cut = new TH2F("X_corr_fib_tofd_fit_nocut", "X corr fib tofd fit no cut", 600, -60., 60., 600, -60., 60.); + fh_X_fib_vs_tofd_no_cut->GetXaxis()->SetTitle("ToFD X pos [cm]"); + fh_X_fib_vs_tofd_no_cut->GetYaxis()->SetTitle("Extrap. X pos [cm]"); + fh_X_fib_vs_tofd_no_cut->GetYaxis()->SetTitleOffset(1.1); + fh_X_fib_vs_tofd_no_cut->GetXaxis()->CenterTitle(true); + fh_X_fib_vs_tofd_no_cut->GetYaxis()->CenterTitle(true); + fh_X_fib_vs_tofd_no_cut->GetXaxis()->SetLabelSize(0.045); + fh_X_fib_vs_tofd_no_cut->GetXaxis()->SetTitleSize(0.045); + fh_X_fib_vs_tofd_no_cut->GetYaxis()->SetLabelSize(0.045); + fh_X_fib_vs_tofd_no_cut->GetYaxis()->SetTitleSize(0.045); + + fh_X_fib_vs_tofd = new TH2F("X_corr_fib_tofd_fit", "X corr fib tofd fit", 600, -60., 60., 600, -60., 60.); + fh_X_fib_vs_tofd->GetXaxis()->SetTitle("ToFD X pos [cm]"); + fh_X_fib_vs_tofd->GetYaxis()->SetTitle("Extrap. X pos [cm]"); + fh_X_fib_vs_tofd->GetYaxis()->SetTitleOffset(1.1); + fh_X_fib_vs_tofd->GetXaxis()->CenterTitle(true); + fh_X_fib_vs_tofd->GetYaxis()->CenterTitle(true); + fh_X_fib_vs_tofd->GetXaxis()->SetLabelSize(0.045); + fh_X_fib_vs_tofd->GetXaxis()->SetTitleSize(0.045); + fh_X_fib_vs_tofd->GetYaxis()->SetLabelSize(0.045); + fh_X_fib_vs_tofd->GetYaxis()->SetTitleSize(0.045); + + fh_X_fib_vs_tofd_ang = new TH2F("X_corr_fib_tofd_ang", "X corr fib tofd ang", 600, -60., 60., 600, -60., 60.); + fh_X_fib_vs_tofd_ang->GetXaxis()->SetTitle("ToFD X pos [cm]"); + fh_X_fib_vs_tofd_ang->GetYaxis()->SetTitle("Extrap. X pos [cm]"); + fh_X_fib_vs_tofd_ang->GetYaxis()->SetTitleOffset(1.1); + fh_X_fib_vs_tofd_ang->GetXaxis()->CenterTitle(true); + fh_X_fib_vs_tofd_ang->GetYaxis()->CenterTitle(true); + fh_X_fib_vs_tofd_ang->GetXaxis()->SetLabelSize(0.045); + fh_X_fib_vs_tofd_ang->GetXaxis()->SetTitleSize(0.045); + fh_X_fib_vs_tofd_ang->GetYaxis()->SetLabelSize(0.045); + fh_X_fib_vs_tofd_ang->GetYaxis()->SetTitleSize(0.045); + + fh_X_fib_vs_tofd_ang_no_cut = new TH2F("X_corr_fib_tofd_ang_no_cut", "X corr fib tofd ang no cut", 600, -60., 60., 600, -60., 60.); + fh_X_fib_vs_tofd_ang_no_cut->GetXaxis()->SetTitle("ToFD X pos [cm]"); + fh_X_fib_vs_tofd_ang_no_cut->GetYaxis()->SetTitle("Extrap. X pos [cm]"); + fh_X_fib_vs_tofd_ang_no_cut->GetYaxis()->SetTitleOffset(1.1); + fh_X_fib_vs_tofd_ang_no_cut->GetXaxis()->CenterTitle(true); + fh_X_fib_vs_tofd_ang_no_cut->GetYaxis()->CenterTitle(true); + fh_X_fib_vs_tofd_ang_no_cut->GetXaxis()->SetLabelSize(0.045); + fh_X_fib_vs_tofd_ang_no_cut->GetXaxis()->SetTitleSize(0.045); + fh_X_fib_vs_tofd_ang_no_cut->GetYaxis()->SetLabelSize(0.045); + fh_X_fib_vs_tofd_ang_no_cut->GetYaxis()->SetTitleSize(0.045); + + cXcorr->Divide(2,2); + cXcorr->cd(1); + gPad->SetLogz(); + fh_X_fib_vs_tofd_no_cut->Draw("colz"); + cXcorr->cd(3); + gPad->SetLogz(); + fh_X_fib_vs_tofd->Draw("colz"); + cXcorr->cd(2); + gPad->SetLogz(); + fh_X_fib_vs_tofd_ang_no_cut->Draw("colz"); + cXcorr->cd(4); + gPad->SetLogz(); + fh_X_fib_vs_tofd_ang->Draw("colz"); + + + cDeltaQ = new TCanvas("ToFDDeltaQ", "ToFDDeltaQ", 10, 10, 800, 700); + cDeltaQ->Divide(3,2); + int g = 0; + for(int i = 0; i<4; i++){ + for(int j=i+1;j<4;j++){ + char strName1[255]; + sprintf(strName1, "tofd_Charge_plane_%d_vs_%d",i+1, j+1); + char strName2[255]; + sprintf(strName2, "tofd Charge plane %d vs %d",i+1, j+1); + fh_tofd_charge_planes[g] = new TH2F(strName1, strName2, 45, 0., 45., 200,-10.,10.); + fh_tofd_charge_planes[g]->GetXaxis()->SetTitle("Bar"); + fh_tofd_charge_planes[g]->GetYaxis()->SetTitle("Delta Q"); + fh_tofd_charge_planes[g]->GetYaxis()->SetTitleOffset(1.); + fh_tofd_charge_planes[g]->GetXaxis()->CenterTitle(true); + fh_tofd_charge_planes[g]->GetYaxis()->CenterTitle(true); + fh_tofd_charge_planes[g]->GetXaxis()->SetLabelSize(0.045); + fh_tofd_charge_planes[g]->GetXaxis()->SetTitleSize(0.045); + fh_tofd_charge_planes[g]->GetYaxis()->SetLabelSize(0.045); + fh_tofd_charge_planes[g]->GetYaxis()->SetTitleSize(0.045); + fh_tofd_charge_planes[g]->SetFillColor(31); + + cDeltaQ->cd(g+1); + gPad->SetLogz(); + fh_tofd_charge_planes[g]->Draw("colz"); + g++; + } + } + + cMult = new TCanvas("Mult", "Mult", 10, 10, 800, 700); + fh_mult_coinc = new TH1D("mult_coinc", "Multiplicity in coincidence", 151, -1, 150.); + fh_mult_coinc->GetXaxis()->SetTitle("Multiplicity"); + fh_mult_coinc->GetYaxis()->SetTitle("Counts"); + fh_mult_coinc->GetYaxis()->SetTitleOffset(1.1); + fh_mult_coinc->GetXaxis()->CenterTitle(true); + fh_mult_coinc->GetYaxis()->CenterTitle(true); + fh_mult_coinc->GetXaxis()->SetLabelSize(0.045); + fh_mult_coinc->GetXaxis()->SetTitleSize(0.045); + fh_mult_coinc->GetYaxis()->SetLabelSize(0.045); + fh_mult_coinc->GetYaxis()->SetTitleSize(0.045); + + fh_mult_minus_los = new TH1D("mult_coinc_LOS", "Multiplicity in coincidence minus LOS", 151, -1., 150.); + fh_mult_minus_los->GetXaxis()->SetTitle("Multiplicity"); + fh_mult_minus_los->GetYaxis()->SetTitle("Counts"); + fh_mult_minus_los->GetYaxis()->SetTitleOffset(1.1); + fh_mult_minus_los->GetXaxis()->CenterTitle(true); + fh_mult_minus_los->GetYaxis()->CenterTitle(true); + fh_mult_minus_los->GetXaxis()->SetLabelSize(0.045); + fh_mult_minus_los->GetXaxis()->SetTitleSize(0.045); + fh_mult_minus_los->GetYaxis()->SetLabelSize(0.045); + fh_mult_minus_los->GetYaxis()->SetTitleSize(0.045); + + fh_mult_fib30 = new TH1D("mult_fib30", "Multiplicity in fib30", 151, -1, 150.); + fh_mult_fib30->GetXaxis()->SetTitle("Multiplicity"); + fh_mult_fib30->GetYaxis()->SetTitle("Counts"); + fh_mult_fib30->GetYaxis()->SetTitleOffset(1.1); + fh_mult_fib30->GetXaxis()->CenterTitle(true); + fh_mult_fib30->GetYaxis()->CenterTitle(true); + fh_mult_fib30->GetXaxis()->SetLabelSize(0.045); + fh_mult_fib30->GetXaxis()->SetTitleSize(0.045); + fh_mult_fib30->GetYaxis()->SetLabelSize(0.045); + fh_mult_fib30->GetYaxis()->SetTitleSize(0.045); + + fh_mult_fib31 = new TH1D("mult_fib31", "Multiplicity in fib31", 151, -1, 150.); + fh_mult_fib31->GetXaxis()->SetTitle("Multiplicity"); + fh_mult_fib31->GetYaxis()->SetTitle("Counts"); + fh_mult_fib31->GetYaxis()->SetTitleOffset(1.1); + fh_mult_fib31->GetXaxis()->CenterTitle(true); + fh_mult_fib31->GetYaxis()->CenterTitle(true); + fh_mult_fib31->GetXaxis()->SetLabelSize(0.045); + fh_mult_fib31->GetXaxis()->SetTitleSize(0.045); + fh_mult_fib31->GetYaxis()->SetLabelSize(0.045); + fh_mult_fib31->GetYaxis()->SetTitleSize(0.045); + + fh_mult_fib32 = new TH1D("mult_fib32", "Multiplicity in fib32", 151, -1, 150.); + fh_mult_fib32->GetXaxis()->SetTitle("Multiplicity"); + fh_mult_fib32->GetYaxis()->SetTitle("Counts"); + fh_mult_fib32->GetYaxis()->SetTitleOffset(1.1); + fh_mult_fib32->GetXaxis()->CenterTitle(true); + fh_mult_fib32->GetYaxis()->CenterTitle(true); + fh_mult_fib32->GetXaxis()->SetLabelSize(0.045); + fh_mult_fib32->GetXaxis()->SetTitleSize(0.045); + fh_mult_fib32->GetYaxis()->SetLabelSize(0.045); + fh_mult_fib32->GetYaxis()->SetTitleSize(0.045); + + fh_mult_fib33 = new TH1D("mult_fib33", "Multiplicity in fib33", 151, -1, 150.); + fh_mult_fib33->GetXaxis()->SetTitle("Multiplicity"); + fh_mult_fib33->GetYaxis()->SetTitle("Counts"); + fh_mult_fib33->GetYaxis()->SetTitleOffset(1.1); + fh_mult_fib33->GetXaxis()->CenterTitle(true); + fh_mult_fib33->GetYaxis()->CenterTitle(true); + fh_mult_fib33->GetXaxis()->SetLabelSize(0.045); + fh_mult_fib33->GetXaxis()->SetTitleSize(0.045); + fh_mult_fib33->GetYaxis()->SetLabelSize(0.045); + fh_mult_fib33->GetYaxis()->SetTitleSize(0.045); + + fh_mult_fibAll = new TH1D("mult_fibAll", "Multiplicity in fibAll", 151, -1, 150.); + fh_mult_fibAll->GetXaxis()->SetTitle("Multiplicity"); + fh_mult_fibAll->GetYaxis()->SetTitle("Counts"); + fh_mult_fibAll->GetYaxis()->SetTitleOffset(1.1); + fh_mult_fibAll->GetXaxis()->CenterTitle(true); + fh_mult_fibAll->GetYaxis()->CenterTitle(true); + fh_mult_fibAll->GetXaxis()->SetLabelSize(0.045); + fh_mult_fibAll->GetXaxis()->SetTitleSize(0.045); + fh_mult_fibAll->GetYaxis()->SetLabelSize(0.045); + fh_mult_fibAll->GetYaxis()->SetTitleSize(0.045); + + fh_ToF_vs_charge = new TH2F("tof_vs_q", "tof vs q",100, 0.,10., 2000,-200.,0.); + fh_ToF_vs_charge->GetXaxis()->SetTitle("Charge"); + fh_ToF_vs_charge->GetYaxis()->SetTitle("ToF ns"); + fh_ToF_vs_charge->GetYaxis()->SetTitleOffset(1.); + fh_ToF_vs_charge->GetXaxis()->CenterTitle(true); + fh_ToF_vs_charge->GetYaxis()->CenterTitle(true); + fh_ToF_vs_charge->GetXaxis()->SetLabelSize(0.045); + fh_ToF_vs_charge->GetXaxis()->SetTitleSize(0.045); + fh_ToF_vs_charge->GetYaxis()->SetLabelSize(0.045); + fh_ToF_vs_charge->GetYaxis()->SetTitleSize(0.045); + + cMult->Divide(3,3); + cMult->cd(1); + gPad->SetLogy(); + fh_mult_coinc->Draw(); + cMult->cd(2); + gPad->SetLogy(); + fh_mult_minus_los->Draw(); + cMult->cd(3); + gPad->SetLogy(); + fh_mult_fib30->Draw(); + cMult->cd(4); + gPad->SetLogy(); + fh_mult_fib31->Draw(); + cMult->cd(5); + gPad->SetLogy(); + fh_mult_fib32->Draw(); + cMult->cd(6); + gPad->SetLogy(); + fh_mult_fib33->Draw(); + cMult->cd(7); + gPad->SetLogy(); + fh_mult_fibAll->Draw(); + cMult->cd(8); + gPad->SetLogz(); + fh_ToF_vs_charge->Draw("colz"); + + + + // MAIN FOLDER + TFolder* mainfol = new TFolder("Tracking_Fibers", "Tracking info"); + mainfol->Add(cTrackingXZ); + mainfol->Add(cAngVsX); + mainfol->Add(cslopeVsX); + mainfol->Add(AngCorrFibToFD); + mainfol->Add(cPID_3D); + mainfol->Add(cToFD_Z); + mainfol->Add(cAngRot); + mainfol->Add(cSlopeRot); + mainfol->Add(cFibToT_Q); + mainfol->Add(cFib30Rel); + mainfol->Add(cFibToTCorr); + mainfol->Add(cFibToFDnocuts); + mainfol->Add(cFibInfo); + mainfol->Add(cXcorr); + mainfol->Add(cDeltaQ); + mainfol->Add(cMult); + + FairRunOnline* run = FairRunOnline::Instance(); + R3BLOG_IF(fatal, NULL == run, "FairRunOnline not found"); + run->GetHttpServer()->Register("", this); + run->AddObject(mainfol); + + // Register command to reset histograms + run->GetHttpServer()->RegisterCommand("Reset_Fiber_Tracking_HIST", Form("/Objects/%s/->Reset_Histo()", GetName())); + + return kSUCCESS; +} + +void R3BFiberTrackingOnlineSpectra::Reset_Histo() +{ + R3BLOG(info, ""); + fh2_fibtracking_planeXZ->Reset(); + fh2_fib_ang_vs_x->Reset(); + fh2_fib_ang_vs_x4->Reset(); + fh2_fib_ang_vs_x5->Reset(); + fh2_fib_ang_vs_x6->Reset(); + fh2_fib_slope_vs_x->Reset(); + fh2_fib_slope_vs_x4->Reset(); + fh2_fib_slope_vs_x5->Reset(); + fh2_fib_slope_vs_x6->Reset(); + // fh2_ZvsBeta->Reset(); + fh2_FibToT_vs_QToFD->Reset(); + fh_ang_vs_pos_rot->Reset(); + fh_slope_vs_pos_rot->Reset(); + fh_pid->Reset(); + fh_pid_slope->Reset(); + fh2_ang_Fib_vs_ToFD->Reset(); + fh2_ang_Fib_vs_ToFD_nocuts->Reset(); + fh_ToFD_Charge->Reset(); + fh_pid_3d->Reset(); +} + +void R3BFiberTrackingOnlineSpectra::Exec(Option_t* option) +{ + + fNEvents += 1; + Bool_t ToFDgate = true; + + fh_mult_fib30->Fill(-1); + fh_mult_fib31->Fill(-1); + fh_mult_fib32->Fill(-1); + fh_mult_fib33->Fill(-1); + fh_mult_fibAll->Fill(-1); + + if ((fTrigger >= 0) && header && (header->GetTrigger() != fTrigger)) return; + bool onspill = true; + //fTpat = 1-16; fTpat_bit = 0-15 + fTpat1 = 16; fTpat2 = 13; + if (fTpat1 > -1 && fTpat2 > -1) + { + Int_t fTpat_bit1 = fTpat1 - 1; + Int_t fTpat_bit2 = fTpat2 - 1; + Int_t tpatbin; + for (int i = 0; i < 16; i++) + { + tpatbin = (header->GetTpat() & (1 << i)); + if (tpatbin != 0 && (i < fTpat_bit1 && i > fTpat_bit2)) + { + onspill = false; + } + } + } + + //cout<<"Tpat "<GetTpat()<<" mask "<GetEntriesFast()<GetTpat() & fTPatMask){ +if(onspill){ + fh_mult_fib30->Fill(0); + fh_mult_fib31->Fill(0); + fh_mult_fib32->Fill(0); + fh_mult_fib33->Fill(0); + fh_mult_fibAll->Fill(0); + if(fFi31HitDataCA && fFi33HitDataCA){ + if( + fFi32HitDataCA->GetEntriesFast() == fFi31HitDataCA->GetEntriesFast() + fFi33HitDataCA->GetEntriesFast() && 4 * fFi32HitDataCA->GetEntriesFast() == fHitTofdItems->GetEntriesFast() + ){ + fh_mult_coinc->Fill(fFi32HitDataCA->GetEntriesFast()); + fh_mult_minus_los->Fill(fCalLos->GetEntriesFast() - fFi32HitDataCA->GetEntriesFast()); + // if(fCalLos->GetEntriesFast() - fFi32HitDataCA->GetEntriesFast() < 0) cout<<"Los "<GetEntriesFast()<<" fib "<GetEntriesFast()<Fill(-1); + } +} + if((header && !(header->GetTpat() & fTPatMask))) return; + if((/*fHitFrs->GetEntriesFast() != fFi32HitDataCA->GetEntriesFast() ||*/fFi32HitDataCA->GetEntriesFast() != 1 ) || (4.*(fFi32HitDataCA->GetEntriesFast()) != fHitTofdItems->GetEntriesFast())) return; + if(!(fFi32HitDataCA->GetEntriesFast() == 1)) return; + if(!(4.*(fFi32HitDataCA->GetEntriesFast()) == fHitTofdItems->GetEntriesFast())) return; + if(fHitTofdItems->GetEntriesFast() != 4) return; + + Double_t xS2=0.,Aq=0.,ToF=0.,Beta=0.,Brho=0.,Q_FRS=0.; + if(fFRSGATE){ + + //FRS HIT DATA + if (fHitFrs && fHitFrs->GetEntriesFast() > 0) + { + //cout<<"Hit FRS"<GetEntriesFast(); + //cout<<"Hit FRS "<At(ihit); + if (!hit) + continue; + //if (hit->GetStaId() != fStaId) + if (hit->GetStaId() != 1) + continue; + xS2 = hit->GetXS2(); + Aq = hit->GetAq(); + ToF = hit->GetTof(); + Beta = hit->GetBeta(); + Brho = hit->GetBrho(); + Q_FRS = hit->GetZ(); + //cout<<"i: "<GetEntriesFast()<GetEntriesFast()!=1) return; + +#if 0 +////Start LOS + Int_t fNofLosDetectors = 1; + Double_t timeTofd = 0; + + Double_t timeLosV[fNofLosDetectors][32]; + Double_t LosTresV[fNofLosDetectors][32]; + Double_t timeLosT[fNofLosDetectors][32]; + Double_t LosTresT[fNofLosDetectors][32]; + Double_t timeLos[fNofLosDetectors][32]; + Double_t totsum[fNofLosDetectors][32]; + Double_t xT_cm[fNofLosDetectors][32]; + Double_t yT_cm[fNofLosDetectors][32]; + Double_t xToT_cm[fNofLosDetectors][32]; + Double_t yToT_cm[fNofLosDetectors][32]; + Double_t xV_cm[fNofLosDetectors][32]; + Double_t yV_cm[fNofLosDetectors][32]; + + Double_t time_V[fNofLosDetectors][32][8]; // [det][multihit][pm] + Double_t time_L[fNofLosDetectors][32][8]; + Double_t time_T[fNofLosDetectors][32][8]; + Double_t tot[fNofLosDetectors][32][8]; + Double_t time_MTDC[32][8] = { 0. }; + Double_t LosTresMTDC[32]; + + for (Int_t idet = 0; idet < fNofLosDetectors; idet++) + { + for (Int_t imult = 0; imult < 32; imult++) + { + timeLosV[idet][imult] = 0.0; + LosTresV[idet][imult] = 0.0 / 0.0; + timeLosT[idet][imult] = 0.0; + LosTresT[idet][imult] = 0.0 / 0.0; + timeLos[idet][imult] = 0.0; + totsum[idet][imult] = 0.0; + xT_cm[idet][imult] = 0.0 / 0.0; + yT_cm[idet][imult] = 0.0 / 0.0; + xToT_cm[idet][imult] = -100000.; + yToT_cm[idet][imult] = -100000.; + xV_cm[idet][imult] = 0.0 / 0.0; + yV_cm[idet][imult] = 0.0 / 0.0; + for (Int_t icha = 0; icha < 8; icha++) + { + time_V[idet][imult][icha] = 0.0 / 0.0; // [det][multihit][pm] + time_L[idet][imult][icha] = 0.0 / 0.0; + time_T[idet][imult][icha] = 0.0 / 0.0; + tot[idet][imult][icha] = 0.0 / 0.0; + } + } + } + Int_t nPartLOS = 0; + Int_t nPartc[fNofLosDetectors]; + for (Int_t d = 0; d < fNofLosDetectors; d++) nPartc[d] = 0; + + Bool_t iLOSType[fNofLosDetectors][32]; + Bool_t iLOSPileUp[fNofLosDetectors][32]; + for (Int_t idet = 0; idet < fNofLosDetectors; idet++) + { + for (Int_t imult = 0; imult < 32; imult++) + { + iLOSType[idet][imult] = false; + iLOSPileUp[idet][imult] = false; + } + } + + nPartLOS = fCalLos->GetEntriesFast(); + + Int_t iDet = 0; + Double_t time_V_LOS1[32][8] = { 0. }; + Double_t time_V_LOS2[32][8] = { 0. }; + + for (Int_t iPart = 0; iPart < nPartLOS; iPart++) + { + /* + * nPart is the number of particle passing through LOS detector in one event + */ + R3BLosCalData* calData = (R3BLosCalData*)fCalLos->At(iPart); + iDet = calData->GetDetector(); + + Double_t sumvtemp = 0, sumltemp = 0, sumttemp = 0; + for (Int_t iCha = 0; iCha < 8; iCha++) + { + sumvtemp += calData->GetTimeV_ns(iCha); + sumltemp += calData->GetTimeL_ns(iCha); + sumttemp += calData->GetTimeT_ns(iCha); + } + if (!(IS_NAN(sumvtemp)) && !(IS_NAN(sumltemp)) && !(IS_NAN(sumltemp))) + { + nPartc[iDet - 1]++; + } + else + { + continue; + } + + for (Int_t iCha = 0; iCha < 8; iCha++) + { + if (iDet == 1) + { + time_V_LOS1[nPartc[iDet - 1] - 1][iCha] = 0. / 0.; + if (!(IS_NAN(calData->GetTimeV_ns(iCha)))) + { // VFTX + time_V_LOS1[nPartc[iDet - 1] - 1][iCha] = calData->GetTimeV_ns(iCha); + } + if (!(IS_NAN(calData->GetTimeL_ns(iCha)))) + { // TAMEX leading + time_L[iDet - 1][nPartc[iDet - 1] - 1][iCha] = calData->GetTimeL_ns(iCha); + } + + if (!(IS_NAN(calData->GetTimeT_ns(iCha)))) + { // TAMEX trailing + time_T[iDet - 1][nPartc[iDet - 1] - 1][iCha] = calData->GetTimeT_ns(iCha); + } + time_MTDC[nPartc[iDet - 1]][iCha] = 0. / 0.; + if (!(IS_NAN(calData->GetTimeM_ns(iCha)))) + { // MTDC + time_MTDC[nPartc[iDet - 1] - 1][iCha] = calData->GetTimeM_ns(iCha); + } + } + if (iDet == 2) + { + time_V_LOS2[nPartc[iDet - 1] - 1][iCha] = 0. / 0.; + if (!(IS_NAN(calData->GetTimeV_ns(iCha)))) + { // VFTX + time_V_LOS2[nPartc[iDet - 1] - 1][iCha] = calData->GetTimeV_ns(iCha); + } + if (!(IS_NAN(calData->GetTimeL_ns(iCha)))) + { // TAMEX leading + time_L[iDet - 1][nPartc[iDet - 1] - 1][iCha] = calData->GetTimeL_ns(iCha); + } + if (!(IS_NAN(calData->GetTimeT_ns(iCha)))) + { // TAMEX trailing + time_T[iDet - 1][nPartc[iDet - 1] - 1][iCha] = calData->GetTimeT_ns(iCha); + } + } + } + + if (!calData) + { + cout << "LOS !calData" << endl; + continue; // can this happen? + } + } + + // Sorting VFTX data: + + // detector 1 + if (nPartc[0] > 0) + { + std::qsort(time_V_LOS1, nPartc[0], sizeof(*time_V_LOS1), [](const void* arg1, const void* arg2) -> int { + double const* lhs = static_cast(arg1); + double const* rhs = static_cast(arg2); + + return (lhs[0] < rhs[0]) ? -1 : ((rhs[0] < lhs[0]) ? 1 : 0); + }); + for (Int_t iPart = 0; iPart < nPartc[0]; iPart++) + { + for (int ipm = 0; ipm < 8; ipm++) + { + time_V[0][iPart][ipm] = time_V_LOS1[iPart][ipm]; + } + } + } + else if(nPartc[0] > 1) return; + + // detector 2 + if (fNofLosDetectors > 1 && nPartc[1] > 0) + { + std::qsort(time_V_LOS2, nPartc[1], sizeof(*time_V_LOS2), [](const void* arg1, const void* arg2) -> int { + double const* lhs = static_cast(arg1); + double const* rhs = static_cast(arg2); + + return (lhs[0] < rhs[0]) ? -1 : ((rhs[0] < lhs[0]) ? 1 : 0); + }); + for (Int_t iPart = 0; iPart < nPartc[1]; iPart++) + { + for (int ipm = 0; ipm < 8; ipm++) + { + time_V[1][iPart][ipm] = time_V_LOS2[iPart][ipm]; + } + } + } + + // End sorting + + std::vector time_first(fNofLosDetectors, -1.); + std::vector time0(fNofLosDetectors, -1.); + std::vector time1(fNofLosDetectors, -1.); + std::vector time_abs(fNofLosDetectors, -1.); + + for (iDet = 1; iDet <= fNofLosDetectors; iDet++) + { + for (Int_t iPart = 0; iPart < nPartc[iDet - 1]; iPart++) + { + Bool_t iLOSTypeMCFD = false; + Bool_t iLOSTypeTAMEX = false; + + if (time_V[iDet - 1][iPart][0] > 0. && !(IS_NAN(time_V[iDet - 1][iPart][0])) && + time_V[iDet - 1][iPart][1] > 0. && !(IS_NAN(time_V[iDet - 1][iPart][1])) && + time_V[iDet - 1][iPart][2] > 0. && !(IS_NAN(time_V[iDet - 1][iPart][2])) && + time_V[iDet - 1][iPart][3] > 0. && !(IS_NAN(time_V[iDet - 1][iPart][3])) && + time_V[iDet - 1][iPart][4] > 0. && !(IS_NAN(time_V[iDet - 1][iPart][4])) && + time_V[iDet - 1][iPart][5] > 0. && !(IS_NAN(time_V[iDet - 1][iPart][5])) && + time_V[iDet - 1][iPart][6] > 0. && !(IS_NAN(time_V[iDet - 1][iPart][6])) && + time_V[iDet - 1][iPart][7] > 0. && !(IS_NAN(time_V[iDet - 1][iPart][7]))) + { + iLOSTypeMCFD = true; // all 8 MCFD times + } + + if (time_L[iDet - 1][iPart][0] > 0. && !(IS_NAN(time_L[iDet - 1][iPart][0])) && + time_L[iDet - 1][iPart][1] > 0. && !(IS_NAN(time_L[iDet - 1][iPart][1])) && + time_L[iDet - 1][iPart][2] > 0. && !(IS_NAN(time_L[iDet - 1][iPart][2])) && + time_L[iDet - 1][iPart][3] > 0. && !(IS_NAN(time_L[iDet - 1][iPart][3])) && + time_L[iDet - 1][iPart][4] > 0. && !(IS_NAN(time_L[iDet - 1][iPart][4])) && + time_L[iDet - 1][iPart][5] > 0. && !(IS_NAN(time_L[iDet - 1][iPart][5])) && + time_L[iDet - 1][iPart][6] > 0. && !(IS_NAN(time_L[iDet - 1][iPart][6])) && + time_L[iDet - 1][iPart][7] > 0. && !(IS_NAN(time_L[iDet - 1][iPart][7])) && + + time_T[iDet - 1][iPart][0] > 0. && !(IS_NAN(time_T[iDet - 1][iPart][0])) && + time_T[iDet - 1][iPart][1] > 0. && !(IS_NAN(time_T[iDet - 1][iPart][1])) && + time_T[iDet - 1][iPart][2] > 0. && !(IS_NAN(time_T[iDet - 1][iPart][2])) && + time_T[iDet - 1][iPart][3] > 0. && !(IS_NAN(time_T[iDet - 1][iPart][3])) && + time_T[iDet - 1][iPart][4] > 0. && !(IS_NAN(time_T[iDet - 1][iPart][4])) && + time_T[iDet - 1][iPart][5] > 0. && !(IS_NAN(time_T[iDet - 1][iPart][5])) && + time_T[iDet - 1][iPart][6] > 0. && !(IS_NAN(time_T[iDet - 1][iPart][6])) && + time_T[iDet - 1][iPart][7] > 0. && !(IS_NAN(time_T[iDet - 1][iPart][7]))) + { + iLOSTypeTAMEX = true; // all 8 leading and trailing times + } + + // We will consider only events in which booth MCFD and TAMEX see same number of channels: + if (iLOSTypeTAMEX && iLOSTypeMCFD) + iLOSType[iDet - 1][iPart] = true; + + if (iLOSType[iDet - 1][iPart]) + { + int nPMT = 0; + int nPMV = 0; + + for (int ipm = 0; ipm < 8; ipm++) + { + + if (time_T[iDet - 1][iPart][ipm] > 0. && time_L[iDet - 1][iPart][ipm] > 0. && + !(IS_NAN(time_T[iDet - 1][iPart][ipm])) && !(IS_NAN(time_L[iDet - 1][iPart][ipm]))) + { + while (time_T[iDet - 1][iPart][ipm] - time_L[iDet - 1][iPart][ipm] <= 0.) + { + time_T[iDet - 1][iPart][ipm] = time_T[iDet - 1][iPart][ipm] + 2048. * 5./*fClockFreq*/; + } + + nPMT = nPMT + 1; + tot[iDet - 1][iPart][ipm] = time_T[iDet - 1][iPart][ipm] - time_L[iDet - 1][iPart][ipm]; + + // pileup rejection + // if (tot[iDet - 1][iPart][ipm] > fEpileup) + // iLOSPileUp[iDet - 1][iPart] = true; + } + + if (tot[iDet - 1][iPart][ipm] != 0. && !(IS_NAN(tot[iDet - 1][iPart][ipm]))) + totsum[iDet - 1][iPart] += tot[iDet - 1][iPart][ipm]; + + if (time_L[iDet - 1][iPart][ipm] > 0. && !(IS_NAN(time_L[iDet - 1][iPart][ipm]))) + timeLosT[iDet - 1][iPart] += time_L[iDet - 1][iPart][ipm]; + + // Calculate detector time + if (time_V[iDet - 1][iPart][ipm] > 0. && !(IS_NAN(time_V[iDet - 1][iPart][ipm]))) + { + timeLosV[iDet - 1][iPart] += time_V[iDet - 1][iPart][ipm]; + nPMV = nPMV + 1; + } + } + + totsum[iDet - 1][iPart] = totsum[iDet - 1][iPart] / nPMT; + + timeLosV[iDet - 1][iPart] = timeLosV[iDet - 1][iPart] / nPMV; + + timeLosT[iDet - 1][iPart] = timeLosT[iDet - 1][iPart] / nPMT; + + timeLos[iDet - 1][iPart] = timeLosV[iDet - 1][iPart]; + + LosTresV[iDet - 1][iPart] = ((time_V[iDet - 1][iPart][0] + time_V[iDet - 1][iPart][2] + + time_V[iDet - 1][iPart][4] + time_V[iDet - 1][iPart][6]) - + (time_V[iDet - 1][iPart][1] + time_V[iDet - 1][iPart][3] + + time_V[iDet - 1][iPart][5] + time_V[iDet - 1][iPart][7])) / + 4.; + + LosTresT[iDet - 1][iPart] = ((time_L[iDet - 1][iPart][0] + time_L[iDet - 1][iPart][2] + + time_L[iDet - 1][iPart][4] + time_L[iDet - 1][iPart][6]) - + (time_L[iDet - 1][iPart][1] + time_L[iDet - 1][iPart][3] + + time_L[iDet - 1][iPart][5] + time_L[iDet - 1][iPart][7])) / + 4.; + + // right koord.-system, Z-axis beam direction: + // Position from tamex: + // xT_cm[iDet - 1][iPart] = (time_L[iDet - 1][iPart][1] + time_L[iDet - 1][iPart][2]) / 2. - + // (time_L[iDet - 1][iPart][5] + time_L[iDet - 1][iPart][6]) / 2.; + // yT_cm[iDet - 1][iPart] = (time_L[iDet - 1][iPart][3] + time_L[iDet - 1][iPart][4]) / 2. - + // (time_L[iDet - 1][iPart][7] + time_L[iDet - 1][iPart][0]) / 2.; + // xT_cm[iDet - 1][iPart] = (xT_cm[iDet - 1][iPart] - flosOffsetXT[iDet - 1]) * flosVeffXT[iDet - 1]; + // yT_cm[iDet - 1][iPart] = (yT_cm[iDet - 1][iPart] - flosOffsetYT[iDet - 1]) * flosVeffYT[iDet - 1]; +// + // // Position from VFTX: + // xV_cm[iDet - 1][iPart] = (time_V[iDet - 1][iPart][1] + time_V[iDet - 1][iPart][2]) / 2. - + // (time_V[iDet - 1][iPart][5] + time_V[iDet - 1][iPart][6]) / 2.; + // yV_cm[iDet - 1][iPart] = (time_V[iDet - 1][iPart][3] + time_V[iDet - 1][iPart][4]) / 2. - + // (time_V[iDet - 1][iPart][7] + time_V[iDet - 1][iPart][0]) / 2.; + // xV_cm[iDet - 1][iPart] = (xV_cm[iDet - 1][iPart] - flosOffsetXV[iDet - 1]) * flosVeffXV[iDet - 1]; + // yV_cm[iDet - 1][iPart] = (yV_cm[iDet - 1][iPart] - flosOffsetYV[iDet - 1]) * flosVeffYV[iDet - 1]; +// + // // Position from ToT: + // if (tot[iDet - 1][iPart][1] > 0. && tot[iDet - 1][iPart][2] > 0. && tot[iDet - 1][iPart][5] > 0. && + // tot[iDet - 1][iPart][6] > 0. && tot[iDet - 1][iPart][0] > 0. && tot[iDet - 1][iPart][3] > 0. && + // tot[iDet - 1][iPart][4] > 0. && tot[iDet - 1][iPart][7] > 0.) + // { + // xToT_cm[iDet - 1][iPart] = (((tot[iDet - 1][iPart][5] + tot[iDet - 1][iPart][6]) / 2. - + // (tot[iDet - 1][iPart][1] + tot[iDet - 1][iPart][2]) / 2.) / + // ((tot[iDet - 1][iPart][1] + tot[iDet - 1][iPart][2] + + // tot[iDet - 1][iPart][5] + tot[iDet - 1][iPart][6]) / + // 4.)); + + // yToT_cm[iDet - 1][iPart] = (((tot[iDet - 1][iPart][0] + tot[iDet - 1][iPart][7]) / 2. - + // (tot[iDet - 1][iPart][3] + tot[iDet - 1][iPart][4]) / 2.) / + // ((tot[iDet - 1][iPart][7] + tot[iDet - 1][iPart][0] + + // tot[iDet - 1][iPart][3] + tot[iDet - 1][iPart][4]) / + // 4.)); + // } +// + // xToT_cm[iDet - 1][iPart] = + // (xToT_cm[iDet - 1][iPart] - flosOffsetXQ[iDet - 1]) * flosVeffXQ[iDet - 1]; + // yToT_cm[iDet - 1][iPart] = + // (yToT_cm[iDet - 1][iPart] - flosOffsetYQ[iDet - 1]) * flosVeffYQ[iDet - 1]; +// + // if (timeLosV[iDet - 1][iPart] > 0. && timeLosV[iDet - 1][iPart] < 8192. * 5. && + // !(IS_NAN(timeLosV[iDet - 1][iPart]))) + // { + // while (timeLosT[iDet - 1][iPart] - timeLosV[iDet - 1][iPart] < 2048. * 5. / 2.) + // { + // timeLosT[iDet - 1][iPart] = timeLosT[iDet - 1][iPart] + 2048. * 5.; + // } + // while (timeLosT[iDet - 1][iPart] - timeLosV[iDet - 1][iPart] > 2048. * 5. / 2.) + // { + // timeLosT[iDet - 1][iPart] = timeLosT[iDet - 1][iPart] - 2048. * 5.; + // } +// + // + // } + //cout<<"los tot "< 70.) return;//C16 Setting + if(totsum[iDet - 1][iPart] < 68. || totsum[iDet - 1][iPart] > 73.6) return;//C12 Setting + + } // if iLosType + } // for iPart + } // for iDet + + // if fCallItems + +////End LOS +#endif + + + + + +///B13 13.9°, C16 13.6°, + if(!fFRSGATE || + (Aq < fAQmax && Aq > fAQmin && Q_FRS > fQmin && Q_FRS < fQmax) + ) + { +//cout<<"FOUND ISOTOPE OF INTEREST"<GetEntriesFast() > 0 && + (fFi31HitDataCA->GetEntriesFast() > 0 || + fFi33HitDataCA->GetEntriesFast() > 0)) + )//First x det (Fi32) needs to be hit. Second row (Fi31 or Fi33) needs at least in 1 det a hit. + { + Double_t xpos32 = 0./0., xpos31 = 0./0., xpos33 = 0./0., xposback = 0./0., xpos30=0./0.; + Double_t ypos32 = 0./0., ypos31 = 0./0., ypos33 = 0./0., ypos30=0./0.; + Double_t tot32 = 0., tot31 = 0., tot33 = 0., tot30 = 0., t32 = 0.; + Int_t nHits32 = fFi32HitDataCA->GetEntriesFast(); + Int_t Fib32 = 0; + + //cout<<"nHits32: "<1) return; + for(Int_t ihit=0; ihitAt(ihit); + if(!hit) + continue; + + //cout<<"ihit: "<GetEloss()<<" fID: "<GetFiberId()<GetEloss() > tot32) + { + tot32 = hit->GetEloss(); + xpos32 = hit->GetX() /*cm*/; + //cout<<"ihit: "<<<GetX()<GetY(); + Fib32 = hit->GetFiberId(); + t32 = hit->GetTime(); + } + } + + Int_t nHits31 = fFi31HitDataCA->GetEntriesFast(); + Int_t nHits33 = fFi33HitDataCA->GetEntriesFast(); + Int_t nHits30 = fFi30HitDataCA->GetEntriesFast(); + + if( !( (nHits31 == 1 && nHits33 == 0) || (nHits31 == 0 && nHits33 == 1) ) ) return; + + //cout<<"before if: n31 "<At(ihit); + if(!hit) + continue; + if (hit->GetEloss() > tot30) + { + tot30 = hit->GetEloss(); + xpos30 = hit->GetX() - 25.6 /*cm*/; + ypos30= hit->GetY(); + } + } + } + + for(Int_t ihit=0; ihitAt(ihit); + if(!hit) + continue; + if (hit->GetEloss() > tot31) + { + tot31 = hit->GetEloss(); + xpos31 = hit->GetX() - 25.6 /*cm*/; + ypos31= hit->GetY(); + } + } + + for(Int_t ihit=0; ihitAt(ihit); + if(!hit) + continue; + if (hit->GetEloss() > tot33) + { + tot33 = hit->GetEloss(); + xpos33 = 25. + hit->GetX() /*cm*/; + ypos33= hit->GetY(); + } + } + if(overlap){ + if(xpos33 > 0.3 || xpos31 < -0.3) return; + } + + if(tot33 > tot31){ + xposback = xpos33; + LeftHit = true; + } + else if(tot33 <= tot31){ + xposback = xpos31; + RightHit = true; + } + + if(!IS_NAN(xposback)){ + + Double_t angX = 0./0.; + Double_t totback = 0.; + if(LeftHit){ + angX = -100.*(xpos32 - xposback) / (160.4) /*cm*/; + totback = tot33; + } + else if(RightHit){ + angX = -100.*(xpos32 - xposback) / (177.3) /*cm*/; + totback = tot31; + } + + + + + Double_t tof = 0./0.; + Double_t tofdq[fNbTofdPlanes]; + Double_t tofdtof[fNbTofdPlanes]; + Int_t bar[fNbTofdPlanes]; + if(!IS_NAN(angX) || 1){ + // fh2_fib_ang_vs_x->Fill(xpos32,angX); + Double_t xpos2[fNbTofdPlanes]; + Double_t q = 0./0.; + Double_t q1= 0./0.; + Double_t xmean = 0.; + if(ToFDgate && fHitTofdItems && fHitTofdItems->GetEntriesFast()>0){ + + Double_t xpos1 = 0. / 0.; + Double_t ypos1 = 0. / 0., ypos2[fNbTofdPlanes]; + //Double_t tot = 0.; + + + for (Int_t i = 0; i < fNbTofdPlanes; i++) + { + xpos2[i] = 0. / 0.; + ypos2[i] = 0. / 0.; + tofdq[i] = 0.; + bar[i] = -1; + } + + Int_t nHits2 = fHitTofdItems->GetEntriesFast(); + //cout<<"nHitsToFD: "<At(ihit); + if (!hitToFD) + continue; + // Looking for the maximum + Int_t iPlane = hitToFD->GetDetId() - 1; + if(iPlane > 1) continue; + if(iPlane==0){ + cout<<"ToFDQ: "<GetEloss()<<" Ch: "<GetBarId()<GetEloss(); + } + if (hitToFD->GetEloss() > tofdq[iPlane]) + { + tofdq[iPlane] = hitToFD->GetEloss(); + xpos2[iPlane] = hitToFD->GetX(); + ypos2[iPlane] = hitToFD->GetY(); + bar[iPlane] = hitToFD->GetBarId(); + tofdtof[iPlane] = hitToFD->GetTof(); + } + } + tof = tofdtof[0]; + + if(bar[3] == 5 || bar[3] == 8){ + //cout<<"bar "< -40 && xpos2[0] < -38.8){ + if(tofdq[0]> 5. || tofdq[1]> 5. || tofdq[2]> 5. ||tofdq[3]> 5.) for(int x = 0; x<4; x++) cout<<"x "< 0. ) return; + + Int_t SetPlane = 2; + //if(q<2. || 1){ + for(int i = 0; i 2)){ //for now, accept +-2 bars across planes + straggler = true; + } + } + + int bararray[4] = {0,0,0,0}; + int badbar_ctr = 0; + if(straggler){ + /// cout<<"straggler"< tofdq[b+1]) { + tmp=tofdq[b+1]; + tofdq[b+1]=tofdq[b]; + tofdq[b]=tmp; + tmpb = bar[b+1]; + bar[b+1] = bar[b]; + bar[b] = tmpb; + tmpt = tofdtof[b+1]; + tofdtof[b+1] = tofdtof[b]; + tofdtof[b] = tmpt; + } + } + } + } + /// for(int i = 0; i<4;i++) cout<<"tofdq["<0.04 && !IS_NAN(q)); + + //for(int i=0; i<4;i++) tofdq[i] = (double)tofdq[i]; + /// cout<<"badcharge_ctr "< -12.5) + + + q = q1; + cout<<"Mean ToFD charge: "<3) q=0.; + //if(q>10.){ //exit(0); + // cout<<"q = "<5.93 && q<6.18 && /*tot32 < 14. &&*/ nHits32 == 1){ + fh_ToT_Fib->Fill(Fib32,tot32); + //cout<<"Fib mult "<0.) + // fh_ToFD_Charge->Fill(q); + + + if(q>4.46 && q<4.94){ + int gh = 0; + for(int l= 0; l<4; l++){ + for(int p = l+1; p<4; p++){ + fh_tofd_charge_planes[gh]->Fill(barV[l], tofdqV[l] - tofdqV[p]); + gh++; + } + } + } + + + } + + if(q>4.){ + if(nHits30 > 0 && tot30 > 14.) fh_mult_fib30->Fill(nHits30); + if(nHits31 > 0 && tot31 > 14.) fh_mult_fib31->Fill(nHits31); + if(nHits32 > 0 && tot32 > 14.) fh_mult_fib32->Fill(nHits32); + if(nHits33 > 0 && tot33 > 14.) fh_mult_fib33->Fill(nHits33); + if(nHits30 > 0 && tot30 > 14. && + nHits32 > 0 && tot32 > 14. && + ((nHits31 > 0 && tot31 > 14.) || + (nHits33 > 0 && tot33 > 14.))) fh_mult_fibAll->Fill(1); + } + + + double xdiff = 0./0.; + double yy = 0./0.; + if(RightHit){ + yy = 611.4; + xdiff = 527.1; + LeftHit = false; + } + if(LeftHit){ + yy = 594.5; + xdiff = 544.; + RightHit = false; + } + + Double_t angFibToFD = ((100*((xposback - xmean) / xdiff)) + (100*((xpos32 - xmean) / 704.4))) / 2.; + Double_t angMean = (-angX + angFibToFD)/2.; + + + double xsum = 0., ysum = 0., x2sum = 0., xysum = 0.,y2sum=0.,r=0.; + double xpoints[3] = {434.1, yy, 1138.5}; + double ypoints[3] = {xpos32, xposback, xmean}; + //double ypoints = 0.; + + for (int i=0;i<3;i++) + { + //if(0==i) ypoints = xpos32; + //if(1==i) ypoints = xposback; + //if(2==i) ypoints = xmean; + + xsum=xsum+xpoints[i]; //calculate sigma(xi) + ysum=ysum+ypoints[i]; //calculate sigma(yi) + x2sum=x2sum+pow(xpoints[i],2); //calculate sigma(x^2i) + y2sum=y2sum+pow(ypoints[i],2); + xysum=xysum+xpoints[i]*ypoints[i]; //calculate sigma(xi*yi) + } + double slope =(2.*xysum-xsum*ysum)/(2.*x2sum-xsum*xsum); //calculate slope + double intercept =(x2sum*ysum-xsum*xysum)/(x2sum*2.-xsum*xsum); //calculate intercept + + //cout<<"linfit: slope = "< 0.993){ + //cout<<"x32 "<GetTpat() & (1 << i)); + if (tpatbin != 0 && (i < 16 || i > -1)) + { + //if(i+1 == 3 || i+1 == 4 || i+1 == 5 || i+1 == 9 || i+1 == 10 || i+1 == 11) return; + if(i+1 == 1 || i+1 == 2 || i+1 == 5 || i+1 == 6 || i+1 == 7 || i+1 == 8 || i+1 == 11 || i+1==12) return; + fh_pos_TPat->Fill(xmean, i+1); + } + } + #endif + + fh_x_vs_tof->Fill(xposback,tofdtof[0]); + fh_x_vs_q->Fill(xposback,q); + + fh_X_fib_vs_tofd_no_cut->Fill(xmean, intercept+slope*1138.5); + fh_X_fib_vs_tofd_ang_no_cut->Fill(xmean, xposback - (angX) * 602.95 /100. + (angX) * 1138.5 /100.); + + /*x1=-3.3 y1=-2.8 | x2=5.1 y2=5.8 + s = (5.8 + 2.8)/(5.1 + 3.3) = 1.0238095 + b= 0.57857155 + x1=-34, y1=-2 | x2=19.40 ,y2=0.7 + s=0.050561798 + b= -0.28089888 + x1=-33.8 , y1=5.2 | x2=29.2 ,y2=-4.6 + s = -0.15555556 + b = -0.057777928 + */ + + //For better PID: Rotate angle vs pos plot via TVector3, project that on "new" x-axis, then plot against Charge ID + //ToDo: Get rotation angle more dynamically. Now precalculated fixed angle. +// TVector3 v1(xposback, angMean -0.0904221,0.); +// //TVector3 v2(intercept, 100*slope + 1.40554,0.); +// //TVector3 v2(xposback, -100*atan(slope) + 1.40554,0.); +// TVector3 v2(xposback, -100*atan(slope) + 0.057777928,0.); +// //double rotAngle2 = -0.079104685; //radiants +// double rotAngle2 = atan(-0.15555556); +// double rotAngle1 = -0.18909835 + 0.050518777; //radiants +// // double rotAngle = -10.834537; //degree +// v1.RotateZ(-rotAngle1); +// v2.RotateZ(-rotAngle2); + if(q>qmin && q5.6 && q<6.5*/){ + fh2_ang_Fib_vs_ToFD_nocuts->Fill(-angX,angFibToFD); //Fill before applying any cuts + fh2_fib_slope_vs_x_nocuts->Fill(xposback, -1000.*(atan(slope)/*+0.2443461014° in rad*/)); + // fh2_fib_ang_vs_x_nocuts->Fill(xposback,angMean); + } + for (int i=0; i<100;i++){ + Double_t zrand = gRandom->Uniform(0., 1230.5/*cm*/);//Umlenkpunkt at 0. Fib32 at 434.1 cm, ToFD at 1138.5cm + //fh2_fibtracking_planeXZ->Fill(zrand, xpos32 - angX * 434.1 /1000. + angX * zrand /1000.); + //if(q>5.93 && q<6.18/*q>5.6 && q<6.5*/) fh2_fibtracking_planeXZ_nocuts->Fill(zrand, intercept + slope * zrand ); + } + + + if((q>4.6 && q<6.5)|| !ToFDgate){ + fh2_fib_ang_vs_x->Fill(xposback,-angX); + fh2_fib_slope_vs_x->Fill(xposback, -1000.*(atan(slope)/*+0.2443461014° in rad*/)); + if(q>=3.8 && q < 4.7){ + fh2_fib_ang_vs_x4->Fill(xposback,angMean); + fh2_fib_slope_vs_x4->Fill(xposback, -1000.*(atan(slope)/*+0.2443461014° in rad*/)); + } + if(q>=4.6 && q < 5.6){ + fh2_fib_ang_vs_x5->Fill(xposback,angMean); + fh2_fib_slope_vs_x5->Fill(xposback, -1000.*(atan(slope)/*+0.2443461014° in rad*/)); + } + if(q>5.93 && q<6.18){ + fh2_fib_ang_vs_x6->Fill(xposback,angMean); + fh2_fib_slope_vs_x6->Fill(xposback, -1000.*(atan(slope)/*+0.2443461014° in rad*/)); + } + } + + //if(1 != check_fit(slope,intercept,434.1, yy, 1138.5,RightHit)) return; + + + //if(r*r < 0.9) return; + + if(!IS_NAN(ypos30) /*&& q>5.93 && q<6.18*/){ + fh_f32X_vs_f30Y->Fill(xpos32,ypos30); + fh_f32Y_vs_f30Y->Fill(ypos32,ypos30); + fh_f32tot_vs_f30Y->Fill(tot32,ypos30); + fh_tofdq_vs_f30Y->Fill(q,ypos30); + } + + + if(!(q>qmin && q5.6 && q<6.5*/)/*r*r < 0.9*/) return; + fh_X_fib_vs_tofd->Fill(xmean, intercept+slope*1138.5); + fh_X_fib_vs_tofd_ang->Fill(xmean, xposback - (angX) * 602.95 /100. + (angX) * 1138.5 /100.); + + if((v2(1) > -0.2 && v2(1) < 0.2) || 1){ + + //double ang_check = (angX * 1.0238095 + 0.57857155 - angFibToFD)/angFibToFD; + double ang_check = angX * 1.0238095 + 0.57857155 - angFibToFD; + if( ( ((tot32 > 12. && totback > 12.) || 1) && !IS_NAN(slope) /*&& ang_check >= -1. && ang_check <= 1. */&& intercept > -4. && intercept <3.) || 1){ + fh2_ang_Fib_vs_ToFD->Fill(-angX,angFibToFD); + for (int i=0; i<100;i++){ + Double_t zrand = gRandom->Uniform(0., 1230.5/*cm*/);//Umlenkpunkt at 0. Fib32 at 434.1 cm, ToFD at 1138.5cm + //fh2_fibtracking_planeXZ->Fill(zrand, xpos32 - angX * 434.1 /1000. + angX * zrand /1000.); + //fh2_fibtracking_planeXZ->Fill(zrand, b + m * zrand); + fh2_fibtracking_planeXZ_nocuts->Fill(zrand, xposback - (angX) * 602.95 /100. + (angX) * zrand /100.); + fh2_fibtracking_planeXZ->Fill(zrand, intercept + slope * zrand ); + } + + } + + if( ( ((tot32 > 12. && totback > 12.) || 1) && intercept > -4. && intercept <3. + && !IS_NAN(slope) /*&& ang_check >= -1. && ang_check <= 1. */) || 1){ + fh2_FibToT_vs_QToFD->Fill(tot32,q); + + fh2_ToT_Fib32_vs_Fib30->Fill(tot32,tot30); + fh2_ToT_Fib32_vs_Fib31->Fill(tot32,tot31); + fh2_ToT_Fib32_vs_Fib33->Fill(tot32,tot33); + + + // cout<<"Q going into cut: "<Fill(v1(0),v1(1)); + fh_slope_vs_pos_rot->Fill(v2(0),v2(1)); + fh_pid->Fill(v1(0),q); + fh_pid_slope->Fill(v2(0),q); + fh_pid_3d->Fill(xposback,angMean,q); + //ntuple->Fill(xposback,angMean,q); + + + + } + + + + + } + } + } + + + + } + + } + +}//End FRS GATE + + //fNEvents += 1;//old +} + +void R3BFiberTrackingOnlineSpectra::FinishEvent() +{ + if (fFi30HitDataCA) + { + fFi30HitDataCA->Clear(); + } + + if (fFi31HitDataCA) + { + fFi31HitDataCA->Clear(); + } + + if (fFi32HitDataCA) + { + fFi32HitDataCA->Clear(); + } + + if (fFi33HitDataCA) + { + fFi33HitDataCA->Clear(); + } + + if(fHitFrs) + { + fHitFrs->Clear(); + } +} + +void R3BFiberTrackingOnlineSpectra::FinishTask() +{ + if (fFi32HitDataCA && (fFi33HitDataCA || fFi31HitDataCA)) + { + cTrackingXZ->Write(); + cAngVsX->Write(); + cslopeVsX->Write(); + cToFD_Z->Write(); + AngCorrFibToFD->Write(); + cAngRot->Write(); + cFibToT_Q->Write(); + cSlopeRot->Write(); + + fh2_fibtracking_planeXZ->Write(); + fh2_fib_ang_vs_x->Write(); + fh2_fib_ang_vs_x4->Write(); + fh2_fib_ang_vs_x5->Write(); + fh2_fib_ang_vs_x6->Write(); + fh2_fib_slope_vs_x->Write(); + fh_pid->Write(); + fh_pid_slope->Write(); + fh2_fib_slope_vs_x->Write(); + fh2_fib_slope_vs_x4->Write(); + fh2_fib_slope_vs_x5->Write(); + fh2_fib_slope_vs_x6->Write(); + fh2_FibToT_vs_QToFD->Write(); + fh_ang_vs_pos_rot->Write(); + fh_slope_vs_pos_rot->Write(); + fh2_ang_Fib_vs_ToFD->Write(); + fh_ToFD_Charge->Write(); + fh_pid_3d->Write(); + + } +} + +ClassImp(R3BFiberTrackingOnlineSpectra); diff --git a/analysis/online/R3BFiberTrackingOnlineSpectra.h b/analysis/online/R3BFiberTrackingOnlineSpectra.h new file mode 100644 index 000000000..3ebeb63fd --- /dev/null +++ b/analysis/online/R3BFiberTrackingOnlineSpectra.h @@ -0,0 +1,292 @@ +/****************************************************************************** + * Copyright (C) 2019 GSI Helmholtzzentrum für Schwerionenforschung GmbH * + * Copyright (C) 2019 Members of R3B Collaboration * + * * + * This software is distributed under the terms of the * + * GNU General Public Licence (GPL) version 3, * + * copied verbatim in the file "LICENSE". * + * * + * In applying this license GSI does not waive the privileges and immunities * + * granted to it by virtue of its status as an Intergovernmental Organization * + * or submit itself to any jurisdiction. * + ******************************************************************************/ + +// ------------------------------------------------------------ +// ----- R3BIncomingTrackingOnlineSpectra ----- +// ----- Created 25/02/22 by J.L. Rodriguez-Sanchez ----- +// ----- Fill tracking online histograms ----- +// ------------------------------------------------------------ + +#ifndef R3BFiberTrackingOnlineSpectra_H +#define R3BFiberTrackingOnlineSpectra_H 1 + +#include "FairTask.h" +#include "TCanvas.h" +#include "TH1.h" +#include "TH2F.h" +#include "TH3D.h" +#include "TNtuple.h" +#include "TMath.h" + +class TClonesArray; +class R3BEventHeader; +class R3BTGeoPar; +class R3BFiberMappingPar; + +class R3BFiberTrackingOnlineSpectra : public FairTask +{ + public: + /** + * Default constructor. + * Creates an instance of the task with default parameters. + */ + R3BFiberTrackingOnlineSpectra(); + + /** + * Standard constructor. + * Creates an instance of the task. + * @param name a name of the task. + * @param iVerbose a verbosity level. + */ + R3BFiberTrackingOnlineSpectra(const TString& name, Int_t iVerbose = 1); + + /** + * Destructor. + * Frees the memory used by the object. + */ + virtual ~R3BFiberTrackingOnlineSpectra(); + + /** + * Method for task initialization. + * This function is called by the framework before + * the event loop. + * @return Initialization status. kSUCCESS, kERROR or kFATAL. + */ + virtual InitStatus Init(); + + /** Virtual method ReInit **/ + virtual InitStatus ReInit(); + + virtual void SetParContainers(); + + /** + * Method for event loop implementation. + * Is called by the framework every time a new event is read. + * @param option an execution option. + */ + virtual void Exec(Option_t* option); + + /** + * A method for finish of processing of an event. + * Is called by the framework for each event after executing + * the tasks. + */ + virtual void FinishEvent(); + + /** + * Method for finish of the task execution. + * Is called by the framework after processing the event loop. + */ + virtual void FinishTask(); + + /** + * Methods to clean histograms. + */ + virtual void Reset_Histo(); + + /** + * Method to set up charge range in the histograms. + */ + inline void Set_Charge_range(Float_t minz, Float_t maxz) + { + fZ_min = minz; + fZ_max = maxz; + } + + void SetFRSGate(Bool_t FG) {fFRSGATE = FG;} + + void SetGateParams(Double_t aqmin, Double_t aqmax, Double_t qmin, Double_t qmax) + { + fAQmin = aqmin; + fAQmax = aqmax; + fQmin = qmin; + fQmax = qmax; + } + + void SetTPatMask(Int_t mask) {fTPatMask = mask;} + inline void SetTpat(Int_t tpat1, Int_t tpat2) + { + fTpat1 = tpat1; + fTpat2 = tpat2; + } + + + + int linreg(int n, const double x[], const double y[], double* m, double* b, double* r){ //Funktion for Linear regression + double sumx = 0.0; /* sum of x */ + double sumx2 = 0.0; /* sum of x**2 */ + double sumxy = 0.0; /* sum of x * y */ + double sumy = 0.0; /* sum of y */ + double sumy2 = 0.0; /* sum of y**2 */ + + for (int i=0;i 25.5) ) sensible_check++; + + if(hittype){ + if((m*fib2 + b < -a) || (m*fib2 + b > 50. + a) ) sensible_check++; + } + else{ + if((m*fib2 + b < -50+a) || (m*fib2 + b > a) ) sensible_check++; + } + + if((m*tofd + b < -60.) || (m*tofd + b > 60.) ) sensible_check++; + + if(sensible_check==0) sensible = 1; + + //if(sensible_check > 0){ + // if(hittype) std::cout<<"check "<*/ + TClonesArray* fFi32HitDataCA; /**< Array with Fiber Hit-input data. >*/ + TClonesArray* fFi33HitDataCA; /**< Array with Fiber Hit-input data. >*/ + TClonesArray* fHitTofdItems; + TClonesArray* fHitFrs; + TClonesArray* fCalLos; + + // Parameters + R3BTGeoPar* fMw0GeoPar; + R3BTGeoPar* fTargetGeoPar; + R3BTGeoPar* fMw1GeoPar; + + // check for trigger should be done globablly (somewhere else) + + Int_t fTpat1; + Int_t fTpat2; + Int_t fTrigger; + Int_t fTPatMask; + R3BEventHeader* header; /**< Event header. */ + Int_t fNEvents; /**< Event counter. */ + Float_t fPosTarget; + Float_t fWidthTarget; + Float_t fDist_acelerator_glad; + Float_t fZ_max, fZ_min; + Bool_t fFRSGATE; + Double_t fAQmin, fAQmax, fQmin, fQmax; + + // Canvas + TCanvas* cTrackingXZ; + TCanvas* cAngVsX; + TCanvas* cslopeVsX; + TCanvas* cToFD_Z; + TCanvas* AngCorrFibToFD; + TCanvas* cAngRot; + TCanvas* cFibToT_Q; + TCanvas* cSlopeRot; + TCanvas* cFibToTCorr; + TCanvas* cFibToFDnocuts; + TCanvas* cPID_3D; + TCanvas* cFibInfo; + TCanvas* cFib30Rel; + TCanvas* cXcorr; + TCanvas* cDeltaQ; + TCanvas* cMult; + + // Histograms for Hit data + + TH2F* fh_f32X_vs_f30Y; + TH2F* fh_f32Y_vs_f30Y; + TH2F* fh_f32tot_vs_f30Y; + TH2F* fh_tofdq_vs_f30Y; + TH2F* fh_ToT_Fib; + TH2F* fh2_FibToT_vs_QToFD; + TH2F* fh_ang_vs_pos_rot; + TH2F* fh_slope_vs_pos_rot; + TH2F* fh_pid; + TH2F* fh_pid_slope; + TH2F* fh2_fibtracking_planeXZ; + TH2F* fh2_fibtracking_planeXZ_nocuts; + TH2F* fh2_ang_Fib_vs_ToFD; + TH2F* fh2_ang_Fib_vs_ToFD_nocuts; + TH2F* fh2_fib_slope_vs_x; + TH2F* fh2_fib_slope_vs_x_nocuts; + TH2F* fh2_fib_slope_vs_x4; + TH2F* fh2_fib_slope_vs_x5; + TH2F* fh2_fib_slope_vs_x6; + TH2F* fh2_fib_ang_vs_x; + TH2F* fh2_fib_ang_vs_x_nocuts; + TH2F* fh2_fib_ang_vs_x4; + TH2F* fh2_fib_ang_vs_x5; + TH2F* fh2_fib_ang_vs_x6; + TH1F* fh_ToFD_Charge; + TH1F* fh_Rsqr; + TH2F* fh2_ToT_Fib32_vs_Fib31; + TH2F* fh2_ToT_Fib32_vs_Fib33; + TH2F* fh2_ToT_Fib32_vs_Fib30; + TH2F* fh_X_fib_vs_tofd; + TH2F* fh_X_fib_vs_tofd_no_cut; + TH2F* fh_X_fib_vs_tofd_ang; + TH2F* fh_X_fib_vs_tofd_ang_no_cut; + TH2F* fh_pos_TPat; + TH2F* fh_x_vs_tof; + TH2F* fh_x_vs_q; + TH2F* fh_tofd_charge_planes[6]; + TH2F* fh_ToF_vs_charge; + TH1D* fh_mult_coinc; + TH1D* fh_mult_minus_los; + TH1D* fh_mult_fib30; + TH1D* fh_mult_fib31; + TH1D* fh_mult_fib32; + TH1D* fh_mult_fib33; + TH1D* fh_mult_fibAll; + TH3* fh_pid_3d; + TNtuple* ntuple; + + public: + ClassDef(R3BFiberTrackingOnlineSpectra, 1) +}; + +#endif /* R3BFiberTrackingOnlineSpectra_H */ diff --git a/fiber/mapmt/R3BFiberMAPMTCal2Hit.cxx b/fiber/mapmt/R3BFiberMAPMTCal2Hit.cxx index abb03070d..606221683 100644 --- a/fiber/mapmt/R3BFiberMAPMTCal2Hit.cxx +++ b/fiber/mapmt/R3BFiberMAPMTCal2Hit.cxx @@ -1,6 +1,6 @@ /****************************************************************************** * Copyright (C) 2019 GSI Helmholtzzentrum für Schwerionenforschung GmbH * - * Copyright (C) 2019-2023 Members of R3B Collaboration * + * Copyright (C) 2019 Members of R3B Collaboration * * * * This software is distributed under the terms of the * * GNU General Public Licence (GPL) version 3, * @@ -73,7 +73,7 @@ R3BFiberMAPMTCal2Hit::R3BFiberMAPMTCal2Hit(const char* a_name, , fnEvents(0) , ftofmin(-1000) , ftofmax(1000) - , fWrite(false) + , fWrite(true) , fGate_ns(100.) , fOnline(kFALSE) , fOrientation(STANDARD) @@ -390,12 +390,29 @@ void R3BFiberMAPMTCal2Hit::Exec(Option_t* option) Double_t t_down = down_tot.lead_ns; Double_t t_up = up_tot.lead_ns; Double_t dtime = fTimeStitch->GetTime(t_up - t_down, "clocktdc", "clocktdc"); - Double_t tof = - fHeader ? fTimeStitch->GetTime((t_up + t_down) / 2. - fHeader->GetTStart(), "vftx", "clocktdc") - : (t_up + t_down) / 2.; - + + Double_t tof0 = fHeader ? fTimeStitch->GetTime(t_up - fHeader->GetTStart(), "vftx", "clocktdc") + : t_up; + Double_t tof1 = fHeader ? fTimeStitch->GetTime(t_down - fHeader->GetTStart(), "vftx", "clocktdc") + : t_down; + + Double_t tof00 = (tof0 + tof1) / 2.; + + Double_t tof01 = fHeader ? fTimeStitch->GetTime((t_up+t_down)/2., "vftx", "clocktdc") + : 0./0.; + Double_t tof11 = fHeader ? fTimeStitch->GetTime(fHeader->GetTStart(), "vftx", "clocktdc") + : 0./0.; + + //Double_t tof = tof01 - tof11; + Double_t tof = fHeader ? fTimeStitch->GetTime((t_up+t_down)/2. - fHeader->GetTStart(), "vftx", "clocktdc") + : 0./0.; + + // if(!TMath::IsNaN(tof)) cout<<"tof "<GetTStart()<<" x "<<(t_up+t_down)/2. - fHeader->GetTStart()<GetTStart(); + //if(tof != tof00) return; // Fill histograms for gain match, offset and sync. - if (fWrite) + + if (fWrite /*&& cal_num == 1*/) { fh_ToT_bottom_Fib_raw->Fill(fiber_id, tot_down); fh_ToT_top_Fib_raw->Fill(fiber_id, tot_up); @@ -431,8 +448,9 @@ void R3BFiberMAPMTCal2Hit::Exec(Option_t* option) // t_up -= offsetUp; // t_up -= tsync; dtime -= offsetDT; + // cout<<"tofb "<= ftofmin && tof <= ftofmax) + if ((tof >= ftofmin && tof <= ftofmax ) || 1) { new ((*fHitItems)[fHitItems->GetEntriesFast()]) R3BFiberMAPMTHitData( fDetId, x, y, eloss, tof, fiber_id, t_down, t_up, tot_down, tot_up); @@ -561,8 +579,7 @@ void R3BFiberMAPMTCal2Hit::FinishTask() if (fIsCalibrator) { Bool_t Redo = false; // Add a redo flag, to redo a calibration, if Maximum falls on noise. - Double_t PercentOfMax = 0.5; - + Double_t PercentOfMax = 0.9; R3BFiberMAPMTHitModulePar* mpar; for (UInt_t i = 1; i <= fNumFibers; i++) @@ -603,7 +620,7 @@ void R3BFiberMAPMTCal2Hit::FinishTask() R3BFiberMAPMTHitModulePar* par = fCalPar->GetModuleParAt(i); par->SetGainDown(proj->GetBinCenter(j)); - R3BLOG(info, fName << " fiberId: " << i << ",gainDown: " << proj->GetBinCenter(j)); + if(i<10){ R3BLOG(info, fName << " fiberId: " << i << ",gainDown: " << proj->GetBinCenter(j));} if (Redo == true) { PartMax = proj->GetMaximum() * PercentOfMax; @@ -640,7 +657,7 @@ void R3BFiberMAPMTCal2Hit::FinishTask() { R3BFiberMAPMTHitModulePar* par1 = fCalPar->GetModuleParAt(i); par1->SetGainUp(proj1->GetBinCenter(j)); - R3BLOG(info, fName << " fiberId: " << i << ",gainUp: " << proj1->GetBinCenter(j)); + if(i<10) {R3BLOG(info, fName << " fiberId: " << i << ",gainUp: " << proj1->GetBinCenter(j));} break; } } @@ -652,19 +669,19 @@ void R3BFiberMAPMTCal2Hit::FinishTask() // par2->SetOffsetUp(0.5 * proj2->GetBinCenter(proj2->GetMaximumBin())); auto Max = proj2->GetXaxis()->GetBinCenter(proj2->GetMaximumBin()); TF1* fGauss = new TF1("fgaus", "gaus", Max - 8., Max + 8.); - proj2->Fit("fgaus", "QR0"); + //proj2->Fit("fgaus", "QR0"); auto offsetdt = fGauss->GetParameter(1); par2->SetOffsetDown(offsetdt); // fGauss->Draw("LSAME"); // proj2->Write(); // tsync - // R3BFiberMAPMTHitModulePar* par3 = fCalPar->GetModuleParAt(i); + R3BFiberMAPMTHitModulePar* par3 = fCalPar->GetModuleParAt(i); auto proj3 = (TH1F*)fh_Fib_ToF_raw->ProjectionY("", i, i, 0); auto sync = proj3->GetXaxis()->GetBinCenter(proj3->GetMaximumBin()); - par2->SetSync(sync); + par3->SetSync(sync); - R3BLOG(info, fName << " fiberId: " << i << ", offset_DT: " << offsetdt << ", sync: " << sync); + if(i<10) {R3BLOG(info, fName << " fiberId: " << i << ", offset_DT: " << offsetdt << ", sync: " << sync);} } fCalPar->setChanged(); diff --git a/fiber/mapmt/R3BFiberMAPMTCal2Hit.h b/fiber/mapmt/R3BFiberMAPMTCal2Hit.h index 16bb714ae..469e10a39 100644 --- a/fiber/mapmt/R3BFiberMAPMTCal2Hit.h +++ b/fiber/mapmt/R3BFiberMAPMTCal2Hit.h @@ -1,6 +1,6 @@ /****************************************************************************** * Copyright (C) 2019 GSI Helmholtzzentrum für Schwerionenforschung GmbH * - * Copyright (C) 2019-2023 Members of R3B Collaboration * + * Copyright (C) 2019 Members of R3B Collaboration * * * * This software is distributed under the terms of the * * GNU General Public Licence (GPL) version 3, * @@ -14,7 +14,7 @@ #ifndef R3BFIBERMAPMTCAL2HIT_H #define R3BFIBERMAPMTCAL2HIT_H 1 -#include +#include "TClonesArray.h" #include "FairTask.h" #include diff --git a/fiber/mapmt/R3BFiberMAPMTMapped2Cal.cxx b/fiber/mapmt/R3BFiberMAPMTMapped2Cal.cxx index 3f0bf1be4..977b86ee5 100644 --- a/fiber/mapmt/R3BFiberMAPMTMapped2Cal.cxx +++ b/fiber/mapmt/R3BFiberMAPMTMapped2Cal.cxx @@ -143,9 +143,9 @@ void R3BFiberMAPMTMapped2Cal::Exec(Option_t* option) Double_t time_ns = -1; if (fine_ns < 0. || fine_ns > fClockFreq) { - R3BLOG(error, - "(" << fName << "): Channel=" << channel << ": Bad CTDC fine time (raw=" << fine_raw - << ",ns=" << fine_ns << ")."); + //R3BLOG(error, + // "(" << fName << "): Channel=" << channel << ": Bad CTDC fine time (raw=" << fine_raw + // << ",ns=" << fine_ns << ")."); continue; } diff --git a/tofd/calibration/R3BTofDCal2Hit.cxx b/tofd/calibration/R3BTofDCal2Hit.cxx index 746b30c21..f30ad7362 100644 --- a/tofd/calibration/R3BTofDCal2Hit.cxx +++ b/tofd/calibration/R3BTofDCal2Hit.cxx @@ -1,6 +1,6 @@ /****************************************************************************** * Copyright (C) 2019 GSI Helmholtzzentrum für Schwerionenforschung GmbH * - * Copyright (C) 2019-2023 Members of R3B Collaboration * + * Copyright (C) 2019 Members of R3B Collaboration * * * * This software is distributed under the terms of the * * GNU General Public Licence (GPL) version 3, * @@ -26,13 +26,14 @@ #include "R3BTofDMappingPar.h" #include "R3BTofdCalData.h" #include "R3BTofdHitData.h" -#include +#include #include "FairLogger.h" #include "FairRuntimeDb.h" #include "TH1F.h" #include "TH2F.h" #include "THnSparse.h" +#include "TFile.h" #include "TClonesArray.h" #include "TMath.h" @@ -155,10 +156,10 @@ R3BTofDCal2Hit::~R3BTofDCal2Hit() void R3BTofDCal2Hit::SetParContainers() { - fMapPar = dynamic_cast(FairRuntimeDb::instance()->getContainer("tofdMappingPar")); - R3BLOG_IF(warn, !fMapPar, "Could not get access to tofdMappingPar container"); + fMapPar = (R3BTofDMappingPar*)FairRuntimeDb::instance()->getContainer("tofdMappingPar"); + R3BLOG_IF(warning, !fMapPar, "Could not get access to tofdMappingPar container"); - fHitPar = dynamic_cast(FairRuntimeDb::instance()->getContainer("tofdHitPar")); + fHitPar = (R3BTofDHitPar*)FairRuntimeDb::instance()->getContainer("tofdHitPar"); if (!fHitPar) { R3BLOG(error, "Could not get access to tofdHitPar container"); @@ -191,13 +192,13 @@ InitStatus R3BTofDCal2Hit::Init() return kFATAL; } - header = dynamic_cast(mgr->GetObject("EventHeader.")); + header = (R3BEventHeader*)mgr->GetObject("EventHeader."); R3BLOG_IF(fatal, NULL == header, "EventHeader. not found"); - fCalItems = dynamic_cast(mgr->GetObject("TofdCal")); + fCalItems = (TClonesArray*)mgr->GetObject("TofdCal"); R3BLOG_IF(fatal, NULL == fCalItems, "TofdCal not found"); - fCalTriggerItems = dynamic_cast(mgr->GetObject("TofdTriggerCal")); + fCalTriggerItems = (TClonesArray*)mgr->GetObject("TofdTriggerCal"); R3BLOG_IF(fatal, NULL == fCalTriggerItems, "TofdTriggerCal not found"); maxevent = mgr->CheckMaxEventNo(); @@ -317,7 +318,7 @@ void R3BTofDCal2Hit::Exec(Option_t* option) // puts("Event"); for (Int_t ihit = 0; ihit < nHits; ihit++) { - auto* hit = dynamic_cast(fCalItems->At(ihit)); + auto* hit = (R3BTofdCalData*)fCalItems->At(ihit); size_t idx = (hit->GetDetectorId() - 1) * fPaddlesPerPlane + (hit->GetBarId() - 1); // std::cout << "Hits: " << hit->GetDetectorId() << ' ' << hit->GetBarId() << ' ' << hit->GetSideId() << ' ' @@ -334,7 +335,7 @@ void R3BTofDCal2Hit::Exec(Option_t* option) std::vector trig_map; for (int i = 0; i < fCalTriggerItems->GetEntriesFast(); ++i) { - auto trig = dynamic_cast(fCalTriggerItems->At(i)); + auto trig = (R3BTofdCalData const*)fCalTriggerItems->At(i); if (trig_map.size() < trig->GetBarId()) { trig_map.resize(trig->GetBarId()); @@ -388,9 +389,9 @@ void R3BTofDCal2Hit::Exec(Option_t* option) { if (!s_was_trig_missing) { - R3BLOG(error, "Missing trigger information!"); - R3BLOG(error, "Top: " << top->GetDetectorId() << ' ' << top->GetSideId() << ' ' << top->GetBarId()); - R3BLOG(error, "Bot: " << bot->GetDetectorId() << ' ' << bot->GetSideId() << ' ' << bot->GetBarId()); + //R3BLOG(error, "Missing trigger information!"); + //R3BLOG(error, "Top: " << top->GetDetectorId() << ' ' << top->GetSideId() << ' ' << top->GetBarId()); + //R3BLOG(error, "Bot: " << bot->GetDetectorId() << ' ' << bot->GetSideId() << ' ' << bot->GetBarId()); s_was_trig_missing = true; } ++n2; @@ -434,10 +435,24 @@ void R3BTofDCal2Hit::Exec(Option_t* option) auto top_tot = fmod(top->GetTimeTrailing_ns() - top->GetTimeLeading_ns() + c_range_ns, c_range_ns); auto bot_tot = fmod(bot->GetTimeTrailing_ns() - bot->GetTimeLeading_ns() + c_range_ns, c_range_ns); - auto THit_raw = (bot->GetTimeLeading_ns() + top->GetTimeLeading_ns()) / 2.; // needed for TOF for ROLUs - - // std::cout<<"ToT: "<GetTimeLeading_ns() + top->GetTimeLeading_ns()) / 2.; // needed for TOF for ROLUs + auto THit_raw = fTimeStitch->GetTime((bot_trig_ns + top_trig_ns) / 2. - header->GetTStart(),"tamex","vftx"); // needed for TOF for ROLUs + + //std::cout<<"ToT: "<GetBarId()<<" BarB: "<GetBarId()<GetBarId() == bot->GetBarId()){ + bool both = false; + nev++; + if(top_tot > c_range_ns/10 && bot_tot > c_range_ns/10){ nev_lmboth++; both = true;} + if(top_tot > c_range_ns/10){ + top_tot = bot_tot; + if(!both) nev_lmt++; + } + if(bot_tot > c_range_ns/10){ + bot_tot = top_tot; + if(!both) nev_lmb++; + } + } // register multi hits vmultihits[iPlane][iBar] += 1; @@ -448,6 +463,8 @@ void R3BTofDCal2Hit::Exec(Option_t* option) continue; } + Double_t Ctop_tot = top_tot * par->GetToTOffset2(); + Double_t Cbot_tot = bot_tot * par->GetToTOffset1(); // walk corrections if (par->GetPar1Walk() == 0. || par->GetPar2Walk() == 0. || par->GetPar3Walk() == 0. || par->GetPar4Walk() == 0. || par->GetPar5Walk() == 0.) @@ -456,13 +473,13 @@ void R3BTofDCal2Hit::Exec(Option_t* option) } else { - auto bot_ns_walk = bot_ns - walk(bot_tot, + auto bot_ns_walk = bot_ns - walk(Cbot_tot, par->GetPar1Walk(), par->GetPar2Walk(), par->GetPar3Walk(), par->GetPar4Walk(), par->GetPar5Walk()); - auto top_ns_walk = top_ns - walk(top_tot, + auto top_ns_walk = top_ns - walk(Ctop_tot, par->GetPar1Walk(), par->GetPar2Walk(), par->GetPar3Walk(), @@ -475,6 +492,7 @@ void R3BTofDCal2Hit::Exec(Option_t* option) // calculate time of hit Double_t THit = (bot_ns + top_ns) / 2. - par->GetSync(); + /// cout<<"GetSynch "<GetSync()<<" tb "<GetSync()<GetLambda() * log((top_tot * par->GetToTOffset2()) / (bot_tot * par->GetToTOffset1())); + par->GetLambda() * log((Ctop_tot) / (Cbot_tot)); + // par->GetLambda() * log((top_tot) / (bot_tot)); + /// cout<<"Lambda "<GetLambda()<<" GetVeff "<GetVeff()< 0) + if (fTofdQ > 0 /*&& 1*/) { if (fTofdTotPos) { @@ -525,8 +545,9 @@ void R3BTofDCal2Hit::Exec(Option_t* option) para[1] = par->GetPolb(); para[2] = par->GetPolc(); para[3] = par->GetPold(); - qb = TMath::Sqrt(top_tot * bot_tot) / - (para[0] + para[1] * pos + para[2] * pow(pos, 2) + para[3] * pow(pos, 3)); + // qb = TMath::Sqrt(Ctop_tot * Cbot_tot) / + qb = TMath::Sqrt(Ctop_tot * Cbot_tot) * + (para[0] + para[1] * pos + para[2] * pow(pos, 2) + para[3] * pow(pos, 3)/*1.*/); qb = qb * fTofdQ; } else @@ -536,13 +557,13 @@ void R3BTofDCal2Hit::Exec(Option_t* option) para[1] = par->GetPar1b(); para[2] = par->GetPar1c(); para[3] = par->GetPar1d(); - auto q1 = bot_tot / + auto q1 = Cbot_tot / (para[0] * (exp(-para[1] * (pos + 100.)) + exp(-para[2] * (pos + 100.))) + para[3]); para[0] = par->GetPar2a(); para[1] = par->GetPar2b(); para[2] = par->GetPar2c(); para[3] = par->GetPar2d(); - auto q2 = top_tot / + auto q2 = Ctop_tot / (para[0] * (exp(-para[1] * (pos + 100.)) + exp(-para[2] * (pos + 100.))) + para[3]); q1 = q1 * fTofdQ; q2 = q2 * fTofdQ; @@ -551,9 +572,12 @@ void R3BTofDCal2Hit::Exec(Option_t* option) } else { - qb = TMath::Sqrt(top_tot * bot_tot); + qb = TMath::Sqrt(Ctop_tot * Cbot_tot); + //qb = fTofdQ*(Ctop_tot + Cbot_tot) / 2.; } - + qb = qb / 200.; + if(qb > 70.) qb = qb / 200.; + // std::cout<<"qb "<GetPar1za(); parz[1] = par->GetPar1zb(); @@ -578,16 +602,25 @@ void R3BTofDCal2Hit::Exec(Option_t* option) LOG(debug) << "y in this event " << pos << " plane " << iPlane << " ibar " << iBar << "\n"; // Tof with respect LOS detector - auto tof = fTimeStitch->GetTime((bot_ns + top_ns) / 2. - header->GetTStart(), "tamex", "vftx"); - // auto tof_corr = par->GetTofSyncOffset() + par->GetTofSyncSlope() * tof; - auto tof_corr = tof - par->GetTofSyncOffset(); + //auto tof = fTimeStitch->GetTime((bot_ns + top_ns) / 2. , "tamex", "vftx") - fTimeStitch->GetTime(header->GetTStart(), "tamex", "vftx");(bot_ns + par->GetOffset1()) - (top_ns + par->GetOffset2()) + auto tof = fTimeStitch->GetTime((bot_ns + par->GetOffset1() + top_ns+par->GetOffset2()) / 2.,"tamex","tamex") - fTimeStitch->GetTime(header->GetTStart(), "tamex", "vftx"); + auto tof_corr = -par->GetTofSyncOffset() + /*par->GetTofSyncSlope() **/ tof; + // cout<<"tof_raw "<GetTofSyncOffset()<<" slope "<GetTofSyncSlope()<<" tof_corr "<GetTofSyncOffset() + par->GetTofSyncSlope() * tof< 0) // { event.push_back( - { parz[0] + parz[1] * qb + parz[2] * qb * qb, THit, xp, pos, iPlane, iBar, THit_raw, tof_corr }); + //{ parz[0] + parz[1] * qb + parz[2] * qb * qb, THit, xp, pos, iPlane, iBar, THit_raw, tof_corr }); + { (parz[0] + (parz[1] * qb)/* - 0.153*/ + parz[2] * qb * qb), THit, xp, pos, iPlane, iBar, THit_raw, tof_corr }); + // { qb, THit, xp, pos, iPlane, iBar, THit_raw, tof_corr }); // } - + //if((parz[0] + (parz[1] * qb) - 0.153 + parz[2] * qb * qb) > 0.){ + // cout<<"p1 "< 0 && parz[2] > 0) { event.push_back( @@ -1063,6 +1096,9 @@ void R3BTofDCal2Hit::FinishTask() sprint << "n1=" << n1 << " n2=" << n2; R3BLOG(info, sprint.str()); + + cout<<"missing leading edges:\n"; + cout<<"nev: "<GetXaxis()->SetTitle("Bar"); fh_tofd_time_los_h2[i]->GetYaxis()->SetTitle("ToF [ns]"); fh_tofd_time_los_h2[i]->GetXaxis()->CenterTitle(true); diff --git a/tofd/pars/R3BTofDCal2HitPar.cxx b/tofd/pars/R3BTofDCal2HitPar.cxx index 1f3ec401a..9d742c1fc 100644 --- a/tofd/pars/R3BTofDCal2HitPar.cxx +++ b/tofd/pars/R3BTofDCal2HitPar.cxx @@ -1,6 +1,6 @@ /****************************************************************************** * Copyright (C) 2019 GSI Helmholtzzentrum für Schwerionenforschung GmbH * - * Copyright (C) 2019-2023 Members of R3B Collaboration * + * Copyright (C) 2019 Members of R3B Collaboration * * * * This software is distributed under the terms of the * * GNU General Public Licence (GPL) version 3, * @@ -21,7 +21,7 @@ #include "R3BTofdCalData.h" #include "FairLogger.h" -#include "FairRootManager.h" +#include #include "FairRtdbRun.h" #include "FairRunIdGenerator.h" #include "FairRuntimeDb.h" @@ -36,7 +36,9 @@ #include "TLine.h" #include "TMath.h" #include "TProfile.h" +#include "TFile.h" #include "TSpectrum.h" +#include "TSystem.h" #include "TVirtualFitter.h" #include @@ -70,8 +72,11 @@ R3BTofDCal2HitPar::R3BTofDCal2HitPar(const char* name, Int_t iVerbose) , fHitPar(NULL) , fTofdY(0.) , fTofdQ(0.) - , fMaxQ(1500) - , fNbZPeaks(1) + , fMinQ(180) + , fMaxQ(300) + , fToTMin(190.) + , fToTMax(210.) + , fNbZPeaks(2) , fZfitType("pol1") , fTofdTotLow(0.) , fTofdTotHigh(200.) @@ -81,123 +86,169 @@ R3BTofDCal2HitPar::R3BTofDCal2HitPar(const char* name, Int_t iVerbose) , fMeanTof(20.) , fHeader(nullptr) , maxevent(0) + , fHistoFile("") { + if(fParameter != 5 && fParameter !=-3 && fParameter!=7){ for (Int_t i = 0; i < fNofPlanes; i++) { - fh_tofd_TotPm[i] = NULL; + fh_tofd_TotPm[i][0] = NULL; + fh_tofd_TotPm[i][1]=NULL; + fh_tofd_TotPm[i][2]=NULL; + fh_tofd_TotPm1[i] = NULL; + fh_tofd_TotPm2[i] = NULL; fhTdiff[i] = NULL; fhTsync[i] = NULL; + fhTofsync[i] = NULL; + //fhlos[i] = NULL; for (Int_t j = 0; j < fPaddlesPerPlane; j++) { + //fhlost[i][j] = NULL; + fhtott[i][j] = NULL; fhLogTot1vsLogTot2[i][j] = NULL; fhSqrtQvsPosToT[i][j] = NULL; fhQvsPos[i][j] = NULL; fhTot1vsPos[i][j] = NULL; fhTot2vsPos[i][j] = NULL; fh1_tofsync[i][j] = NULL; - // fhTot1vsTot2[i][j] = NULL; - } + //fhTot1vsTot2[i][j] = NULL; + fh1_walk_b[i][j] = NULL; + fh1_walk_t[i][j] = NULL; + //fhZoffset[i][j]=NULL; + } } + } } R3BTofDCal2HitPar::~R3BTofDCal2HitPar() { - for (Int_t i = 0; i < fNofPlanes; i++) - { - if (fh_tofd_TotPm[i]) - delete fh_tofd_TotPm[i]; - if (fhTdiff[i]) - delete fhTdiff[i]; - if (fhTsync[i]) - delete fhTsync[i]; - for (Int_t j = 0; j < fPaddlesPerPlane; j++) - { - if (fhLogTot1vsLogTot2[i][j]) - delete fhLogTot1vsLogTot2[i][j]; - if (fhSqrtQvsPosToT[i][j]) - delete fhSqrtQvsPosToT[i][j]; - if (fhQvsPos[i][j]) - delete fhQvsPos[i][j]; - if (fhTot1vsPos[i][j]) - delete fhTot1vsPos[i][j]; - if (fhTot2vsPos[i][j]) - delete fhTot2vsPos[i][j]; - if (fh1_tofsync[i][j]) - delete fh1_tofsync[i][j]; - /* - if (fhTot1vsTot2[i][j]) - delete fhTot1vsTot2[i][j]; - */ - } - } - if (fHitPar) - { - delete fHitPar; - } + if(fParameter != 5 && fParameter != -3 && fParameter!=7){ + for (Int_t i = 0; i < fNofPlanes; i++) + { + //if(fhlos[i]) + //delete fhlos[i]; + if(fhTofsync[i]) + delete fhTofsync[i]; + if (fh_tofd_TotPm[i][0]) + delete fh_tofd_TotPm[i][0]; + if (fh_tofd_TotPm[i][1]) + delete fh_tofd_TotPm[i][1]; + if (fh_tofd_TotPm[i][2]) + delete fh_tofd_TotPm[i][2]; + if (fhTdiff[i]) + delete fhTdiff[i]; + if (fhTsync[i]) + delete fhTsync[i]; + for (Int_t j = 0; j < fPaddlesPerPlane; j++) + { + //if(fhlost[i][j]) + //delete fhlost[i][j]; + if(fhtott[i][j]) + delete fhtott[i][j]; + if (fhLogTot1vsLogTot2[i][j]) + delete fhLogTot1vsLogTot2[i][j]; + if (fhSqrtQvsPosToT[i][j]) + delete fhSqrtQvsPosToT[i][j]; + if (fhQvsPos[i][j]) + delete fhQvsPos[i][j]; + if (fhTot1vsPos[i][j]) + delete fhTot1vsPos[i][j]; + if (fhTot2vsPos[i][j]) + delete fhTot2vsPos[i][j]; + if (fh1_tofsync[i][j]) + delete fh1_tofsync[i][j]; + if( fh1_walk_t[i][j]) + delete fh1_walk_t[i][j]; + if( fh1_walk_b[i][j]) + delete fh1_walk_b[i][j]; + + // if (fhZoffset[i][j]) + // delete fhZoffset[i][j]; + + } + } + if (fHitPar) + { + delete fHitPar; + } + } } InitStatus R3BTofDCal2HitPar::Init() { - R3BLOG(info, ""); - FairRootManager* rm = FairRootManager::Instance(); - if (!rm) - { - R3BLOG(fatal, "FairRootManager not found"); - return kFATAL; - } - - fHeader = dynamic_cast(rm->GetObject("EventHeader.")); - R3BLOG_IF(fatal, NULL == fHeader, "EventHeader. not found"); - - fCalData = dynamic_cast(rm->GetObject("TofdCal")); - R3BLOG_IF(fatal, NULL == fCalData, "TofdCal not found"); - - fCalTriggerItems = dynamic_cast(rm->GetObject("TofdTriggerCal")); - R3BLOG_IF(fatal, NULL == fCalTriggerItems, "TofdTriggerCal not found"); - - fHitPar = dynamic_cast(FairRuntimeDb::instance()->getContainer("tofdHitPar")); - if (!fHitPar) - { - R3BLOG(error, "Could not get access to tofdHitPar container"); - return kFATAL; - } - - for (Int_t i = 1; i <= fNofPlanes; i++) - for (Int_t j = 1; j <= fPaddlesPerPlane; j++) - CreateHistograms(i, j); - - // Definition of a time stich object to correlate times coming from different systems - fTimeStitch = new R3BCoarseTimeStitch(); - - return kSUCCESS; + R3BLOG(info, ""); + FairRootManager* rm = FairRootManager::Instance(); + if (!rm) + { + R3BLOG(fatal, "FairRootManager not found"); + return kFATAL; + } + + fHeader = (R3BEventHeader*)rm->GetObject("EventHeader."); + R3BLOG_IF(fatal, NULL == fHeader, "EventHeader. not found"); + + fCalData = (TClonesArray*)rm->GetObject("TofdCal"); + R3BLOG_IF(fatal, NULL == fCalData, "TofdCal not found"); + + fCalTriggerItems = (TClonesArray*)rm->GetObject("TofdTriggerCal"); + R3BLOG_IF(fatal, NULL == fCalTriggerItems, "TofdTriggerCal not found"); + + fHitPar = (R3BTofDHitPar*)FairRuntimeDb::instance()->getContainer("tofdHitPar"); + if (!fHitPar) + { + R3BLOG(error, "Could not get access to tofdHitPar container"); + return kFATAL; + } + if(fParameter != 5 && fParameter != -3){ + for (Int_t i = 1; i <= fNofPlanes; i++){ + for (Int_t j = 1; j <= fPaddlesPerPlane; j++){ + SetHistogramsNULL(i, j); + } + } + + for (Int_t i = 1; i <= fNofPlanes; i++){ + for (Int_t j = 1; j <= fPaddlesPerPlane; j++){ + if(j==1) SetHistogramsNULL(i, j); + CreateHistograms(i, j); + } + } + } + // Definition of a time stich object to correlate times coming from different systems + fTimeStitch = new R3BCoarseTimeStitch(); + + + + return kSUCCESS; } void R3BTofDCal2HitPar::SetParContainers() { - fMapPar = dynamic_cast(FairRuntimeDb::instance()->getContainer("tofdMappingPar")); - R3BLOG_IF(warn, !fMapPar, "Could not get access to tofdMappingPar container"); -} + fMapPar = (R3BTofDMappingPar*)FairRuntimeDb::instance()->getContainer("tofdMappingPar"); + R3BLOG_IF(warning, !fMapPar, "Could not get access to tofdMappingPar container"); + } namespace { - uint64_t n1, n2; + uint64_t n1, n2; }; void R3BTofDCal2HitPar::Exec(Option_t* option) { - // test for requested trigger (if possible) - if ((fTrigger >= 0) && (fHeader) && (fHeader->GetTrigger() != fTrigger)) - return; + //auto parX = fHitPar->GetModuleParAt(1,1); + if(fParameter == 5 || fParameter == -3 || fParameter==7) return; + // test for requested trigger (if possible) + if ((fTrigger >= 0) && (fHeader) && (fHeader->GetTrigger() != fTrigger)) + return; + + // fTpat = 1-16; fTpat_bit = 0-15 + Int_t fTpat_bit = fTpat - 1; + if (fTpat_bit >= 0) + { + Int_t itpat = fHeader->GetTpat(); + Int_t tpatvalue = (itpat & (1 << fTpat_bit)) >> fTpat_bit; + if ((fHeader) && (tpatvalue == 0)) + return; + } - // fTpat = 1-16; fTpat_bit = 0-15 - Int_t fTpat_bit = fTpat - 1; - if (fTpat_bit >= 0) - { - Int_t itpat = fHeader->GetTpat(); - Int_t tpatvalue = (itpat & (1 << fTpat_bit)) >> fTpat_bit; - if ((fHeader) && (tpatvalue == 0)) - return; - } // ToFD detector std::vector> multihits(fNofPlanes, std::vector(fPaddlesPerPlane)); @@ -206,6 +257,7 @@ void R3BTofDCal2HitPar::Exec(Option_t* option) multihits[i][j] = 0; UInt_t nHits = fCalData->GetEntriesFast(); + // Organize cals into bars. struct Entry { @@ -215,7 +267,7 @@ void R3BTofDCal2HitPar::Exec(Option_t* option) std::map bar_map; for (Int_t ihit = 0; ihit < nHits; ihit++) { - auto* hit = dynamic_cast(fCalData->At(ihit)); + auto* hit = (R3BTofdCalData*)fCalData->At(ihit); size_t idx = (hit->GetDetectorId() - 1) * fPaddlesPerPlane + (hit->GetBarId() - 1); auto ret = bar_map.insert(std::pair(idx, Entry())); auto& vec = 1 == hit->GetSideId() ? ret.first->second.top : ret.first->second.bot; @@ -247,8 +299,8 @@ void R3BTofDCal2HitPar::Exec(Option_t* option) Double_t top_trig_ns = 0, bot_trig_ns = 0; if (top_trig_i < trig_num && bot_trig_i < trig_num) { - auto top_trig = dynamic_cast(fCalTriggerItems->At(top_trig_i)); - auto bot_trig = dynamic_cast(fCalTriggerItems->At(bot_trig_i)); + auto top_trig = (R3BTofdCalData const*)fCalTriggerItems->At(top_trig_i); + auto bot_trig = (R3BTofdCalData const*)fCalTriggerItems->At(bot_trig_i); top_trig_ns = top_trig->GetTimeLeading_ns(); bot_trig_ns = bot_trig->GetTimeLeading_ns(); ++n1; @@ -265,8 +317,11 @@ void R3BTofDCal2HitPar::Exec(Option_t* option) // Shift the cyclic difference window by half a window-length and move it back, // this way the trigger time will be at 0. - auto top_ns = fTimeStitch->GetTime(top->GetTimeLeading_ns() - top_trig_ns); - auto bot_ns = fTimeStitch->GetTime(bot->GetTimeLeading_ns() - bot_trig_ns); + auto top_ns = + fmod(top->GetTimeLeading_ns() - top_trig_ns + c_range_ns + c_range_ns / 2, c_range_ns) - c_range_ns / 2; + auto bot_ns = + fmod(bot->GetTimeLeading_ns() - bot_trig_ns + c_range_ns + c_range_ns / 2, c_range_ns) - c_range_ns / 2; + auto dt = top_ns - bot_ns; // Handle wrap-around. auto dt_mod = fmod(dt + c_range_ns, c_range_ns); @@ -297,1099 +352,2668 @@ void R3BTofDCal2HitPar::Exec(Option_t* option) auto bot_tot = fmod(bot->GetTimeTrailing_ns() - bot->GetTimeLeading_ns() + c_range_ns, c_range_ns); // register multi hits multihits[iPlane - 1][iBar - 1] += 1; + + Double_t Ctop_tot = 0./0.; + Double_t Cbot_tot = 0./0.; - if (fParameter == 1) + Bool_t ToTexeption = false; + + if (fParameter == 1) { // calculate tdiff auto tdiff = (bot_ns) - (top_ns); - +/* + if(iPlane == 2 && iBar == 7){ + if(sqrt(top_tot*bot_tot) > fToTMin && sqrt(top_tot*bot_tot) < 500.) ToTexeption = true; + } + if(iPlane == 2 && iBar == 1){ + if(sqrt(top_tot*bot_tot)> fToTMin && sqrt(top_tot*bot_tot) < 350.) ToTexeption = true; + } + if(iPlane == 2 && iBar == 9){ + if(sqrt(top_tot*bot_tot)>fToTMin && sqrt(top_tot*bot_tot) < 500.) ToTexeption = true; + } +*/ // fill control histograms fhLogTot1vsLogTot2[iPlane - 1][iBar - 1]->Fill(TMath::Log(top_tot), TMath::Log(bot_tot)); - fh_tofd_TotPm[iPlane - 1]->Fill(iBar, top_tot); - fh_tofd_TotPm[iPlane - 1]->Fill(-iBar - 1, bot_tot); + if((sqrt(top_tot*bot_tot)>fToTMin && sqrt(top_tot*bot_tot)Fill(iBar, top_tot); + fh_tofd_TotPm[iPlane - 1][0]->Fill(-iBar - 1, bot_tot); + fh_tofd_TotPm1[iPlane - 1]->Fill(iBar, top_tot); + fh_tofd_TotPm1[iPlane - 1]->Fill(-iBar - 1, bot_tot); + } + if(/*top_tot>10. && */top_tot<90.){ + fh_tofd_TotPm2[iPlane - 1]->Fill(iBar, top_tot); + cout<<"Fill Top ToT "<10. && */bot_tot<90.){ + fh_tofd_TotPm2[iPlane - 1]->Fill(-iBar - 1, bot_tot); + cout<<"Fill Bot ToT "<Fill(iBar, tdiff); - + if((sqrt(top_tot*bot_tot)>fToTMin &&sqrt(top_tot*bot_tot)Fill(iBar, tdiff); + } auto posToT = TMath::Log(top_tot / bot_tot); - fhSqrtQvsPosToT[iPlane - 1][iBar - 1]->Fill(posToT, sqrt(top_tot * bot_tot)); - + if((sqrt(top_tot*bot_tot)>fToTMin &&sqrt(top_tot*bot_tot)Fill(posToT, sqrt(top_tot * bot_tot)); + } // Time of hit auto THit = (top_ns + bot_ns) / 2.; - fhTsync[iPlane - 1]->Fill(iBar, THit); - - // Tof with respect LOS detector - auto tof = fTimeStitch->GetTime((top_ns + bot_ns) / 2. - fHeader->GetTStart(), "tamex", "vftx"); - fh1_tofsync[iPlane - 1][iBar - 1]->Fill(tof); - // std::cout << "top" << top_ns << " bot"<GetTStart() << std::endl; - } - - // prepare offset and sync calculation - if (fTofdQ == 0.0 && fParameter > 1) - { - LOG(debug) << "Fill histo for offset and sync calculation plane " << iPlane << " bar " << iBar; - auto par = fHitPar->GetModuleParAt(iPlane, iBar); - if (!par) - { - R3BLOG(error, - "Hit par not found, Plane: " << top->GetDetectorId() << ", Bar: " << top->GetBarId()); - continue; - } - // calculate tdiff - auto tdiff = (bot_ns + par->GetOffset1()) - (top_ns + par->GetOffset2()); - - // fill control histograms - // fhTot1vsTot2[iPlane - 1][iBar - 1]->Fill(top_tot, bot_tot); - fhLogTot1vsLogTot2[iPlane - 1][iBar - 1]->Fill(TMath::Log(top_tot), TMath::Log(bot_tot)); - fh_tofd_TotPm[iPlane - 1]->Fill(iBar, top_tot); - fh_tofd_TotPm[iPlane - 1]->Fill(-iBar - 1, bot_tot); - - // Time differences of one paddle; offset histo - fhTdiff[iPlane - 1]->Fill(iBar, tdiff); - - // offset histo via ToT - /*if (fTofdY == 0.) - { - R3BLOG(debug, "Will prepare for offset and sync calculation"); - posToT = TMath::Log(top_tot / bot_tot); - } - else - {*/ - auto posToT = TMath::Log((top_tot + par->GetToTOffset2()) / (bot_tot + par->GetToTOffset1())); - //} - fhSqrtQvsPosToT[iPlane - 1][iBar - 1]->Fill(posToT, sqrt(top_tot * bot_tot)); - // fhTot1vsPos[iPlane - 1][iBar - 1]->Fill(posToT, bot_tot); - // fhTot2vsPos[iPlane - 1][iBar - 1]->Fill(posToT, top_tot); - - // Time of hit - auto THit = (top_ns + bot_ns) / 2. - par->GetSync(); - fhTsync[iPlane - 1]->Fill(iBar, THit); - - // Tof with respect LOS detector - auto tof = fTimeStitch->GetTime((top_ns + bot_ns) / 2. - fHeader->GetTStart(), "tamex", "vftx"); - fh1_tofsync[iPlane - 1][iBar - 1]->Fill(tof - par->GetTofSyncOffset()); - } - else if (fTofdQ > 0 && fParameter > 1) - { - // get parameter - auto para = fHitPar->GetModuleParAt(iPlane, iBar); - if (!para) - { - R3BLOG(error, - "Hit par not found, Plane: " << top->GetDetectorId() << ", Bar: " << top->GetBarId()); - continue; - } - - // calculate tdiff with offest - auto tdiff = (bot_ns + para->GetOffset1()) - (top_ns + para->GetOffset2()); - /* - // walk corrections - if (para->GetPar1Walk() == 0. || para->GetPar2Walk() == 0. || - para->GetPar3Walk() == 0. || para->GetPar4Walk() == 0. || para->GetPar5Walk() == 0.){ - R3BLOG(warn, "Walk correction not found!"); - }else{ - - bot_ns = bot_ns - walk(bot_tot, - para->GetPar1Walk(), - para->GetPar2Walk(), - para->GetPar3Walk(), - para->GetPar4Walk(), - para->GetPar5Walk()); - top_ns = top_ns - walk(top_tot, - para->GetPar1Walk(), - para->GetPar2Walk(), - para->GetPar3Walk(), - para->GetPar4Walk(), - para->GetPar5Walk()); - }*/ - - auto posToT = - para->GetLambda() * log((top_tot + para->GetToTOffset2()) / (bot_tot + para->GetToTOffset1())); - - // fill control histograms - // fhTot1vsTot2[iPlane - 1][iBar - 1]->Fill(top_tot, bot_tot); - fhSqrtQvsPosToT[iPlane - 1][iBar - 1]->Fill(posToT, sqrt(top_tot * bot_tot)); - - // Time differences of one paddle - fhTdiff[iPlane - 1]->Fill(iBar, tdiff); - - // Time of hit - auto THit = (bot_ns + top_ns) / 2. - para->GetSync(); - - // Sync of one plane - fhTsync[iPlane - 1]->Fill(iBar, THit); - } - - // prepare double exponential fit - if (!fTofdSmiley && fTofdQ > 0.1 && fParameter > 1) - { - LOG(debug) << "Prepare histo for double exponential fit"; - auto par = fHitPar->GetModuleParAt(iPlane, iBar); - if (!par) - { - R3BLOG(error, - "Hit par not found, Plane: " << top->GetDetectorId() << ", Bar: " << top->GetBarId()); - continue; - } - - // calculate y position - auto pos = ((bot_ns + par->GetOffset1()) - (top_ns + par->GetOffset2())) * par->GetVeff(); - - // fill fitting histograms and smiley histogram - fhTot1vsPos[iPlane - 1][iBar - 1]->Fill(pos, bot_tot); - fhTot2vsPos[iPlane - 1][iBar - 1]->Fill(pos, top_tot); - } - - // std::cout<< " hola 10" << fTofdZ << std::endl; - - // prepare charge fit / quench correction - if (fTofdZ == true && fParameter > 1) - { - // std::cout<< " hola 20" << std::endl; - LOG(debug) << "Prepare histo for quenching correction"; - // get parameter - auto par = fHitPar->GetModuleParAt(iPlane, iBar); - if (!par) - { - R3BLOG(error, - "Hit par not found, Plane: " << top->GetDetectorId() << ", Bar: " << top->GetBarId()); - continue; - } - - // Calculate y position - auto posToT = - par->GetLambda() * log((top_tot + par->GetToTOffset2()) / (bot_tot + par->GetToTOffset1())); - - Double_t parq[4]; - parq[0] = par->GetPar1a(); - parq[1] = par->GetPar1b(); - parq[2] = par->GetPar1c(); - parq[3] = par->GetPar1d(); - - // Calculate charge - Double_t qb = 0.; - if (fTofdSmiley) - { - qb = TMath::Sqrt(top_tot * bot_tot); // - //(parq[0] + parq[1] * posToT + parq[2] * pow(posToT, 2) + - // parq[3] * pow(posToT, 3)); - qb = qb * - fTofdQ; // theory says: dE ~ Z^2 but we see quenching -> just use linear and fit the rest - // std::cout<< " hola30 " << qb << " , " << fTofdQ << std::endl; - } - else - { - // double exponential - auto pos = ((bot_ns + par->GetOffset1()) - (top_ns + par->GetOffset2())) * par->GetVeff(); - auto q1 = bot_tot / - (parq[0] * (exp(-parq[1] * (pos + 100.)) + exp(-parq[2] * (pos + 100.))) + parq[3]); - parq[0] = par->GetPar2a(); - parq[1] = par->GetPar2b(); - parq[2] = par->GetPar2c(); - parq[3] = par->GetPar2d(); - auto q2 = top_tot / - (parq[0] * (exp(-parq[1] * (pos + 100.)) + exp(-parq[2] * (pos + 100.))) + parq[3]); - q1 = q1 * - fTofdQ; // theory says: dE ~ Z^2 but we see quenching -> just use linear and fit the rest - q2 = q2 * fTofdQ; - qb = (q1 + q2) / 2.; - } - - // Fill control histograms and Q vs Pos without multihits - if (multihits[iPlane - 1][iBar - 1] < 2 && (qb > 0.)) - { - // std::cout<< " h: " << posToT << " , " << qb << std::endl; - fhQvsPos[iPlane - 1][iBar - 1]->Fill(posToT, qb); - } - } + if((sqrt(top_tot*bot_tot)>fToTMin &&sqrt(top_tot*bot_tot)Fill(iBar, THit); + } + // Tof with respect LOS detector + auto tof = fTimeStitch->GetTime((top_ns + bot_ns) / 2. - fHeader->GetTStart(), "tamex", "vftx"); + if((sqrt(top_tot*bot_tot)>fToTMin &&sqrt(top_tot*bot_tot)Fill(tof); + fhTofsync[iPlane-1]->Fill(iBar,tof); + //cout<<"top_ns--> "<fToTMin && sqrt(top_tot*bot_tot)Fill(top_tot,top_ns - fHeader->GetTStart()); + fh1_walk_b[iPlane-1][iBar-1]->Fill(bot_tot,bot_ns - fHeader->GetTStart()); + + } + //auto lost = fHeader->GetTStart(); + // if(tof>-153.){ + // cout<<"sqrtToT"<Fill(iBar,lost); + // if(iPlane==1 && iBar==11){ + //fhlost[1][11]->Fill(THit,lost);} + //fhtott[iPlane-1][iBar-1]->Fill(THit,top_tot); + } + + if(fTofdQ == /*0.0*/ 8. && fParameter == 10){ + auto par = fHitPar->GetModuleParAt(iPlane, iBar); + if (!par) + { + R3BLOG(error, + "Hit par not found, Plane: " << top->GetDetectorId() << ", Bar: " << top->GetBarId()); + continue; + } + + //Ctop_tot = top_tot * par->GetToTOffset3() - par->GetToTOffset4(); + //Cbot_tot = bot_tot * par->GetToTOffset1() - par->GetToTOffset2(); + Ctop_tot = top_tot * par->GetToTOffset2(); + Cbot_tot = bot_tot * par->GetToTOffset1(); + + if(Ctop_tot > 180.){ + fh_tofd_TotPm[iPlane - 1][0]->Fill(iBar, Ctop_tot); + fh_tofd_TotPm1[iPlane - 1]->Fill(iBar, Ctop_tot); + cout<<"Fill Top ToT high"< 180.){ + fh_tofd_TotPm[iPlane - 1][0]->Fill(-iBar - 1, Cbot_tot); + fh_tofd_TotPm1[iPlane - 1]->Fill(-iBar - 1, Cbot_tot); + cout<<"Fill Bot ToT high"<10. && */Ctop_tot<70.){ + fh_tofd_TotPm2[iPlane - 1]->Fill(iBar, Ctop_tot); + cout<<"Fill Top ToT low"<10. && */Cbot_tot<70.){ + fh_tofd_TotPm2[iPlane - 1]->Fill(-iBar - 1, Cbot_tot); + cout<<"Fill Bot ToT low"<Fill(iBar, top_tot * par->GetToTOffset2()); + fh_tofd_TotPm[iPlane - 1][1]->Fill(-iBar - 1, bot_tot * par->GetToTOffset1()); + } + + // prepare offset and sync calculation + if (fTofdQ > 0. && fParameter > 1 && fParameter != 10) + { + + LOG(debug) << "Fill histo for offset and sync calculation plane " << iPlane << " bar " << iBar; + auto par = fHitPar->GetModuleParAt(iPlane, iBar); + if (!par) + { + R3BLOG(error, + "Hit par not found, Plane: " << top->GetDetectorId() << ", Bar: " << top->GetBarId()); + continue; + } + + // calculate tdiff + auto tdiff = (bot_ns + par->GetOffset1()) - (top_ns + par->GetOffset2()); + //Ctop_tot = top_tot * par->GetToTOffset3() - par->GetToTOffset4(); + //Cbot_tot = bot_tot * par->GetToTOffset1() - par->GetToTOffset2(); + Ctop_tot = top_tot * par->GetToTOffset2(); + Cbot_tot = bot_tot * par->GetToTOffset1(); + + // fill control histograms + // fhTot1vsTot2[iPlane - 1][iBar - 1]->Fill(top_tot, bot_tot); + fhLogTot1vsLogTot2[iPlane - 1][iBar - 1]->Fill(TMath::Log(top_tot), TMath::Log(bot_tot)); + fh_tofd_TotPm[iPlane - 1][0]->Fill(iBar, top_tot); + fh_tofd_TotPm[iPlane - 1][0]->Fill(-iBar - 1, bot_tot); + fh_tofd_TotPm[iPlane - 1][1]->Fill(iBar, top_tot * par->GetToTOffset2()); + fh_tofd_TotPm[iPlane - 1][1]->Fill(-iBar - 1, bot_tot * par->GetToTOffset1()); + fh_tofd_TotPm[iPlane - 1][2]->Fill(iBar, top_tot / par->GetToTOffset2()); + fh_tofd_TotPm[iPlane - 1][2]->Fill(-iBar - 1, bot_tot / par->GetToTOffset1()); + //cout<<"ToT_t "<GetToTOffset3()<<" par_b "<GetToTOffset4()<<" | ToT_b "<GetToTOffset1()<<" par_b "<GetToTOffset2()<fToTMin && sqrt(Ctop_tot*Cbot_tot)Fill(iBar, tdiff); + } + // offset histo via ToT + /*if (fTofdY == 0.) + { + R3BLOG(debug, "Will prepare for offset and sync calculation"); + posToT = TMath::Log(top_tot / bot_tot); + } + else + {*/ + auto posToT = TMath::Log((Ctop_tot) / (Cbot_tot)); + + //} + + cout<fToTMin && sqrt(Ctop_tot*Cbot_tot)Fill(posToT, sqrt(Ctop_tot * Cbot_tot)); + //fhTot1vsPos[iPlane - 1][iBar - 1]->Fill(posToT, Cbot_tot); + //fhTot2vsPos[iPlane - 1][iBar - 1]->Fill(posToT, Ctop_tot); + } + // Time of hit + if((sqrt(Ctop_tot*Cbot_tot)>fToTMin && sqrt(Ctop_tot*Cbot_tot)GetSync(); + fhTsync[iPlane - 1]->Fill(iBar, THit); + } + // Tof with respect LOS detector + if((sqrt(Ctop_tot*Cbot_tot)>fToTMin && sqrt(Ctop_tot*Cbot_tot)GetTime((top_ns + bot_ns) / 2. - fHeader->GetTStart(), "tamex", "vftx"); + if(fParameter == 2){ + fhTofsync[iPlane-1]->Fill(iBar,tof); + fh1_tofsync[iPlane - 1][iBar - 1]->Fill(tof); + } + else{ + fhTofsync[iPlane-1]->Fill(iBar,tof - par->GetTofSyncOffset()); + fh1_tofsync[iPlane - 1][iBar - 1]->Fill(tof - par->GetTofSyncOffset()); + } + } + } + if(fParameter==0){ + + if((sqrt(Ctop_tot*Cbot_tot)>fToTMin && sqrt(Ctop_tot*Cbot_tot)Fill(Ctop_tot,top_ns); + fh1_walk_b[iPlane-1][iBar-1]->Fill(Cbot_tot,bot_ns); + } + } + + + else if ((fTofdQ > 0 && fParameter > 1 && fParameter != 10) || fParameter == 3) + { + + // Double_t qb = 0.; + // get parameter + auto para = fHitPar->GetModuleParAt(iPlane, iBar); + if (!para) + { + R3BLOG(error, + "Hit par not found, Plane: " << top->GetDetectorId() << ", Bar: " << top->GetBarId()); + continue; + } + + // calculate tdiff with offest + auto tdiff = (bot_ns + para->GetOffset1()) - (top_ns + para->GetOffset2()); + //Ctop_tot = top_tot * para->GetToTOffset3() - para->GetToTOffset4(); + //Cbot_tot = bot_tot * para->GetToTOffset1() - para->GetToTOffset2(); + Ctop_tot = top_tot * para->GetToTOffset2(); + Cbot_tot = bot_tot * para->GetToTOffset1(); + + // walk corrections + if (para->GetPar1Walk() == 0. || para->GetPar2Walk() == 0. || + para->GetPar3Walk() == 0. || para->GetPar4Walk() == 0. || para->GetPar5Walk() == 0.){ + R3BLOG(warning, "Walk correction not found!"); + }else{ + + bot_ns = bot_ns - walk(Cbot_tot, + para->GetPar1Walk(), + para->GetPar2Walk(), + para->GetPar3Walk(), + para->GetPar4Walk(), + para->GetPar5Walk()); + top_ns = top_ns - walk(Ctop_tot, + para->GetPar1Walk(), + para->GetPar2Walk(), + para->GetPar3Walk(), + para->GetPar4Walk(), + para->GetPar5Walk()); + } + auto tof = fTimeStitch->GetTime((top_ns + bot_ns) / 2. - fHeader->GetTStart(), "tamex", "vftx"); + if((sqrt(top_tot*bot_tot)>fToTMin && sqrt(top_tot*bot_tot)Fill(iBar,tof - para->GetTofSyncOffset()); + // fh1_tofsync[iPlane - 1][iBar - 1]->Fill(tof - para->GetTofSyncOffset()); + } + + auto posToT = para->GetLambda() * log((Ctop_tot) / (Cbot_tot)); + if((sqrt(Ctop_tot*Cbot_tot)>fToTMin && sqrt(Ctop_tot*Cbot_tot)Fill(top_tot, bot_tot); + if((sqrt(Ctop_tot*Cbot_tot)>fToTMin && sqrt(Ctop_tot*Cbot_tot)Fill(posToT, sqrt(Ctop_tot * Cbot_tot)); + // std::cout<< " h0000000000: " << posToT << " , " << sqrt(Ctop_tot * Cbot_tot) << std::endl; + } + // Time differences of one paddle + if((sqrt(Ctop_tot*Cbot_tot)>fToTMin && sqrt(Ctop_tot*Cbot_tot)Fill(iBar, tdiff); + } + // Time of hit + auto THit = (bot_ns + top_ns) / 2. - para->GetSync(); + if((sqrt(Ctop_tot*Cbot_tot)>fToTMin && sqrt(Ctop_tot*Cbot_tot)fToTMin && sqrt(Ctop_tot*Cbot_tot)Fill(iBar, THit); + } + /* if (fTofdSmiley) + { + qb = TMath::Sqrt(top_tot * bot_tot); + // /(parq[0] + parq[1] * posToT + parq[2] * pow(posToT, 2) + parq[3] * pow(posToT, 3)); + + qb = qb * fTofdQ; + }*/ + + } + // prepare double exponential fit + if ((!fTofdSmiley && fTofdQ > 0.1 && fParameter > 1 && fParameter != 10)) + { + LOG(debug) << "Prepare histo for double exponential fit"; + auto par = fHitPar->GetModuleParAt(iPlane, iBar); + if (!par) + { + R3BLOG(error, + "Hit par not found, Plane: " << top->GetDetectorId() << ", Bar: " << top->GetBarId()); + continue; + } + //Ctop_tot = top_tot * par->GetToTOffset3() - par->GetToTOffset4(); + //Cbot_tot = bot_tot * par->GetToTOffset1() - par->GetToTOffset2(); + Ctop_tot = top_tot * par->GetToTOffset2(); + Cbot_tot = bot_tot * par->GetToTOffset1(); + + // calculate y position + auto pos = ((bot_ns + par->GetOffset1()) - (top_ns + par->GetOffset2())) * par->GetVeff(); + + // fill fitting histograms and smiley histogram + fhTot1vsPos[iPlane - 1][iBar - 1]->Fill(pos, Cbot_tot); + fhTot2vsPos[iPlane - 1][iBar - 1]->Fill(pos, Ctop_tot); + + } + + // std::cout<< " hola 10" << fTofdZ << std::endl; + + // prepare charge fit / quench correction + if (fTofdZ == true && fParameter > 1 && fParameter != 10) + { + + // std::cout<< " hola 20" << std::endl; + LOG(debug) << "Prepare histo for quenching correction"; + // get parameter + auto par = fHitPar->GetModuleParAt(iPlane, iBar); + if (!par) + { + R3BLOG(error, + "Hit par not found, Plane: " << top->GetDetectorId() << ", Bar: " << top->GetBarId()); + continue; + } + + // Calculate y position + auto posToT = + par->GetLambda() * log((Ctop_tot) / (Cbot_tot)); + + + Double_t parf[4]; + parf[0] = par->GetPar1a(); + parf[1] = par->GetPar1b(); + parf[2] = par->GetPar1c(); + parf[3] = par->GetPar1d(); + + + Double_t parq[4]; + parq[0] = par->GetPola(); + parq[1] = par->GetPolb(); + parq[2] = par->GetPolc(); + parq[3] = par->GetPold(); + + // Calculate charge + Double_t qb = 0.; + Double_t qbx = 0.; + if (fTofdSmiley) + { + if(parq[0]==0. && parq[1]==0. && parq[2]==0. && parq[3]==0.) parq[0] = 1.; + + qb = TMath::Sqrt(Ctop_tot * Cbot_tot);///(parq[0] + parq[1] * posToT + parq[2] * pow(posToT, 2) + parq[3] * pow(posToT, 3)); + qbx = TMath::Sqrt(Ctop_tot * Cbot_tot)*(parq[0] + parq[1] * posToT + parq[2] * pow(posToT, 2) + parq[3] * pow(posToT, 3)); + qb = qb * + fTofdQ / 200.; // theory says: dE ~ Z^2 but we see quenching -> just use linear and fit the rest + qbx = qbx * fTofdQ / 200.; + std::cout<< " hola30 " << qb << " , " << fTofdQ << " , "< 70.){ + qbx = qbx / 200.; + cout<<"XXXXXXX "<< qb << " , " << fTofdQ << " , "<GetOffset1()) - (top_ns + par->GetOffset2())) * par->GetVeff(); + auto q1 = Cbot_tot / + (parf[0] * (exp(-parf[1] * (pos + 100.)) + exp(-parf[2] * (pos + 100.))) + parf[3]*1.); + parq[0] = par->GetPar2a(); + parq[1] = par->GetPar2b(); + parq[2] = par->GetPar2c(); + parq[3] = par->GetPar2d(); + for(int i=0;i<4;i++){ + cout<<"parf["< just use linear and fit the rest + q2 = q2 * fTofdQ / 200.; + cout<<"q1 "<GetPar1za(); + parz[1] = par->GetPar1zb(); + parz[2] = par->GetPar1zc(); + + qbx = parz[0] + parz[1]*qbx + parz[2]*qbx*qbx; + } + + if (multihits[iPlane - 1][iBar - 1] < 2 /*&& qb>190. */ || 1) + { + fhQvsPos[iPlane - 1][iBar - 1]->Fill(posToT, qb); + fhQXvsPos[iPlane - 1][iBar - 1]->Fill(posToT, qbx); + } + + } + /* if(fParameter==5) + { + auto par = fHitPar->GetModuleParAt(iPlane, iBar); + + Double_t parz[3]; + parz[0] = par->GetPar1za(); + parz[1] = par->GetPar1zb(); + parz[2] = par->GetPar1zc(); + //parq[3] = par->GetPold(); + + + + Double_t qb = 0.; + if (fTofdSmiley) + { + qb = TMath::Sqrt(top_tot * bot_tot);///(parq[0] + parq[1] * posToT + parq[2] * pow(posToT, 2) + parq[3] * pow(posToT, 3)); + qb = qb * + fTofdQ; // theory says: dE ~ Z^2 but we see quenching -> just use linear and fit the rest + // std::cout<< " hola30 " << qb << " , " << fTofdQ << std::endl; + } + + // Fill control histograms and Q vs Pos without multihits + if (multihits[iPlane - 1][iBar - 1] < 2 && (qb > 5.6 && qb < 6.7.)) + { + // std::cout<< " h: " << posToT << " , " << qb << std::endl; + fhZoffset[iPlane - 1][iBar - 1]->Fill(parz[0]+parz[1]*qb); + } + + + + } +*/ + ++top_i; + ++bot_i; + + // Increment events + fNEvents += 1; + } + else if (dt < 0 && dt > -c_range_ns / 2) + { + ++top_i; + } + else + { + ++bot_i; + } + } + } + } + void R3BTofDCal2HitPar::SetHistogramsNULL(Int_t iPlane, Int_t iBar){ + fhTdiff[iPlane - 1] = NULL; + fhTsync[iPlane - 1] = NULL; + fhTofsync[iPlane - 1] = NULL; + fh_tofd_TotPm[iPlane - 1][0] = NULL; + fh_tofd_TotPm[iPlane - 1][1] = NULL; + fh_tofd_TotPm[iPlane - 1][2] = NULL; + fh_tofd_TotPm1[iPlane - 1] = NULL; + fh_tofd_TotPm2[iPlane - 1] = NULL; + fhLogTot1vsLogTot2[iPlane - 1][iBar - 1] = NULL; + fhSqrtQvsPosToT[iPlane - 1][iBar - 1] = NULL; + fhQvsPos[iPlane - 1][iBar - 1] = NULL; + fhQXvsPos[iPlane - 1][iBar - 1] = NULL; + fhTot1vsPos[iPlane - 1][iBar - 1] = NULL; + fhTot2vsPos[iPlane - 1][iBar - 1] = NULL; + fh1_tofsync[iPlane - 1][iBar - 1] = NULL; + fh1_walk_t[iPlane - 1][iBar - 1] = NULL; + fh1_walk_b[iPlane - 1][iBar - 1] = NULL; + + } + void R3BTofDCal2HitPar::CreateHistograms(Int_t iPlane, Int_t iBar) + { + Double_t max_charge = fMaxQ; + + if (NULL == fhTdiff[iPlane - 1]) + { + char *strName1; + strName1 = (char*)malloc(255); + char *strName2; + strName2 = (char*)malloc(255); + sprintf(strName1, "Time_Diff_Plane_%d", iPlane); + sprintf(strName2, "Time Diff Plane %d", iPlane); + fhTdiff[iPlane - 1] = new TH2F(strName1, strName2, 50, 0, 50, 4000, -20., 20.); + fhTdiff[iPlane - 1]->GetXaxis()->SetTitle("Bar "); + fhTdiff[iPlane - 1]->GetYaxis()->SetTitle("Time difference (PM1 - PM2) in ns"); + //cout<<"sprint1"<GetXaxis()->SetTitle("Bar "); + fhTsync[iPlane - 1]->GetYaxis()->SetTitle("THit in ns"); + //cout<<"sprint2"<GetXaxis()->SetTitle("Bar "); + fhTofsync[iPlane - 1]->GetYaxis()->SetTitle("Tof in ns"); + //cout<<"sprint3"<GetXaxis()->SetTitle("Bar "); + fhlos[iPlane - 1]->GetYaxis()->SetTitle("los t in ns"); + } + */ + if (NULL == fh_tofd_TotPm[iPlane - 1][0]) + { + char *strName; + strName = (char*)malloc(255); + char *strName2; + strName2 = (char*)malloc(255); + sprintf(strName, "Tofd_ToT_plane_%d_bc", iPlane); + sprintf(strName2, "Tofd ToT plane %d bc", iPlane); + fh_tofd_TotPm[iPlane - 1][0] = new TH2F(strName, strName2, 90, -45, 45, 500, 0., 500.); + fh_tofd_TotPm[iPlane - 1][0]->GetXaxis()->SetTitle("Bar "); + fh_tofd_TotPm[iPlane - 1][0]->GetYaxis()->SetTitle("ToT / ns"); + //cout<<"sprint4"<GetXaxis()->SetTitle("Bar "); + fh_tofd_TotPm1[iPlane - 1]->GetYaxis()->SetTitle("ToT / ns"); + //cout<<"sprint4"<GetXaxis()->SetTitle("Bar "); + fh_tofd_TotPm2[iPlane - 1]->GetYaxis()->SetTitle("ToT / ns"); + //cout<<"sprint4"<GetXaxis()->SetTitle("Bar "); + fh_tofd_TotPm[iPlane - 1][1]->GetYaxis()->SetTitle("ToT / ns"); + //cout<<"sprint4"<GetXaxis()->SetTitle("Bar "); + fh_tofd_TotPm[iPlane - 1][2]->GetYaxis()->SetTitle("ToT / ns"); + //cout<<"sprint4"<GetXaxis()->SetTitle("Log(ToT) of PM2"); + fhLogTot1vsLogTot2[iPlane - 1][iBar - 1]->GetYaxis()->SetTitle("Log(ToT) of PM1"); + free(strName); + free(strName2); + } + if (NULL == fhSqrtQvsPosToT[iPlane - 1][iBar - 1]) + { + char *strName; + strName = (char*)malloc(255); + char *strName2; + strName2 = (char*)malloc(255); + sprintf(strName, "SqrtQ_vs_PosToT_Plane_%d_Bar_%d", iPlane, iBar); + sprintf(strName2, "SqrtQ vs PosToT Plane %d Bar %d", iPlane, iBar); + //cout<<"sprint6 "<GetYaxis()->SetTitle("sqrt(PM1*PM2)"); + fhSqrtQvsPosToT[iPlane - 1][iBar - 1]->GetXaxis()->SetTitle("Position from ToT in cm"); + free(strName); + free(strName2); + } + if (NULL == fhQvsPos[iPlane - 1][iBar - 1]) + { + char *strName; + strName = (char*)malloc(255); + char *strName2; + strName2 = (char*)malloc(255); + sprintf(strName, "Q_vs_Pos_Plane_%d_Bar_%d", iPlane, iBar); + sprintf(strName2, "Q vs Pos Plane %d Bar %d", iPlane, iBar); + //cout<<"sprint7 "<GetYaxis()->SetTitle("Charge"); + fhQvsPos[iPlane - 1][iBar - 1]->GetXaxis()->SetTitle("Position in cm"); + free(strName); + free(strName2); + } + if (NULL == fhQXvsPos[iPlane - 1][iBar - 1]) + { + char *strName; + strName = (char*)malloc(255); + char *strName2; + strName2 = (char*)malloc(255); + sprintf(strName, "QX_vs_Pos_Plane_%d_Bar_%d", iPlane, iBar); + sprintf(strName2, "QX vs Pos Plane %d Bar %d", iPlane, iBar); + //cout<<"sprint7 "<GetYaxis()->SetTitle("Charge"); + fhQXvsPos[iPlane - 1][iBar - 1]->GetXaxis()->SetTitle("Position in cm"); + free(strName); + free(strName2); + } + /* if (NULL == fhZoffset[iPlane - 1][iBar - 1]) + { + char strName[255]; + sprintf(strName, "Zoffset_Plane_%d_Bar_%d", iPlane, iBar); + fhZoffset[iPlane - 1][iBar - 1] = new TH2F(strName, "", 500,0,10); + fhZoffset[iPlane - 1][iBar - 1]->GetYaxis()->SetTitle("counts"); + fhZoffset[iPlane - 1][iBar - 1]->GetXaxis()->SetTitle("Z"); + } + + if (NULL == fhlost[i][j]) + { + char strName[255]; + sprintf(strName, "los_vs_thit_Plane_%d_Bar_%d", 1, 11); + fhlost[i][j] = new TH2F(strName, "", 20000, -4000, -3000, 50000, -5000, 5000); + fhlost[i][j]->GetYaxis()->SetTitle("thit"); + fhlost[i][j]->GetXaxis()->SetTitle("los"); + } + */ + /* if (NULL == fhtott[iPlane - 1][iBar - 1]) + { + char strName[255]; + sprintf(strName, "tot_vs_thit_Plane_%d_Bar_%d", iPlane, iBar); + fhtott[iPlane - 1][iBar - 1] = new TH2F(strName, "", 20000, -4000, -2000, 1000, 0., 500); + fhtott[iPlane - 1][iBar - 1]->GetYaxis()->SetTitle("thit"); + fhtott[iPlane - 1][iBar - 1]->GetXaxis()->SetTitle("tot"); + }*/ + if (NULL == fhTot1vsPos[iPlane - 1][iBar - 1] /*&& !fTofdSmiley*/) + { + char *strName; + strName = (char*)malloc(255); + char *strName2; + strName2 = (char*)malloc(255); + sprintf(strName, "Tot1_vs_Pos_Plane_%d_Bar_%d", iPlane, iBar); + sprintf(strName2, "Tot1 vs Pos Plane %d Bar %d", iPlane, iBar); + fhTot1vsPos[iPlane - 1][iBar - 1] = new TH2F(strName, strName2, 300, -100, 100, 400, 0., 300.); + fhTot1vsPos[iPlane - 1][iBar - 1]->GetXaxis()->SetTitle("Pos in cm"); + fhTot1vsPos[iPlane - 1][iBar - 1]->GetYaxis()->SetTitle("ToT of PM1 in ns"); + //cout<<"sprint8"<GetXaxis()->SetTitle("Pos in cm"); + fhTot2vsPos[iPlane - 1][iBar - 1]->GetYaxis()->SetTitle("ToT of PM2 in ns"); + //cout<<"sprint9"<GetXaxis()->SetTitle("ToF [ns]"); + //cout<<"sprint10"<GetYaxis()->SetTitle("PMT time ns"); + fh1_walk_t[iPlane - 1][iBar - 1]->GetXaxis()->SetTitle("ToT of PMT"); + //cout<<"sprint11"<GetYaxis()->SetTitle("PMT time ns"); + fh1_walk_b[iPlane - 1][iBar - 1]->GetXaxis()->SetTitle("ToT of PMT"); + //cout<<"sprint12"< -c_range_ns / 2) - { - ++top_i; - } - else - { - ++bot_i; - } - } } -} -void R3BTofDCal2HitPar::CreateHistograms(Int_t iPlane, Int_t iBar) -{ - Double_t max_charge = fMaxQ; - if (NULL == fhTdiff[iPlane - 1]) + void R3BTofDCal2HitPar::calcOffset() { - char strName1[255]; - char strName2[255]; - sprintf(strName1, "Time_Diff_Plane_%d", iPlane); - sprintf(strName2, "Time Diff Plane %d", iPlane); - fhTdiff[iPlane - 1] = new TH2F(strName1, strName2, 50, 0, 50, 4000, -20., 20.); - fhTdiff[iPlane - 1]->GetXaxis()->SetTitle("Bar "); - fhTdiff[iPlane - 1]->GetYaxis()->SetTitle("Time difference (PM1 - PM2) in ns"); + TCanvas* cOffset = new TCanvas("cOffset", "cOffset", 10, 10, 1000, 900); + cOffset->Divide(2, 2); + R3BTofDHitModulePar* mpar; + for (Int_t i = 0; i < fNofPlanes; i++) + { + if (fhTdiff[i]) + { + LOG(warning) << "Found histo Time_Diff_Plane_" << i + 1; + // auto* h = (TH2F*)fhTdiff[i]->Clone(); + for (Int_t j = 0; j < fPaddlesPerPlane; j++) + { + int s=0; + mpar = new R3BTofDHitModulePar(); + Double_t offset = 0.; + cOffset->cd(i + 1); + fhTdiff[i]->Draw("colz"); + TH1F* histo_py = (TH1F*)fhTdiff[i]->ProjectionY("histo_py", j + 2, j + 2, ""); + histo_py->Rebin(4); + Int_t binmax = histo_py->GetMaximumBin(); + Double_t Max = histo_py->GetXaxis()->GetBinCenter(binmax); + TF1* fgaus = new TF1("fgaus", "gaus(0)", Max - 0.3, Max + 0.3); + histo_py->Fit("fgaus", "QR0"); + offset = fgaus->GetParameter(1); // histo_py->GetXaxis()->GetBinCenter(binmax); + LOG(warning) << " Plane " << i + 1 << " Bar " << j + 1 << " Offset " << offset; + mpar->SetPlane(i + 1); + mpar->SetPaddle(j + 1); + mpar->SetOffset1(-offset / 2.); + mpar->SetOffset2(offset / 2.); + fHitPar->AddModulePar(mpar); + //cin>>s; + } + fhTdiff[i]->Write(); + } + } + fHitPar->setChanged(); } - if (NULL == fhTsync[iPlane - 1]) + void R3BTofDCal2HitPar::calcToTOffset(Double_t totLow, Double_t totHigh) { - char strName[255]; - sprintf(strName, "Time_Sync_Plane_%d", iPlane); - fhTsync[iPlane - 1] = new TH2F(strName, strName, 45, 0, 45, 10000, -3500, 1000.); - fhTsync[iPlane - 1]->GetXaxis()->SetTitle("Bar "); - fhTsync[iPlane - 1]->GetYaxis()->SetTitle("THit in ns"); + TCanvas* cToTOffset1 = new TCanvas("cToTOffset1", "cToTOffset1", 10, 10, 1000, 900); + cToTOffset1->Divide(1, 2); + TCanvas* cToTOffset2 = new TCanvas("cToTOffset2", "cToTOffset2", 10, 10, 1000, 900); + cToTOffset2->Divide(1, 2); + Double_t fitrange = 25.; + if(fParameter == -1) fitrange = 10.; + for (Int_t i = 0; i < fNofPlanes; i++) + { + R3BTofDHitModulePar* par; + for (Int_t j = 0; j < (fPaddlesPerPlane); j++) + { + if(j==45 || j == 46) continue; + Double_t offset_bot_high = 0.; + Double_t offset_top_high = 0.; + Double_t offset_bot_low = 0.; + Double_t offset_top_low = 0.; + Double_t offset = 0.; + Double_t a_bot = 1.; + Double_t a_top = 1.; + Double_t b_bot = 0.; + Double_t b_top = 0.; + Double_t Max1 = 1.; + Double_t Max2 = 1.; + Double_t Max3 = 0.; + Double_t Max4 = 0.; + + if(j<44) par = fHitPar->GetModuleParAt(i + 1, 1 + j); + else par = fHitPar->GetModuleParAt(i + 1, j - fPaddlesPerPlane - 2); + if (fhSqrtQvsPosToT[i][j]) + { + LOG(info) << "Found histo SqrtQ_vs_PosToT_Plane_" << i + 1 << "_Bar_" << j + 1; + // auto* h = fhSqrtQvsPosToT[i][j]->Clone(); + cToTOffset1->cd(1); + fhSqrtQvsPosToT[i][j]->Draw("colz"); + auto histo_py = (TH2F*)fhSqrtQvsPosToT[i][j]->ProjectionX("histo_py", totLow, totHigh, ""); + cToTOffset1->cd(2); + histo_py->Rebin(2); + histo_py->Draw(); + Int_t binmax = histo_py->GetMaximumBin(); + Double_t Max = histo_py->GetXaxis()->GetBinCenter(binmax); + TF1* fgaus = new TF1( + "fgaus", "gaus(0)", Max - 0.2, Max + 0.2); // new TF1("fgaus", "gaus(0)", Max - 0.06, Max + 0.06); + histo_py->Fit("fgaus", "QR0"); + offset = fgaus->GetParameter(1); + fgaus->Draw("SAME"); + histo_py->SetAxisRange(Max - .5, Max + .5, "X"); + fhSqrtQvsPosToT[i][j]->SetAxisRange(Max - .5, Max + .5, "X"); + fhSqrtQvsPosToT[i][j]->SetAxisRange(totLow, totHigh, "Y"); + cToTOffset1->Update(); + delete fgaus; + fhSqrtQvsPosToT[i][j]->Write(); + } + LOG(warning) << " Plane " << i + 1 << " Bar " << j + 1 << " ToT Offset " << offset << "\n"; + LOG(info) <<"TOP " << 1. / sqrt(exp(offset)) << " BOT " <SetToTOffset1(sqrt(exp(offset))); + par->SetToTOffset2(1. / sqrt(exp(offset))); + + +/* if (fh_tofd_TotPm1[i]) + { + cToTOffset1->cd(1); + fh_tofd_TotPm1[i]->Draw("colz"); + auto histo_py1 = (TH1F*)fh_tofd_TotPm1[i]->ProjectionY("histo_py1", fPaddlesPerPlane - j, fPaddlesPerPlane - j, "");//bot high max + auto histo_py2 = (TH1F*)fh_tofd_TotPm1[i]->ProjectionY("histo_py2", fPaddlesPerPlane+j+3,fPaddlesPerPlane + j+3, "");//top high max + cToTOffset1->cd(2); + histo_py1->Rebin(2); + histo_py1->Draw(); + Int_t binmax1 = histo_py1->GetMaximumBin(); + Max1 = histo_py1->GetXaxis()->GetBinCenter(binmax1); + Int_t binmax2 = histo_py2->GetMaximumBin(); + Max2 = histo_py2->GetXaxis()->GetBinCenter(binmax2); + TF1* fgaus1 = new TF1( + "fgaus", "gaus(0)", Max1 - fitrange, Max1 + fitrange); // new TF1("fgaus", "gaus(0)", Max - 0.06, Max + 0.06); + histo_py1->Fit("fgaus", "QR0"); + offset_bot_high = fgaus1->GetParameter(1); + fgaus1->Draw("SAME"); + histo_py1->SetAxisRange(0., 350., "Y"); + cToTOffset1->Update(); + delete fgaus1; + histo_py2->Rebin(2); + histo_py2->Draw(); + TF1* fgaus2 = new TF1( + "fgaus", "gaus(0)", Max2 - fitrange, Max2 + fitrange); // new TF1("fgaus", "gaus(0)", Max - 0.06, Max + 0.06); + histo_py2->Fit("fgaus", "QR0"); + offset_top_high = fgaus2->GetParameter(1); + fgaus2->Draw("SAME"); + histo_py2->SetAxisRange(0., 350., "Y"); + cToTOffset1->Update(); + delete fgaus2; + } + if (fh_tofd_TotPm2[i]) + { + cToTOffset2->cd(1); + fh_tofd_TotPm2[i]->Draw("colz"); + auto histo_py3 = (TH1F*)fh_tofd_TotPm2[i]->ProjectionY("histo_py3", fPaddlesPerPlane - j, fPaddlesPerPlane - j, "");//bot low max + auto histo_py4 = (TH1F*)fh_tofd_TotPm2[i]->ProjectionY("histo_py4", fPaddlesPerPlane+j+3,fPaddlesPerPlane + j+3, "");//top low max + cToTOffset2->cd(2); + Int_t binmax3 = histo_py3->GetMaximumBin(); + Max3 = histo_py3->GetXaxis()->GetBinCenter(binmax3); + Int_t binmax4 = histo_py4->GetMaximumBin(); + Max4 = histo_py4->GetXaxis()->GetBinCenter(binmax4); + histo_py3->Rebin(2); + histo_py3->Draw(); + TF1* fgaus3 = new TF1( + "fgaus", "gaus(0)", Max3 - 25., Max3 + 25.); // new TF1("fgaus", "gaus(0)", Max - 0.06, Max + 0.06); + histo_py3->Fit("fgaus", "QR0"); + offset_bot_low = fgaus3->GetParameter(1); + fgaus3->Draw("SAME"); + histo_py3->SetAxisRange(0., 100., "Y"); + cToTOffset2->Update(); + delete fgaus3; + histo_py4->Rebin(2); + histo_py4->Draw(); + TF1* fgaus4 = new TF1( + "fgaus", "gaus(0)", Max4 - 25., Max4 + 25.); // new TF1("fgaus", "gaus(0)", Max - 0.06, Max + 0.06); + histo_py4->Fit("fgaus", "QR0"); + offset_top_low = fgaus4->GetParameter(1); + fgaus4->Draw("SAME"); + histo_py4->SetAxisRange(0., 100., "Y"); + cToTOffset2->Update(); + delete fgaus4; + + } + + + if((offset_bot_high != 0)){ + //a_bot = 160. / (offset_bot_high - offset_bot_low); + //b_bot = a_bot * offset_bot_high - 200.; + a_bot = 200. / offset_bot_high; + } + if(offset_top_high != 0){ + //a_top = 160. / (offset_top_high - offset_top_low); + //b_top = a_top * offset_top_high - 200.; + a_top = 200. / offset_top_high; + } + if(fParameter == 10){ + a_bot = 160. / (Max1 - Max3); + b_bot = a_bot * Max1 - 200.; + a_top = 160. / (Max2 - Max4); + b_top = a_top * Max2 - 200.; + + //a_bot = a_bot * par->GetToTOffset1(); + //b_bot = b_bot + a_bot * par->GetToTOffset2(); + //a_top = a_top * par->GetToTOffset3(); + //b_top = b_top + a_top * par->GetToTOffset4(); + } + cout<<"debugG "<<"P "<SetToTOffset1(a_bot);//Bottom factor + // par->SetToTOffset2(b_bot);//Bottom offset + par->SetToTOffset2(a_top);//Top factor + // par->SetToTOffset4(b_top);//Top offset + + */ + } + } + fHitPar->setChanged(); } - if (NULL == fh_tofd_TotPm[iPlane - 1]) - { - char strName[255]; - sprintf(strName, "Tofd_ToT_plane_%d", iPlane); - char strName2[255]; - sprintf(strName2, "Tofd ToT plane %d", iPlane); - fh_tofd_TotPm[iPlane - 1] = new TH2F(strName, strName2, 90, -45, 45, 300, 0., 300.); - fh_tofd_TotPm[iPlane - 1]->GetXaxis()->SetTitle("Bar "); - fh_tofd_TotPm[iPlane - 1]->GetYaxis()->SetTitle("ToT / ns"); - } - if (NULL == fhLogTot1vsLogTot2[iPlane - 1][iBar - 1]) + + void R3BTofDCal2HitPar::calcSync() { - char strName[255]; - sprintf(strName, "Plane_%d_Bar_%d_LogToT1vsLogToT2", iPlane, iBar); - fhLogTot1vsLogTot2[iPlane - 1][iBar - 1] = new TH2F(strName, "", 400, 2., 6., 400, 2., 6.); - fhLogTot1vsLogTot2[iPlane - 1][iBar - 1]->GetXaxis()->SetTitle("Log(ToT) of PM2"); - fhLogTot1vsLogTot2[iPlane - 1][iBar - 1]->GetYaxis()->SetTitle("Log(ToT) of PM1"); + TCanvas* cSync = new TCanvas("cSync", "cSync", 10, 10, 1000, 900); + cSync->Divide(2, 2); + for (Int_t i = 0; i < fNofPlanes; i++) + { + if (fhTsync[i]) + { + LOG(info) << "Found histo Time_Sync_Plane_" << i + 1; + // auto h = fhTsync[i]->Clone(); + for (Int_t j = 0; j < fPaddlesPerPlane; j++) + { + cSync->cd(i + 1); + fhTsync[i]->Draw("colz"); + auto* histo_py = (TH1F*)fhTsync[i]->ProjectionY("histo_py", j + 2, j + 2, ""); + R3BTofDHitModulePar* par = fHitPar->GetModuleParAt(i + 1, j + 1); + Int_t binmax = histo_py->GetMaximumBin(); + Double_t Max = histo_py->GetXaxis()->GetBinCenter(binmax); + Double_t MaxEntry = histo_py->GetBinContent(binmax); + TF1* fgaus = new TF1("fgaus", "gaus(0)", Max - 10., Max + 10.); + fgaus->SetParameters(MaxEntry, Max, 20); + histo_py->Fit("fgaus", "QR0"); + Double_t sync = fgaus->GetParameter(1); // histo_py->GetXaxis()->GetBinCenter(binmax); + par->SetSync(sync); + LOG(info) << " Plane " << i + 1 << " Bar " << j + 1 << " Sync " << sync; + } + fhTsync[i]->Write(); + } + } + fHitPar->setChanged(); } - if (NULL == fhSqrtQvsPosToT[iPlane - 1][iBar - 1]) + + int R3BTofDCal2HitPar::zcorr(TH2F* histo, Double_t min, Double_t max, Double_t* pars, Int_t pl, Int_t b, Int_t index) { - char strName[255]; - sprintf(strName, "SqrtQ_vs_PosToT_Plane_%d_Bar_%d", iPlane, iBar); - fhSqrtQvsPosToT[iPlane - 1][iBar - 1] = - new TH2F(strName, "", 2000, -100, 100, max_charge * 4, 0., max_charge * 4); - fhSqrtQvsPosToT[iPlane - 1][iBar - 1]->GetYaxis()->SetTitle("sqrt(PM1*PM2)"); - fhSqrtQvsPosToT[iPlane - 1][iBar - 1]->GetXaxis()->SetTitle("Position from ToT in cm"); + if (histo->GetEntries() < 10) + { + R3BLOG(warning, "Nb of events below 10 with "<< histo->GetEntries() << " entries for histo with index" << index); + + //return; + } + + Double_t par[3000] = { 0 }; + Double_t x[3000] = { 0 }; + char strName[255]; + sprintf(strName, "canvas_%d", index); + TCanvas* c1 = new TCanvas(strName, "", 100, 100, 800, 800); + c1->Divide(1, 3); + c1->cd(1); + auto* h = (TH2F*)histo->Clone(); + h->Draw("colz"); + h->SetAxisRange(min, max, "Y"); + // Projection of charge axis + auto* h1 = h->ProjectionY("p_y", 1, 2000); + + c1->cd(2); + //gPad->SetLogy(); + h1->Draw(); + // Use TSpectrum to find the peak candidates + Double_t threshhold = 0.001; + Int_t sigma = 10; + Double_t sqrtf = 1.; + Int_t nPeaks = fNbZPeaks; +doagainfit: + TSpectrum* s = new TSpectrum(nPeaks); + Int_t nfound = s->Search(h1, sigma, "nobackground", threshhold); + std::cout << "Found " << nfound << " candidate peaks to fit\n"; + + if (nfound == 0) + { + delete s; + delete c1; + delete h; + delete h1; + return -1; + } + + c1->Update(); + // Eliminate background peaks + Int_t mPeaks = 0; + Double_t* xpeaks = s->GetPositionX(); + for (Int_t p = 0; p <= nfound; p++) + { + Float_t xp = xpeaks[p]; + Int_t bin = h1->GetXaxis()->FindBin(xp); + Float_t yp = h1->GetBinContent(bin); + if (yp - sqrtf * TMath::Sqrt(yp) < 1.){ + cout<<"jump "<Update(); + + if (mPeaks < 2 && 0) + { + pars[0] = 0.; + //pars[1] = fMaxQ / peaks[0]; + pars[1] = 8. / peaks[0]; + pars[2] = 0.; + + delete s; + delete c1; + delete h; + delete h1; + return -1; + } + cout<<"mPeaks "<> peakscheck; + if(peakscheck == "n"){ + delete s; + std::cout << "New sigma: "; + std::cin >> sigma; + std::cout << "New threshhold: "; + std::cin >> threshhold; + std::cout << "New Sqrt factor: "; + std::cin >> sqrtf; + std::cout << "New # of peaks: "; + std::cin >> nPeaks; + goto doagainfit; + } + if(peakscheck == "goback"){ + int a=0; + std::cout << "How many bars? "; + std::cin >> a; + return a; + } + while(peakscheck == 'c'){ + cout<< "We see "<> peakscheck; + } + + nfp = 0; + for (Int_t i = 0; i < mPeaks; i++) + { + std::cout << "Peak @ " << peaks[i]; + while ((std::cout << " corresponds to Z=") && !(std::cin >> z)) + { + std::cout << "That's not a number;"; + std::cin.clear(); + std::cin.ignore(std::numeric_limits::max(), '\n'); + } + if (z == 0) + continue; + x[nfp] = (Double_t)z; + zpeaks[nfp] = peaks[i]; + // std::cout<<"z real " << x[nfp] << " z found " << zpeaks[nfp] <<"\n"; + nfp++; + } + doagain = 'n'; +// std::cout << "Do again? (y/n) "; +// std::cin >> doagain; + } while (doagain != "n"); + + if (mPeaks < 2) + { + pars[0] = 0.; + //pars[1] = fMaxQ / peaks[0]; + pars[1] = z / peaks[0]; + pars[2] = 0.; + + delete s; + delete c1; + delete h; + delete h1; + return -1; + } + + // fit charge axis + std::cout << "Selected " << nfp << " useful peaks to fit\nStart fitting...\n"; + c1->cd(3)->Clear(); + c1->Update(); + c1->cd(3); + auto gr1 = new TGraph(nfp, zpeaks, x); + gr1->SetMarkerColor(4); + gr1->SetMarkerSize(1.5); + gr1->SetMarkerStyle(20); + gr1->Draw("AP"); + + fZfitType = "pol2"; + + // TF1* fitz = new TF1("fitz", "[0]*TMath::Power(x,[2])+[1]", min, max); + if (fZfitType != "pol1" && fZfitType != "pol2") + { + R3BLOG(error, "Fit " << fZfitType << " is not allowed, use pol1 or pol2 "); + return -1; + } + auto fitz = new TF1("fitz", fZfitType, min, max); + fitz->SetLineColor(2); + fitz->SetLineWidth(2); + fitz->SetLineStyle(1); + // fitz->SetParameters(1.5, 2., .1); + gr1->Fit("fitz", "Q"); + fitz->Draw("lsame"); + c1->Update(); + // write parameters + std::cout << "Is OK? (y/n) "; + std::cin >> doagain; + if (doagain == "n"){ + std::cout << "New sigma: "; + std::cin >> sigma; + std::cout << "New threshhold: "; + std::cin >> threshhold; + delete s; + goto doagainfit; + } + + int nbpars = 2; + if (fZfitType == "pol2") + nbpars = 3; + + for (Int_t j = 0; j < nbpars; j++) + { + pars[j] = fitz->GetParameter(j); + std::cout<Write(); + delete s; + delete h; + delete h1; + delete gr1; + delete c1; + delete fitz; + + return -1; } - if (NULL == fhQvsPos[iPlane - 1][iBar - 1]) + + void R3BTofDCal2HitPar::calcVeff() { - char strName[255]; - sprintf(strName, "Q_vs_Pos_Plane_%d_Bar_%d", iPlane, iBar); - fhQvsPos[iPlane - 1][iBar - 1] = new TH2F(strName, "", 2000, -100, 100, max_charge / 2, 0., max_charge); - fhQvsPos[iPlane - 1][iBar - 1]->GetYaxis()->SetTitle("Charge"); - fhQvsPos[iPlane - 1][iBar - 1]->GetXaxis()->SetTitle("Position in cm"); + TCanvas* cVeff = new TCanvas("cVeff", "cVeff", 10, 10, 1000, 900); + cVeff->Divide(2, 2); + for (Int_t i = 0; i < fNofPlanes; i++) + { + for (Int_t j = 0; j < fPaddlesPerPlane; j++) + { + Double_t max = 0.; + Double_t veff = 7.; + if (fhTdiff[i]) + { + LOG(info) << "Found histo Time_Diff_Plane_" << i + 1; + auto par = fHitPar->GetModuleParAt(i + 1, j + 1); + if (!par) + { + LOG(warning) << "Hit par not found, Plane: " << i + 1 << ", Bar: " << j + 1; + continue; + } + // auto* h = (TH2F*)histofilename->Get(Form("Time_Diff_Plane_%i", i + 1))->Clone(); + cVeff->cd(i + 1); + // h->Draw("colz"); + TH1F* histo_py = (TH1F*)fhTdiff[i]->ProjectionY("histo_py", j + 2, j + 2, ""); + Int_t binmax = histo_py->GetMaximumBin(); + max = histo_py->GetXaxis()->GetBinCenter(binmax); + Double_t maxEntry = histo_py->GetBinContent(binmax); + auto* fgaus = new TF1("fgaus", "gaus(0)", max - 0.3, max + 0.3); /// TODO: find best range + fgaus->SetParameters(maxEntry, max, 20); + histo_py->Fit("fgaus", "QR0"); + Double_t offset1 = par->GetOffset1(); + Double_t offset2 = par->GetOffset2(); + max = fgaus->GetParameter(1) + offset1 - offset2; /// TODO: needs to be tested + // max = max+offset1-offset2; + veff = fTofdY / max; // effective speed of light in [cm/s] + LOG(info) << "Plane " << i + 1 << " Bar " << j + 1 << " offset " << par->GetOffset1(); + LOG(info) << "Plane " << i + 1 << " Bar " << j + 1 << " max " << max; + LOG(info) << "Plane " << i + 1 << " Bar " << j + 1 << " veff " << veff; + par->SetVeff(veff); + } + } + } + fHitPar->setChanged(); } - if (NULL == fhTot1vsPos[iPlane - 1][iBar - 1] && !fTofdSmiley) + void R3BTofDCal2HitPar::calcLambda(Double_t totLow, Double_t totHigh) { - char strName[255]; - sprintf(strName, "Tot1_vs_Pos_Plane_%d_Bar_%d", iPlane, iBar); - fhTot1vsPos[iPlane - 1][iBar - 1] = new TH2F(strName, "", 200, -100, 100, 400, 0., 200.); - fhTot1vsPos[iPlane - 1][iBar - 1]->GetXaxis()->SetTitle("Pos in cm"); - fhTot1vsPos[iPlane - 1][iBar - 1]->GetYaxis()->SetTitle("ToT of PM1 in ns"); + TCanvas* cToTOffset = new TCanvas("cLambda", "cLambda", 10, 10, 1000, 900); + cToTOffset->Divide(1, 2); + for (Int_t i = 0; i < fNofPlanes; i++) + { + for (Int_t j = 0; j < fPaddlesPerPlane; j++) + { + Double_t offset = 0.; + auto par = fHitPar->GetModuleParAt(i + 1, j + 1); + if (fhSqrtQvsPosToT[i][j]) + { + LOG(info) << "Found histo SqrtQ_vs_PosToT_Plane_" << i + 1 << "_Bar_" << j + 1; + // auto* h = (TH2F*)histofilename->Get(Form("SqrtQ_vs_PosToT_Plane_%i_Bar_%i", i + 1, j + 1))->Clone(); + cToTOffset->cd(1); + fhSqrtQvsPosToT[i][j]->Draw("colz"); + auto* histo_py = (TH2F*)fhSqrtQvsPosToT[i][j]->ProjectionX("histo_py", totLow, totHigh, ""); + cToTOffset->cd(2); + histo_py->Draw(); + Int_t binmax = histo_py->GetMaximumBin(); + Double_t Max = histo_py->GetXaxis()->GetBinCenter(binmax); + TF1* fgaus = new TF1("fgaus", "gaus(0)", Max - 0.6, Max + 0.6); + histo_py->Fit("fgaus", "QR0"); + offset = fgaus->GetParameter(1); + fgaus->Draw("SAME"); + histo_py->SetAxisRange(Max - .5, Max + .5, "X"); + fhSqrtQvsPosToT[i][j]->SetAxisRange(Max - .5, Max + .5, "X"); + fhSqrtQvsPosToT[i][j]->SetAxisRange(totLow, totHigh, "Y"); + cToTOffset->Update(); + delete fgaus; + delete histo_py; + } + else + LOG(error) << "Missing histo plane " << i + 1 << " bar " << j + 1; + Double_t lambda = fTofdY / offset; + LOG(info) << " Plane " << i + 1 << " Bar " << j + 1 << " ToT Offset " << offset << " Lambda " << lambda + << "\n"; + par->SetLambda(lambda); + } + } + fHitPar->setChanged(); } - if (NULL == fhTot2vsPos[iPlane - 1][iBar - 1] && !fTofdSmiley) + /* void R3BTofDCal2HitPar::calcwalk(Double_t totLow, Double_t totHigh) + { + for (Int_t i = 0; i < fNofPlanes; i++) + { + for (Int_t j = 0; j < fPaddlesPerPlane; j++) + { + auto par = fHitPar->GetModuleParAt(i + 1, j + 1); + Double_t y[1000], x[1000]; + Int_t n = 0; + + TGraph* gr1 = new TGraph(); + TGraph* gr2 = new TGraph(); + TCanvas* cfit_walk = new TCanvas("cfit_walk", "fit walk", 100, 100, 800, 800); + cfit_walk->Clear(); + cfit_walk->Divide(1, 4); + cfit_walk->cd(1); + TH2F* histowb1 = (TH2F*)histowb->Clone(); + histowb1->Draw("colz"); + TH2F* histowb2 = (TH2F*)histowb->Clone(); + histowb2->RebinX(50); + histowb2->GetYaxis()->SetRangeUser(fTofdTotLow, fTofdTotHigh); + // histo2->SetAxisRange(fTofdTotLow,fTofdTotHigh,"Y"); + cfit_walk->cd(2); + histowb2->Draw("colz"); + std::cout << "Searching for points to fit...\n"; + for (Int_t i = 1; i < histowb2->GetNbinsX(); i++) { - char strName[255]; - sprintf(strName, "Tot2_vs_Pos_Plane_%d_Bar_%d", iPlane, iBar); - fhTot2vsPos[iPlane - 1][iBar - 1] = new TH2F(strName, "", 200, -100, 100, 400, 0., 200.); - fhTot2vsPos[iPlane - 1][iBar - 1]->GetXaxis()->SetTitle("Pos in cm"); - fhTot2vsPos[iPlane - 1][iBar - 1]->GetYaxis()->SetTitle("ToT of PM2 in ns"); - } - if (fh1_tofsync[iPlane - 1][iBar - 1] == NULL) + std::cout<<"Bin "<GetNbinsX()<<" with cut: "<cd(2); + TLine* lw = new TLine( + histowb2->GetXaxis()->GetBinCenter(i), fTofdTotLow, histowb2->GetXaxis()->GetBinCenter(i), fTofdTotHigh); + lw->SetLineColor(kRed); + lw->SetLineWidth(2.); + lw->Draw(); + cfit_smiley->cd(3); + TH1F* histo_wby = (TH1F*)histowb2->ProjectionY("histo_wby", i, i, ""); + histo_wby->Draw(); + // cfit_smiley->Update(); + x[n] = histowb2->GetXaxis()->GetBinCenter(i); + Int_t binmaxwb = histo_wby->GetMaximumBin(); + y[n] = histo_wby->GetXaxis()->GetBinCenter(binmax); + + if ((x[n] < min || x[n] > max) || (y[n] < fTofdTotLow || y[n] > fTofdTotHigh)) { - char strName[255]; - sprintf(strName, "tofdiff_plane_%d_bar_%d", iPlane, iBar); - fh1_tofsync[iPlane - 1][iBar - 1] = new TH1F(strName, strName, 25000, 2950, 3150); - fh1_tofsync[iPlane - 1][iBar - 1]->GetXaxis()->SetTitle("ToF [ns]"); + delete histo_wby; + continue; } -} - -void R3BTofDCal2HitPar::calcOffset() -{ - TCanvas* cOffset = new TCanvas("cOffset", "cOffset", 10, 10, 1000, 900); - cOffset->Divide(2, 2); - R3BTofDHitModulePar* mpar; - for (Int_t i = 0; i < fNofPlanes; i++) + if (histo_wby->GetMaximum() > 5) { - if (fhTdiff[i]) - { - LOG(warn) << "Found histo Time_Diff_Plane_" << i + 1; - // auto* h = (TH2F*)fhTdiff[i]->Clone(); - for (Int_t j = 0; j < fPaddlesPerPlane; j++) - { - mpar = new R3BTofDHitModulePar(); - Double_t offset = 0.; - cOffset->cd(i + 1); - fhTdiff[i]->Draw("colz"); - TH1F* histo_py = (TH1F*)fhTdiff[i]->ProjectionY("histo_py", j + 2, j + 2, ""); - histo_py->Rebin(4); - Int_t binmax = histo_py->GetMaximumBin(); - Double_t Max = histo_py->GetXaxis()->GetBinCenter(binmax); - TF1* fgaus = new TF1("fgaus", "gaus(0)", Max - 0.3, Max + 0.3); - histo_py->Fit("fgaus", "QR0"); - offset = fgaus->GetParameter(1); // histo_py->GetXaxis()->GetBinCenter(binmax); - LOG(warn) << " Plane " << i + 1 << " Bar " << j + 1 << " Offset " << offset; - mpar->SetPlane(i + 1); - mpar->SetPaddle(j + 1); - mpar->SetOffset1(-offset / 2.); - mpar->SetOffset2(offset / 2.); - fHitPar->AddModulePar(mpar); - } - } + n++; + delete lw; } - fHitPar->setChanged(); -} -void R3BTofDCal2HitPar::calcToTOffset(Double_t totLow, Double_t totHigh) -{ - TCanvas* cToTOffset = new TCanvas("cToTOffset", "cToTOffset", 10, 10, 1000, 900); - cToTOffset->Divide(1, 2); - for (Int_t i = 0; i < fNofPlanes; i++) - { - for (Int_t j = 0; j < fPaddlesPerPlane; j++) - { - Double_t offset = 0.; - R3BTofDHitModulePar* par = fHitPar->GetModuleParAt(i + 1, j + 1); - if (fhSqrtQvsPosToT[i][j]) - { - LOG(info) << "Found histo SqrtQ_vs_PosToT_Plane_" << i + 1 << "_Bar_" << j + 1; - // auto* h = fhSqrtQvsPosToT[i][j]->Clone(); - cToTOffset->cd(1); - fhSqrtQvsPosToT[i][j]->Draw("colz"); - auto histo_py = (TH2F*)fhSqrtQvsPosToT[i][j]->ProjectionX("histo_py", totLow, totHigh, ""); - cToTOffset->cd(2); - histo_py->Rebin(2); - histo_py->Draw(); - Int_t binmax = histo_py->GetMaximumBin(); - Double_t Max = histo_py->GetXaxis()->GetBinCenter(binmax); - TF1* fgaus = new TF1( - "fgaus", "gaus(0)", Max - 0.2, Max + 0.2); // new TF1("fgaus", "gaus(0)", Max - 0.06, Max + 0.06); - histo_py->Fit("fgaus", "QR0"); - offset = fgaus->GetParameter(1); - fgaus->Draw("SAME"); - histo_py->SetAxisRange(Max - .5, Max + .5, "X"); - fhSqrtQvsPosToT[i][j]->SetAxisRange(Max - .5, Max + .5, "X"); - fhSqrtQvsPosToT[i][j]->SetAxisRange(totLow, totHigh, "Y"); - cToTOffset->Update(); - delete fgaus; - } - LOG(warn) << " Plane " << i + 1 << " Bar " << j + 1 << " ToT Offset " << offset << "\n"; - par->SetToTOffset1(sqrt(exp(offset))); - par->SetToTOffset2(1. / sqrt(exp(offset))); - } + delete histo_wby; } - fHitPar->setChanged(); -} + gr1 = new TGraph(n, x, y); + gr1->SetTitle("Points found for fitting; bot time in ns; Bot ToT"); + gr1->Draw("A*"); + std::cout << "Start fitting\n"; + //Double_t y = 0; + //y = -30.2 + par1 * TMath::Power(Q, par2) + par3 / Q + par4 * Q + par5 * Q * Q; + TF1* f1 = new TF1("f1", "", min, max); + f1->SetLineColor(2); + gr1->Fit("f1", "Q", -30.2 + par[1] * TMath::Power(x, par[2]) + par[3] / x + par[4] * x + par[5] * x * x; + , min, max); + //for (Int_t j = 0; j <= 3; j++) -void R3BTofDCal2HitPar::calcSync() -{ - TCanvas* cSync = new TCanvas("cSync", "cSync", 10, 10, 1000, 900); - cSync->Divide(2, 2); - for (Int_t i = 0; i < fNofPlanes; i++) - { - if (fhTsync[i]) - { - LOG(info) << "Found histo Time_Sync_Plane_" << i + 1; - // auto h = fhTsync[i]->Clone(); - for (Int_t j = 0; j < fPaddlesPerPlane; j++) - { - cSync->cd(i + 1); - fhTsync[i]->Draw("colz"); - auto* histo_py = (TH1F*)fhTsync[i]->ProjectionY("histo_py", j + 2, j + 2, ""); - R3BTofDHitModulePar* par = fHitPar->GetModuleParAt(i + 1, j + 1); - Int_t binmax = histo_py->GetMaximumBin(); - Double_t Max = histo_py->GetXaxis()->GetBinCenter(binmax); - Double_t MaxEntry = histo_py->GetBinContent(binmax); - TF1* fgaus = new TF1("fgaus", "gaus(0)", Max - 10., Max + 10.); - fgaus->SetParameters(MaxEntry, Max, 20); - histo_py->Fit("fgaus", "QR0"); - Double_t sync = fgaus->GetParameter(1); // histo_py->GetXaxis()->GetBinCenter(binmax); - par->SetSync(sync); - LOG(info) << " Plane " << i + 1 << " Bar " << j + 1 << " Sync " << sync; - } - } - } - fHitPar->setChanged(); -} + para1 = f1->GetParameter(par[1]); + para2 = f1->GetParameter(par[2]); + para3 = f1->GetParameter(par[3]); + para4 = f1->GetParameter(par[4]); + para5 = f1->GetParameter(par[5]); -void R3BTofDCal2HitPar::zcorr(TH2F* histo, Int_t min, Int_t max, Double_t* pars, Int_t index) -{ - if (histo->GetEntries() < 100) - { - R3BLOG(warn, "Nb of events below 100 for histo with index" << index); - return; - } - Double_t par[3000] = { 0 }; - Int_t nPeaks = fNbZPeaks; - Double_t x[3000] = { 0 }; - char strName[255]; - sprintf(strName, "canvas_%d", index); - TCanvas* c1 = new TCanvas(strName, "", 100, 100, 800, 800); - c1->Divide(1, 3); - c1->cd(1); - auto* h = dynamic_cast(histo->Clone()); - h->Draw("colz"); - h->SetAxisRange(min, max, "Y"); - // Projection of charge axis - auto* h1 = h->ProjectionY("p_y"); - c1->cd(2); - h1->Draw(); - // Use TSpectrum to find the peak candidates - TSpectrum* s = new TSpectrum(nPeaks); - Int_t nfound = s->Search(h1, 10, "", 0.001); - std::cout << "Found " << nfound << " candidate peaks to fit\n"; - - if (nfound == 0) - { - delete s; - delete c1; - delete h; - delete h1; - return; - } - c1->Update(); - // Eliminate background peaks - nPeaks = 0; - Double_t* xpeaks = s->GetPositionX(); - for (Int_t p = 0; p <= nfound; p++) - { - Float_t xp = xpeaks[p]; - Int_t bin = h1->GetXaxis()->FindBin(xp); - Float_t yp = h1->GetBinContent(bin); - if (yp - TMath::Sqrt(yp) < 1.) - continue; - par[2 * nPeaks] = yp; - par[2 * nPeaks + 1] = xp; - nPeaks++; - } - Double_t peaks[nPeaks]; - for (Int_t i = 0; i < nPeaks; i++) - { - // printf("Found peak @ %f\n",xpeaks[i]); - peaks[i] = par[2 * i + 1]; - } - c1->Update(); + std::cout << "Parameter: " << para << "\n"; + par->SetPar1Walk(para1); + par->SetPar2Walk(para2); + par->SetPar3Walk(para3); + par->SetPar4Walk(para4); + par->SetPar5Walk(para5); - if (nPeaks < 2) - { - pars[0] = 0.; - pars[1] = fTofdQ / peaks[0]; - pars[2] = 0.; - - delete s; - delete c1; - delete h; - delete h1; - return; - } - // select useful peaks - sort(peaks, peaks + nPeaks); - Double_t zpeaks[3000] = { 0 }; - string doagain = "y"; - Int_t nfp; -doagainfit: - do - { - nfp = 0; - for (Int_t i = 0; i < nPeaks; i++) - { - std::cout << "Peak @ " << peaks[i]; - Int_t z = 0; - while ((std::cout << " corresponds to Z=") && !(std::cin >> z)) - { - std::cout << "That's not a number;"; - std::cin.clear(); - std::cin.ignore(std::numeric_limits::max(), '\n'); - } - if (z == 0) - continue; - x[nfp] = (Double_t)z; - zpeaks[nfp] = peaks[i]; - // std::cout<<"z real " << x[nfp] << " z found " << zpeaks[nfp] <<"\n"; - nfp++; - } - std::cout << "Do again? (y/n) "; - std::cin >> doagain; - } while (doagain != "n"); - - // fit charge axis - std::cout << "Selected " << nfp << " useful peaks to fit\nStart fitting...\n"; - c1->cd(3)->Clear(); - c1->Update(); - c1->cd(3); - auto gr1 = new TGraph(nfp, zpeaks, x); - gr1->SetMarkerColor(4); - gr1->SetMarkerSize(1.5); - gr1->SetMarkerStyle(20); - gr1->Draw("AP"); - // TF1* fitz = new TF1("fitz", "[0]*TMath::Power(x,[2])+[1]", min, max); - if (fZfitType != "pol1" && fZfitType != "pol2") - { - R3BLOG(error, "Fit " << fZfitType << " is not allowed, use pol1 or pol2 "); - return; - } - auto fitz = new TF1("fitz", fZfitType, min, max); - fitz->SetLineColor(2); - fitz->SetLineWidth(2); - fitz->SetLineStyle(1); - // fitz->SetParameters(1.5, 2., .1); - gr1->Fit("fitz", "Q"); - fitz->Draw("lsame"); - c1->Update(); - // write parameters - std::cout << "Is OK? (y/n) "; - std::cin >> doagain; - if (doagain == "n") - goto doagainfit; - - int nbpars = 2; - if (fZfitType == "pol2") - nbpars = 3; - - for (Int_t j = 0; j < nbpars; j++) - { - pars[j] = fitz->GetParameter(j); - // std::cout<Divide(2, 2); - for (Int_t i = 0; i < fNofPlanes; i++) - { - for (Int_t j = 0; j < fPaddlesPerPlane; j++) - { - Double_t max = 0.; - Double_t veff = 7.; - if (fhTdiff[i]) - { - LOG(info) << "Found histo Time_Diff_Plane_" << i + 1; - auto par = fHitPar->GetModuleParAt(i + 1, j + 1); - if (!par) - { - LOG(warn) << "Hit par not found, Plane: " << i + 1 << ", Bar: " << j + 1; - continue; - } - // auto* h = (TH2F*)histofilename->Get(Form("Time_Diff_Plane_%i", i + 1))->Clone(); - cVeff->cd(i + 1); - // h->Draw("colz"); - TH1F* histo_py = (TH1F*)fhTdiff[i]->ProjectionY("histo_py", j + 2, j + 2, ""); - Int_t binmax = histo_py->GetMaximumBin(); - max = histo_py->GetXaxis()->GetBinCenter(binmax); - Double_t maxEntry = histo_py->GetBinContent(binmax); - auto* fgaus = new TF1("fgaus", "gaus(0)", max - 0.3, max + 0.3); /// TODO: find best range - fgaus->SetParameters(maxEntry, max, 20); - histo_py->Fit("fgaus", "QR0"); - Double_t offset1 = par->GetOffset1(); - Double_t offset2 = par->GetOffset2(); - max = fgaus->GetParameter(1) + offset1 - offset2; /// TODO: needs to be tested - // max = max+offset1-offset2; - veff = fTofdY / max; // effective speed of light in [cm/s] - LOG(info) << "Plane " << i + 1 << " Bar " << j + 1 << " offset " << par->GetOffset1(); - LOG(info) << "Plane " << i + 1 << " Bar " << j + 1 << " max " << max; - LOG(info) << "Plane " << i + 1 << " Bar " << j + 1 << " veff " << veff; - par->SetVeff(veff); - } - } - } - fHitPar->setChanged(); } - -void R3BTofDCal2HitPar::calcLambda(Double_t totLow, Double_t totHigh) -{ - TCanvas* cToTOffset = new TCanvas("cLambda", "cLambda", 10, 10, 1000, 900); - cToTOffset->Divide(1, 2); - for (Int_t i = 0; i < fNofPlanes; i++) - { - for (Int_t j = 0; j < fPaddlesPerPlane; j++) - { - Double_t offset = 0.; - auto par = fHitPar->GetModuleParAt(i + 1, j + 1); - if (fhSqrtQvsPosToT[i][j]) - { - LOG(info) << "Found histo SqrtQ_vs_PosToT_Plane_" << i + 1 << "_Bar_" << j + 1; - // auto* h = (TH2F*)histofilename->Get(Form("SqrtQ_vs_PosToT_Plane_%i_Bar_%i", i + 1, j + 1))->Clone(); - cToTOffset->cd(1); - fhSqrtQvsPosToT[i][j]->Draw("colz"); - auto* histo_py = (TH2F*)fhSqrtQvsPosToT[i][j]->ProjectionX("histo_py", totLow, totHigh, ""); - cToTOffset->cd(2); - histo_py->Draw(); - Int_t binmax = histo_py->GetMaximumBin(); - Double_t Max = histo_py->GetXaxis()->GetBinCenter(binmax); - TF1* fgaus = new TF1("fgaus", "gaus(0)", Max - 0.06, Max + 0.06); - histo_py->Fit("fgaus", "QR0"); - offset = fgaus->GetParameter(1); - fgaus->Draw("SAME"); - histo_py->SetAxisRange(Max - .5, Max + .5, "X"); - fhSqrtQvsPosToT[i][j]->SetAxisRange(Max - .5, Max + .5, "X"); - fhSqrtQvsPosToT[i][j]->SetAxisRange(totLow, totHigh, "Y"); - cToTOffset->Update(); - delete fgaus; - delete histo_py; - } - else - LOG(error) << "Missing histo plane " << i + 1 << " bar " << j + 1; - Double_t lambda = fTofdY / offset; - LOG(info) << " Plane " << i + 1 << " Bar " << j + 1 << " ToT Offset " << offset << " Lambda " << lambda - << "\n"; - par->SetLambda(lambda); - } - } - fHitPar->setChanged(); } - +fHitPar->setChanged(); +}*/ void R3BTofDCal2HitPar::smiley(TH2F* histo, Double_t min, Double_t max, Double_t* para) { - // This fits the smiley: Sqrt(q1*q2) returns position dependent charge, we fit that via pol3 and try to correct - Double_t y[1000], x[1000]; - Int_t n = 0; - for (Int_t j = 0; j <= 3; j++) - { - para[j] = 0; - } - TGraph* gr1 = new TGraph(); - TGraph* gr2 = new TGraph(); - TCanvas* cfit_smiley = new TCanvas("cfit_smiley", "fit smiley", 100, 100, 800, 800); - cfit_smiley->Clear(); - cfit_smiley->Divide(1, 4); - cfit_smiley->cd(1); - TH2F* histo1 = dynamic_cast(histo->Clone()); - histo1->Draw("colz"); - TH2F* histo2 = dynamic_cast(histo->Clone()); - histo2->RebinX(50); - histo2->GetYaxis()->SetRangeUser(fTofdTotLow, fTofdTotHigh); - // histo2->SetAxisRange(fTofdTotLow,fTofdTotHigh,"Y"); - cfit_smiley->cd(2); - histo2->Draw("colz"); - std::cout << "Searching for points to fit...\n"; - for (Int_t i = 1; i < histo2->GetNbinsX(); i++) - { - // std::cout<<"Bin "<GetNbinsX()<<" with cut: "<cd(2); - TLine* l = new TLine( - histo2->GetXaxis()->GetBinCenter(i), fTofdTotLow, histo2->GetXaxis()->GetBinCenter(i), fTofdTotHigh); - l->SetLineColor(kRed); - l->SetLineWidth(2.); - l->Draw(); - cfit_smiley->cd(3); - TH1F* histo_py = (TH1F*)histo2->ProjectionY("histo_py", i, i, ""); - histo_py->Draw(); - // cfit_smiley->Update(); - x[n] = histo2->GetXaxis()->GetBinCenter(i); - Int_t binmax = histo_py->GetMaximumBin(); - y[n] = histo_py->GetXaxis()->GetBinCenter(binmax); - - if ((x[n] < min || x[n] > max) || (y[n] < fTofdTotLow || y[n] > fTofdTotHigh)) - { - delete histo_py; - continue; - } - if (histo_py->GetMaximum() > 5) - { - n++; - delete l; - } - delete histo_py; - } - gr1 = new TGraph(n, x, y); - gr1->SetTitle("Points found for fitting; x position in cm; sqrt(tot1*tot2)"); - gr1->Draw("A*"); - std::cout << "Start fitting\n"; - TF1* f1 = new TF1("f1", "pol3", min, max); - f1->SetLineColor(2); - gr1->Fit("f1", "Q", "", min, max); - for (Int_t j = 0; j <= 3; j++) - { - para[j] = f1->GetParameter(j); - std::cout << "Parameter: " << para[j] << "\n"; - } - // fit again but with more information and better cuts - std::cout << "Fit again with more information\n"; - n = 0; - cfit_smiley->cd(4); - for (Int_t i = 1; i < histo2->GetNbinsX(); i++) - { - Double_t pos = histo2->GetXaxis()->GetBinCenter(i); - Double_t ymean = f1->Eval(pos); - histo2->SetAxisRange(ymean - 5., ymean + 5., "Y"); - histo2->Draw("colz"); - TH1F* histo_py = (TH1F*)histo2->ProjectionY("histo_py", i, i, ""); - histo_py->Draw(); - x[n] = histo2->GetXaxis()->GetBinCenter(i); - Int_t binmax = histo_py->GetMaximumBin(); - y[n] = histo_py->GetXaxis()->GetBinCenter(binmax); - if (histo_py->GetMaximum() > 2) - n++; - delete histo_py; - } - gr2 = new TGraph(n, x, y); - gr2->SetTitle("More information;x position in cm;sqrt(q1*q2)"); - gr2->Draw("A*"); - f1->DrawCopy("SAME"); - TF1* f2 = new TF1("f2", "pol3", min, max); - f2->SetParameters(para[0], para[1], para[2], para[3]); - f2->SetLineColor(3); - gr2->Fit("f2", "0Q", "", min, max); - f2->Draw("SAME"); - std::cout << "Will write:\n"; - for (Int_t j = 0; j <= 3; j++) - { - para[j] = f2->GetParameter(j); - std::cout << "Parameter: " << para[j] << "\n"; - } - histo2->GetYaxis()->SetRangeUser(fTofdTotLow, fTofdTotHigh); - auto legend = new TLegend(.9, 0.7, .99, 0.9); - legend->AddEntry("f1", "First Fit", "l"); - legend->AddEntry("f2", "Second Fit", "l"); - legend->Draw(); - cfit_smiley->Update(); - // gPad->WaitPrimitive(); - // gSystem->Sleep(3000); - delete histo1; - delete histo2; - delete gr1; - delete gr2; - delete f1; - delete f2; - delete cfit_smiley; + cout<<"Enter smiley function"<Clear(); + cfit_smiley->Divide(1, 4); + cfit_smiley->cd(1); + TH2F* histo1 = (TH2F*)histo->Clone(); + histo1->Draw("colz"); + TH2F* histo2 = (TH2F*)histo->Clone(); + histo2->RebinX(50); + histo2->GetYaxis()->SetRangeUser(fTofdTotLow, fTofdTotHigh); + // histo2->SetAxisRange(fTofdTotLow,fTofdTotHigh,"Y"); + cfit_smiley->cd(2); + histo2->Draw("colz"); + std::cout << "Searching for points to fit...\n"; + for (Int_t i = 1; i < histo2->GetNbinsX(); i++) + { + std::cout<<"Bin "<GetNbinsX()<<" with cut: "<cd(2); + TLine* l = new TLine( + histo2->GetXaxis()->GetBinCenter(i), fTofdTotLow, histo2->GetXaxis()->GetBinCenter(i), fTofdTotHigh); + l->SetLineColor(kRed); + l->SetLineWidth(2.); + l->Draw(); + cfit_smiley->cd(3); + TH1F* histo_py = (TH1F*)histo2->ProjectionY("histo_py", i, i, ""); + //histo_py->Rebin(rebin1); + histo_py->Draw(); + cfit_smiley->Update(); + x[n] = histo2->GetXaxis()->GetBinCenter(i); + Int_t binmax = histo_py->GetMaximumBin(); + Double_t y_check = histo_py->GetXaxis()->GetBinCenter(binmax); + + if ((x[n] < min || x[n] > max) || (y_check < fTofdTotLow || y_check > fTofdTotHigh)) + { + delete histo_py; + continue; + } + if(x[n] == 0.) cout<<"x=0 "; + if (histo_py->GetMaximum() > coarse_threshold) + { + Double_t fitrange = 5.; + TF1* fgaus = new TF1("fgausX", "gaus(0)", /*histo_py->GetXaxis()->GetBinCenter(binmax)*/y_check - fitrange, /*histo_py->GetXaxis()->GetBinCenter(binmax)*/y_check + fitrange); + fgaus->SetParameter(1, y_check); + fgaus->SetParLimits(1,y_check - (fitrange*1.25),y_check + (fitrange*1.25)); + histo_py->Fit(fgaus, "QR"); + Double_t fix = fgaus->GetParameter(1); + //y[n] = y_check; + y[n] = fix; + if(x[n] == 0.) cout<<"y "<Draw(); + cout<<"y "<SetTitle("Points found for fitting; x position in cm; sqrt(tot1*tot2)"); + gr1->Draw("A*"); + std::cout << "Start fitting\n"; + TF1* f1 = new TF1("f1", "pol3", min, max); + f1->SetLineColor(2); + gr1->Fit("f1", "Q", "", min, max); + for (Int_t j = 0; j <= 3; j++) + { + para[j] = f1->GetParameter(j); + std::cout << "Parameter: " << para[j] << "\n"; + } + // fit again but with more information and better cuts + std::cout << "Fit again with more information\n"; + n = 0; + cfit_smiley->cd(4); + for (Int_t i = 1; i < histo2->GetNbinsX(); i++) + { + Double_t pos = histo2->GetXaxis()->GetBinCenter(i); + Double_t ymean = f1->Eval(pos); + histo2->SetAxisRange(ymean - 5., ymean + 5., "Y"); + histo2->Draw("colz"); + TH1F* histo_py = (TH1F*)histo2->ProjectionY("histo_py", i, i, ""); + //histo_py->Rebin(rebin2); + histo_py->Draw(); + x[n] = histo2->GetXaxis()->GetBinCenter(i); + Int_t binmax = histo_py->GetMaximumBin(); + Double_t y_check = histo_py->GetXaxis()->GetBinCenter(binmax); + //y[n] = histo_py->GetXaxis()->GetBinCenter(binmax); + + if ((x[n] < min || x[n] > max) || (y_check < fTofdTotLow || y_check > fTofdTotHigh)) + { + delete histo_py; + continue; + } + if (histo_py->GetMaximum() > fine_threshold){ + Double_t fitrange = 5.; + TF1* fgaus2 = new TF1("fgaus2", "gaus(0)", y_check - (fitrange), y_check + (fitrange)); + fgaus2->SetParameter(1, y_check); + fgaus2->SetParLimits(1,y_check - (fitrange*1.25),y_check + (fitrange*1.25)); + histo_py->Fit(fgaus2, "QR"); + y[n] = fgaus2->GetParameter(1); + //y[n] = y_check; + cout<<"fine y "<Draw(); + n++; + } + delete histo_py; + } + for(int c = 0; c<=n;c++){ + cout<<"n_fine "<SetTitle("More information;x position in cm;sqrt(q1*q2)"); + gr2->Draw("A*"); + f1->DrawCopy("SAME"); + TF1* f2 = new TF1("f2", "pol3", min, max); + f2->SetParameters(para[0], para[1], para[2], para[3]); + f2->SetLineColor(3); + gr2->Fit("f2", "0Q", "", min, max); + f2->Draw("SAME"); + std::cout << "Will write:\n"; + for (Int_t j = 0; j <= 3; j++) + { + para[j] = f2->GetParameter(j); + std::cout << "Parameter"<GetYaxis()->SetRangeUser(fTofdTotLow, fTofdTotHigh); + auto legend = new TLegend(.9, 0.7, .99, 0.9); + legend->AddEntry("f1", "First Fit", "l"); + legend->AddEntry("f2", "Second Fit", "l"); + legend->Draw(); + cfit_smiley->Update(); + //gPad->WaitPrimitive(); + // gSystem->Sleep(3000); + char cCheck = 'y'; + while(!(cCheck == 'y' || cCheck == 'n')){ + cout<< "Continue? (y/n) "; + cin >> cCheck; + } + cout<<"1"<> rebin1; + cout<<"New coarse threshold: "; + cin>> coarse_threshold; + cout<<"New Rebin2: "; + cin>> rebin2; + cout<<"New fine threshold: "; + cin>> fine_threshold; + goto Redo; + } + cout<<"2"<Clear(); - cfit_exp->Divide(1, 3); - cfit_exp->cd(1); - TH2F* histo1 = dynamic_cast(histo->Clone()); - TH2F* histo2 = dynamic_cast(histo->Clone()); - histo1->Draw("colz"); - cfit_exp->cd(2); - for (Int_t i = 1; i < histo1->GetNbinsX() - 1; i++) - { - TH1F* histo_py = (TH1F*)histo1->ProjectionY("histo_py", i, i, ""); - histo_py->Draw(); - x[n] = histo1->GetXaxis()->GetBinCenter(i); - Int_t binmax = histo_py->GetMaximumBin(); - y[n] = histo_py->GetXaxis()->GetBinCenter(binmax); - if ((x[n] < -40. || x[n] > 40.) || y[n] < 50.) - { - delete histo_py; - continue; - } - if (histo_py->GetMaximum() > 5) - n++; - delete histo_py; - } - gr1 = new TGraph(n, x, y); - gr1->Draw("A*"); - TF1* f1 = new TF1("f1", "[0]*(exp(-[1]*(x+100.))+exp(-[2]*(x+100.)))+[3]", min, max); - f1->SetParameters(520., 0.001, 17234, -485.); - f1->SetLineColor(2); - gr1->Fit("f1", "", "", min, max); - for (Int_t j = 0; j <= 3; j++) - { - para[j] = f1->GetParameter(j); - std::cout << "Parameter: " << para[j] << "\n"; - } - // fit again but with more information and better cuts - n = 0; - cfit_exp->cd(3); - for (Int_t i = 1; i < histo2->GetNbinsX(); i++) - { - Double_t pos = histo2->GetXaxis()->GetBinCenter(i); - Double_t ymean = para[0] * (exp(-para[1] * (pos + 100.)) + exp(-para[2] * (pos + 100.))) + para[3]; - histo2->SetAxisRange(ymean - 5., ymean + 5., "Y"); - histo2->Draw("colz"); - TH1F* histo_py = (TH1F*)histo2->ProjectionY("histo_py", i, i, ""); - histo_py->Draw(); - x[n] = histo2->GetXaxis()->GetBinCenter(i); - Int_t binmax = histo_py->GetMaximumBin(); - y[n] = histo_py->GetXaxis()->GetBinCenter(binmax); - if (histo_py->GetMaximum() > 2) - n++; - delete histo_py; - } - gr2 = new TGraph(n, x, y); - gr2->Draw("A*"); - TF1* f2 = new TF1("f2", "[0]*(exp(-[1]*(x+100.))+exp(-[2]*(x+100.)))+[3]", min, max); - f2->SetParameters(para[0], para[1], para[2], para[3]); - f2->SetLineColor(2); - gr2->Fit("f2", "", "", min, max); - for (Int_t j = 0; j <= 3; j++) - { - para[j] = f2->GetParameter(j); - std::cout << "Parameter: " << para[j] << "\n"; - } - cfit_exp->Update(); - // gPad->WaitPrimitive(); - // gSystem->Sleep(3000); - delete gr1; - delete gr2; - delete f1; - delete f2; + // This fits the exponential decay of the light in a paddle. The 2 PMTs are fit with the same function but one + // side will deliver negative attenuation parameters and the other positive. + Double_t y[1000], x[1000]; + Int_t n = 0; + for (Int_t j = 0; j <= 3; j++) + { + para[j] = 0; + } + TGraph* gr1 = new TGraph(); + TGraph* gr2 = new TGraph(); + TCanvas* cfit_exp = new TCanvas("cfit_exp", "fit exponential", 100, 100, 800, 800); + cfit_exp->Clear(); + cfit_exp->Divide(1, 3); + cfit_exp->cd(1); + TH2F* histo1 = (TH2F*)histo->Clone(); + TH2F* histo2 = (TH2F*)histo->Clone(); + histo1->Draw("colz"); + cfit_exp->cd(2); + for (Int_t i = 1; i < histo1->GetNbinsX() - 1; i++) + { + TH1F* histo_py = (TH1F*)histo1->ProjectionY("histo_py", i, i, ""); + histo_py->Draw(); + x[n] = histo1->GetXaxis()->GetBinCenter(i); + Int_t binmax = histo_py->GetMaximumBin(); + y[n] = histo_py->GetXaxis()->GetBinCenter(binmax); + if ((x[n] < -40. || x[n] > 40.) || y[n] < 50.) + { + delete histo_py; + continue; + } + if (histo_py->GetMaximum() > 5) + n++; + delete histo_py; + } + gr1 = new TGraph(n, x, y); + gr1->Draw("A*"); + TF1* f1 = new TF1("f1", "[0]*(exp(-[1]*(x+100.))+exp(-[2]*(x+100.)))+[3]", min, max); + f1->SetParameters(520., 0.001, 17234, -485.); + f1->SetLineColor(2); + gr1->Fit("f1", "", "", min, max); + for (Int_t j = 0; j <= 3; j++) + { + para[j] = f1->GetParameter(j); + std::cout << "Parameter: " << para[j] << "\n"; + } + // fit again but with more information and better cuts + n = 0; + cfit_exp->cd(3); + for (Int_t i = 1; i < histo2->GetNbinsX(); i++) + { + Double_t pos = histo2->GetXaxis()->GetBinCenter(i); + Double_t ymean = para[0] * (exp(-para[1] * (pos + 100.)) + exp(-para[2] * (pos + 100.))) + para[3]; + histo2->SetAxisRange(ymean - 5., ymean + 5., "Y"); + histo2->Draw("colz"); + TH1F* histo_py = (TH1F*)histo2->ProjectionY("histo_py", i, i, ""); + histo_py->Draw(); + x[n] = histo2->GetXaxis()->GetBinCenter(i); + Int_t binmax = histo_py->GetMaximumBin(); + y[n] = histo_py->GetXaxis()->GetBinCenter(binmax); + if (histo_py->GetMaximum() > 2) + n++; + delete histo_py; + } + gr2 = new TGraph(n, x, y); + gr2->Draw("A*"); + TF1* f2 = new TF1("f2", "[0]*(exp(-[1]*(x+100.))+exp(-[2]*(x+100.)))+[3]", min, max); + f2->SetParameters(para[0], para[1], para[2], para[3]); + f2->SetLineColor(2); + gr2->Fit("f2", "", "", min, max); + for (Int_t j = 0; j <= 3; j++) + { + para[j] = f2->GetParameter(j); + std::cout << "Parameter: " << para[j] << "\n"; + } + cfit_exp->Update(); + // gPad->WaitPrimitive(); + // gSystem->Sleep(3000); + delete gr1; + delete gr2; + delete f1; + delete f2; } - -void R3BTofDCal2HitPar::FinishTask() -{ - if (fParameter == 1) - { - // Determine time offset of the 2 PMTs of one paddle. This procedure - // assumes a sweep run in the middle of the ToF wall horizontally. - // Since all paddles are mounted vertically one can determine the offset. - // Half of the offset is added to PM1 and half to PM2. - LOG(info) << "Calling function calcOffset"; - calcOffset(); - // Determine ToT offset between top and bottom PMT - LOG(info) << "Calling function calcToTOffset"; - calcToTOffset(fTofdTotLow, fTofdTotHigh); - // Determine sync offset between paddles - LOG(info) << "Calling function calcSync"; - calcSync(); - LOG(error) << "Call walk correction before next step!"; - - for (Int_t i = 0; i < fNofPlanes; i++) - { - for (Int_t j = 0; j < fPaddlesPerPlane; j++) - { - auto par = fHitPar->GetModuleParAt(i + 1, j + 1); - Int_t binmax = fh1_tofsync[i][j]->GetMaximumBin(); - auto tofsync = fh1_tofsync[i][j]->GetXaxis()->GetBinCenter(binmax); - - TF1* fgauss = new TF1("fgaus", "gaus(0)", tofsync - 0.25, tofsync + 0.25); - fh1_tofsync[i][j]->Fit("fgaus", "QR"); - auto tof_offset = fgauss->GetParameter(1); - - par->SetTofSyncOffset(tof_offset - fMeanTof); - LOG(info) << " Plane " << i + 1 << " Bar " << j + 1 << " Tof-Sync " << tof_offset; - } - } - } - else if (fParameter == 2) - { - // Determine effective speed of light in [cm/s] for each paddle - LOG(info) << "Calling function"; - calcVeff(); - // Determine light attenuation lambda for each paddle - LOG(info) << "Calling function calcLambda"; - calcLambda(fTofdTotLow, fTofdTotHigh); - } - else if (fParameter == 3) +/* void R3BTofDCal2HitPar::calcZoffset() { - // calculation of position dependend charge - if (fTofdSmiley) - { - LOG(info) << "Calling function smiley"; - Double_t para2[4]; - for (Int_t i = 0; i < 4; i++) - para2[i] = 0.; - Double_t min2 = -40.; // -40 effective bar length - Double_t max2 = 40.; // 40 effective bar length = 80 cm - // we will use 50 here for some fit safety margin - for (Int_t i = 0; i < fNofPlanes; i++) - { - for (Int_t j = 0; j < fPaddlesPerPlane; j++) - { - if (fhSqrtQvsPosToT[i][j]) - { - LOG(info) << "Calling Plane " << i + 1 << " Bar " << j + 1; - auto par = fHitPar->GetModuleParAt(i + 1, j + 1); - smiley(fhSqrtQvsPosToT[i][j], min2, max2, para2); - par->SetPola(para2[0]); - par->SetPolb(para2[1]); - par->SetPolc(para2[2]); - par->SetPold(para2[3]); - } - } - } - fHitPar->setChanged(); - } - else - { - LOG(info) << "Calling function doubleExp"; - Double_t para[4]; - for (Int_t i = 0; i < 4; i++) - para[i] = 0.; - Double_t min = -40.; // effective bar length - Double_t max = 40.; // effective bar length = 80 cm - - for (Int_t i = 0; i < fNofPlanes; i++) - { - for (Int_t j = 0; j < fPaddlesPerPlane; j++) - { - if (fhTot1vsPos[i][j]) - { - auto par = fHitPar->GetModuleParAt(i + 1, j + 1); - doubleExp(fhTot1vsPos[i][j], min, max, para); - Double_t offset1 = par->GetOffset1(); - Double_t offset2 = par->GetOffset2(); - Double_t veff = par->GetVeff(); - Double_t sync = par->GetSync(); - par->SetPar1a(para[0]); - par->SetPar1b(para[1]); - par->SetPar1c(para[2]); - par->SetPar1d(para[3]); - } - if (fhTot2vsPos[i][j]) - { - auto par = fHitPar->GetModuleParAt(i + 1, j + 1); - doubleExp(fhTot2vsPos[i][j], min, max, para); - Double_t offset1 = par->GetOffset1(); - Double_t offset2 = par->GetOffset2(); - Double_t veff = par->GetVeff(); - Double_t sync = par->GetSync(); - par->SetPar2a(para[0]); - par->SetPar2b(para[1]); - par->SetPar2c(para[2]); - par->SetPar2d(para[3]); - } - } - } - fHitPar->setChanged(); - } + TCanvas* cZoff = new TCanvas("cZoff", "cZoff", 10, 10, 1000, 900); + cZoff->Divide(2, 2); + for (Int_t i = 0; i < fNofPlanes; i++) + { + for (Int_t j = 0; j < fPaddlesPerPlane; j++) + { + //Double_t max = 0.; + //Double_t veff = 7.; + if (fZoffset[i][j]) + { + LOG(info) << "Found histo Time_Diff_Plane_" << i + 1; + auto par = fHitPar->GetModuleParAt(i + 1, j + 1); + if (!par) + { + LOG(warning) << "Hit par not found, Plane: " << i + 1 << ", Bar: " << j + 1; + continue; + } + // auto* h = (TH2F*)histofilename->Get(Form("Time_Diff_Plane_%i", i + 1))->Clone(); + cZoff->cd(i + 1); + // h->Draw("colz"); + TH1F* histo_zy = (TH1F*)fhZoffset[i][j]->ProjectionY("histo_zy", j + 2, j + 2, ""); + Int_t binmax = histo_zy->GetMaximumBin(); + max = histo_zy->GetXaxis()->GetBinCenter(binmax); + Double_t maxEntry = histo_zy->GetBinContent(binmax); + auto* fgaus = new TF1("fgaus", "gaus(0)", max - 0.4, max + 0.4); /// TODO: find best range + fgaus->SetParameters(maxEntry, max, 20); + histo_zy->Fit("fgaus", "QR0"); + //Double_t offset1 = par->GetOffset1(); + //Double_t offset2 = par->GetOffset2(); + Double_t cent = fgaus->GetParameter(1) /// TODO: needs to be tested + // max = max+offset1-offset2; + Zoffset = fMaxQ - cent; // effective speed of light in [cm/s] + //LOG(info) << "Plane " << i + 1 << " Bar " << j + 1 << " offset " << par->GetOffset1(); + //LOG(info) << "Plane " << i + 1 << " Bar " << j + 1 << " max " << max; + //LOG(info) << "Plane " << i + 1 << " Bar " << j + 1 << " veff " << veff; + par->SetZoffset(Zoffset); + } + } + } + fHitPar->setChanged(); } +*/ - if (fParameter == 4) - { - // Z correction for each plane - LOG(warn) << "Calling function zcorr"; - Double_t para[8]; - Double_t pars[3]; - pars[0] = 0.; - pars[1] = 0.; - pars[2] = 0.; - Int_t min = 0, max = fMaxQ; // select range for peak search - for (Int_t i = 0; i < fNofPlanes; i++) - { - for (Int_t j = 0; j < fPaddlesPerPlane; j++) - { - if (fhQvsPos[i][j]) - { - auto par = fHitPar->GetModuleParAt(i + 1, j + 1); - std::cout << "Calling Plane: " << i + 1 << " Bar " << j + 1 << "\n"; - Int_t index = i * fPaddlesPerPlane + j; - zcorr(fhQvsPos[i][j], min, max, pars, index); - std::cout << "Write parameter: " << pars[0] << " " << pars[1] << " " << pars[2] << "\n"; - par->SetPar1za(pars[0]); - par->SetPar1zb(pars[1]); - par->SetPar1zc(pars[2]); - } - } - } - fHitPar->setChanged(); - } - for (Int_t i = 0; i < N_TOFD_HIT_PLANE_MAX; i++) - { - if (fh_tofd_TotPm[i]) - fh_tofd_TotPm[i]->Write(); // control histogram for ToT - if (fhTsync[i]) - fhTsync[i]->Write(); // histogram for sync calculation - if (fhTdiff[i]) - fhTdiff[i]->Write(); // histogram for offset and veff calculation - for (Int_t j = 0; j < N_TOFD_HIT_PADDLE_MAX; j++) - { - if (fh1_tofsync[i][j]) - fh1_tofsync[i][j]->Write(); // histogram for ToF sync calculation - if (fhLogTot1vsLogTot2[i][j]) - fhLogTot1vsLogTot2[i][j]->Write(); // control histogram Log(ToT) Pm1 vs Log(ToT) Pm2 - if (fhSqrtQvsPosToT[i][j]) - fhSqrtQvsPosToT[i][j]->Write(); // histogram for ToT offset calculation - if (fhQvsPos[i][j]) - fhQvsPos[i][j]->Write(); // histogram for charge fit - if (!fTofdSmiley) - { - if (fhTot1vsPos[i][j]) - fhTot1vsPos[i][j]->Write(); // histogram for position dependence of charge 1 - if (fhTot2vsPos[i][j]) - fhTot2vsPos[i][j]->Write(); // histogram for position dependence of charge 2 - } - /* - if (fhTot1vsTot2[i][j]) - fhTot1vsTot2[i][j]->Write(); // control histogram ToT Pm1 vs ToT Pm2 - */ - } - } +void R3BTofDCal2HitPar::FinishTask() +{ + + if (fParameter == 10) + { + // Determine time offset of the 2 PMTs of one paddle. This procedure + // assumes a sweep run in the middle of the ToF wall horizontally. + // Since all paddles are mounted vertically one can determine the offset. + // Half of the offset is added to PM1 and half to PM2. + // Determine ToT offset between top and bottom PMT + LOG(info) << "Calling function calcToTOffset"; + calcToTOffset(fTofdTotLow, fTofdTotHigh); + // Determine sync offset between paddles + + for (Int_t i = 0; i < fNofPlanes; i++) + { + fhTofsync[i]->Write(); + fh_tofd_TotPm[i][0]->Write(); + fh_tofd_TotPm[i][1]->Write(); + fh_tofd_TotPm[i][2]->Write(); + } + } + if (fParameter == 1) + { + // Determine time offset of the 2 PMTs of one paddle. This procedure + // assumes a sweep run in the middle of the ToF wall horizontally. + // Since all paddles are mounted vertically one can determine the offset. + // Half of the offset is added to PM1 and half to PM2. + LOG(info) << "Calling function calcOffset"; + calcOffset(); + // Determine ToT offset between top and bottom PMT + LOG(info) << "Calling function calcToTOffset"; + calcToTOffset(fTofdTotLow, fTofdTotHigh); + // Determine sync offset between paddles + LOG(info) << "Calling function calcSync"; + calcSync(); + LOG(error) << "Call walk correction before next step!"; + + /* for (Int_t i = 0; i < fNofPlanes; i++) + { + for (Int_t j = 0; j < fPaddlesPerPlane; j++) + { + auto par = fHitPar->GetModuleParAt(i + 1, j + 1); + Int_t binmax = fh1_tofsync[i][j]->GetMaximumBin(); + auto tofsync = fh1_tofsync[i][j]->GetXaxis()->GetBinCenter(binmax); + par->SetTofSyncOffset(tofsync - fMeanTof); + LOG(info) << " Plane " << i + 1 << " Bar " << j + 1 << " Tof-Sync " << tofsync; + fh1_tofsync[i][j]->Write(); + } + fhTofsync[i]->Write(); + fh_tofd_TotPm[i][0]->Write(); + fh_tofd_TotPm[i][1]->Write(); + fh_tofd_TotPm[i][2]->Write(); + }*/ + } + else if(fParameter == 0) + { + // Calculate walk corrections + LOG(info) << "Calling function"; + //calcwalk(); + + + } + else if (fParameter == 2) + { + //Calculate ToF + LOG(info) << "Calculate ToF"; + for (Int_t i = 0; i < fNofPlanes; i++) + { + for (Int_t j = 0; j < fPaddlesPerPlane; j++) + { + auto par = fHitPar->GetModuleParAt(i + 1, j + 1); + Int_t binmax = fh1_tofsync[i][j]->GetMaximumBin(); + auto tofsync = fh1_tofsync[i][j]->GetXaxis()->GetBinCenter(binmax); + par->SetTofSyncOffset(tofsync - fMeanTof); + LOG(info) << " Plane " << i + 1 << " Bar " << j + 1 << " Tof-Sync " << tofsync; + fh1_tofsync[i][j]->Write(); + } + fhTofsync[i]->Write(); + fh_tofd_TotPm[i][0]->Write(); + fh_tofd_TotPm[i][1]->Write(); + fh_tofd_TotPm[i][2]->Write(); + } + // Determine effective speed of light in [cm/s] for each paddle + LOG(info) << "Calling function"; + calcVeff(); + // Determine light attenuation lambda for each paddle + LOG(info) << "Calling function calcLambda"; + calcLambda(fTofdTotLow, fTofdTotHigh); + + for(int i=0; i<4; i++){ + fhTofsync[i]->Write(); + fhTdiff[i]->Write(); + fh_tofd_TotPm[i][1]->Write(); + fh_tofd_TotPm[i][0]->Write(); + fh_tofd_TotPm[i][2]->Write(); + for(int j=0;j<44;j++){ + fhSqrtQvsPosToT[i][j]->Write(); + } + } + } + else if (fParameter == 3) + { + // calculation of position dependend charge + if (fTofdSmiley) + { + + // LOG(info) << "Calling function smiley"; + // Double_t para2[4]; + // for (Int_t i = 0; i < 4; i++) + // para2[i] = 0.; + // Double_t min2 = -40.; // -40 effective bar length + // Double_t max2 = 40.; // 40 effective bar length = 80 cm + // // we will use 50 here for some fit safety margin + // for (Int_t i = 0; i < /*fNofPlanes*/1; i++) + // { + // for (Int_t j = 0; j < fPaddlesPerPlane; j++) + // { + // if (fhSqrtQvsPosToT[i][j]) + // { + // LOG(info) << "Calling Plane " << i + 1 << " Bar " << j + 1; + // auto par = fHitPar->GetModuleParAt(i + 1, j + 1); + // smiley(fhSqrtQvsPosToT[i][j], min2, max2, para2); +// + // cout<<"Plane "<GetToTOffset1(); + // Double_t offset2 = par->GetToTOffset2(); + // Double_t veff = par->GetVeff(); + // Double_t sync = par->GetSync(); + // par->SetPola(para2[0]); + // par->SetPolb(para2[1]); + // par->SetPolc(para2[2]); + // par->SetPold(para2[3]); + // + // fhSqrtQvsPosToT[i][j]->Write(); + // } + // } + // } + + for(int i=0;iWrite(); + fhTot2vsPos[i][j]->Write(); + fhTot1vsPos[i][j]->Write(); + } + } + //fHitPar->setChanged(); + } + else + { + //LOG(info) << "Calling function doubleExp"; + //Double_t para[4]; + //for (Int_t i = 0; i < 4; i++) + // para[i] = 0.; + //Double_t min = -40.; // effective bar length + //Double_t max = 40.; // effective bar length = 80 cm +// + //for (Int_t i = 0; i < fNofPlanes; i++) + //{ + // for (Int_t j = 0; j < fPaddlesPerPlane; j++) + // { + // if (fhTot1vsPos[i][j]) + // { + // auto par = fHitPar->GetModuleParAt(i + 1, j + 1); + // doubleExp(fhTot1vsPos[i][j], min, max, para); + // Double_t offset1 = par->GetOffset1(); + // Double_t offset2 = par->GetOffset2(); + // Double_t veff = par->GetVeff(); + // Double_t sync = par->GetSync(); + // par->SetPar1a(para[0]); + // par->SetPar1b(para[1]); + // par->SetPar1c(para[2]); + // par->SetPar1d(para[3]); + // } + // if (fhTot2vsPos[i][j]) + // { + // auto par = fHitPar->GetModuleParAt(i + 1, j + 1); + // doubleExp(fhTot2vsPos[i][j], min, max, para); + // Double_t offset1 = par->GetOffset1(); + // Double_t offset2 = par->GetOffset2(); + // Double_t veff = par->GetVeff(); + // Double_t sync = par->GetSync(); + // par->SetPar2a(para[0]); + // par->SetPar2b(para[1]); + // par->SetPar2c(para[2]); + // par->SetPar2d(para[3]); + // } + // } + //} + for(int i=0;iWrite(); + fhTot1vsPos[i][j]->Write(); + } + } + //fHitPar->setChanged(); + } + } + + if(fParameter == -3){ + if (fTofdSmiley) + { + + LOG(info) << "Calling function smiley"; + TFile* histofilename = TFile::Open(fHistoFile); + Double_t para2[4]; + for (Int_t i = 0; i < 4; i++) + para2[i] = 0.; + Double_t min2 = -40.; // -40 effective bar length + Double_t max2 = 40.; // 40 effective bar length = 80 cm + // we will use 50 here for some fit safety margin + for (Int_t i = 0; i < fNofPlanes; i++) + { + for (Int_t j = 0; j < fPaddlesPerPlane; j++) + { + //if (fhSqrtQvsPosToT[i][j]) + if (histofilename->Get(Form("SqrtQ_vs_PosToT_Plane_%i_Bar_%i", i + 1, j + 1))) + { + + for (Int_t l = 0; l < 4; l++) + para2[l] = 0.; + LOG(info) << "Calling Plane " << i + 1 << " Bar " << j + 1; + auto par = fHitPar->GetModuleParAt(i + 1, j + 1); + smiley(/*fhSqrtQvsPosToT[i][j]*/(TH2F*)histofilename->Get(Form("SqrtQ_vs_PosToT_Plane_%i_Bar_%i", i + 1, j + 1)), min2, max2, para2); + + cout<<"Plane "<GetToTOffset1(); + Double_t offset2 = par->GetToTOffset2(); + Double_t veff = par->GetVeff(); + Double_t sync = par->GetSync(); + cout<<"3.5"<SetPola(para2[0]); + cout<<"3.6"<SetPolb(para2[1]); + cout<<"3.7"<SetPolc(para2[2]); + cout<<"3.8"<SetPold(para2[3]); + cout<<"4"<Write(); + } + } + } + fHitPar->setChanged(); + + } + else + { + LOG(info) << "Calling function doubleExp"; + Double_t para[4]; + for (Int_t i = 0; i < 4; i++) + para[i] = 0.; + Double_t min = -40.; // effective bar length + Double_t max = 40.; // effective bar length = 80 cm + + for (Int_t i = 0; i < fNofPlanes; i++) + { + for (Int_t j = 0; j < fPaddlesPerPlane; j++) + { + if (fhTot1vsPos[i][j]) + { + auto par = fHitPar->GetModuleParAt(i + 1, j + 1); + doubleExp(fhTot1vsPos[i][j], min, max, para); + Double_t offset1 = par->GetOffset1(); + Double_t offset2 = par->GetOffset2(); + Double_t veff = par->GetVeff(); + Double_t sync = par->GetSync(); + par->SetPar1a(para[0]); + par->SetPar1b(para[1]); + par->SetPar1c(para[2]); + par->SetPar1d(para[3]); + fhTot1vsPos[i][j]->Write(); + } + if (fhTot2vsPos[i][j]) + { + auto par = fHitPar->GetModuleParAt(i + 1, j + 1); + doubleExp(fhTot2vsPos[i][j], min, max, para); + Double_t offset1 = par->GetOffset1(); + Double_t offset2 = par->GetOffset2(); + Double_t veff = par->GetVeff(); + Double_t sync = par->GetSync(); + par->SetPar2a(para[0]); + par->SetPar2b(para[1]); + par->SetPar2c(para[2]); + par->SetPar2d(para[3]); + fhTot2vsPos[i][j]->Write(); + } + } + } + fHitPar->setChanged(); + } + } + + if (fParameter == 4) + { + //Write Histograms for Charge correction. + //Call fParameter = 5 to read these Histogramms. + for (Int_t i = 0; i < fNofPlanes; i++) + { + + //if(i > 0) continue; + + + for (Int_t j = 0; j < fPaddlesPerPlane; j++) + { + if (fhQvsPos[i][j]) + { + /*auto par = fHitPar->GetModuleParAt(i + 1, j + 1); + std::cout << "Calling Plane: " << i + 1 << " Bar " << j + 1 << "\n"; + Int_t index = i * fPaddlesPerPlane + j; + int a = zcorr(fhQvsPos[i][j], min, max, pars, index); + std::cout << "Write parameter: " << pars[0] << " " << pars[1] << " " << pars[2] << "\n"; + + if(a > 0) j = j - a; + + Double_t offset1 = par->GetOffset1(); + Double_t offset2 = par->GetOffset2(); + Double_t veff = par->GetVeff(); + Double_t sync = par->GetSync(); + par->SetPar1za(pars[0]); + par->SetPar1zb(pars[1]); + par->SetPar1zc(pars[2]); + */ + fhQvsPos[i][j]->Write(); + } + if(fhQXvsPos[i][j]){ + fhQXvsPos[i][j]->Write(); + } + } + } + fHitPar->setChanged(); + } + + if(fParameter == 5){ + // Z correction for each plane + LOG(warning) << "Calling function zcorr"; + TFile* histofilename = TFile::Open(fHistoFile); + Double_t para[8]; + Double_t pars[3]; + pars[0] = 0.; + pars[1] = 0.; + pars[2] = 0.; + //Int_t min = fMinQ, max = fMaxQ; // select range for peak search + Double_t min = 0., max = 16.; // select range for peak search + cout<<"1"< 0) continue; + cout<<"2"<Get(Form("QX_vs_Pos_Plane_%i_Bar_%i", i + 1, j + 1))) + { + cout<<"4"<GetModuleParAt(i + 1, j + 1); + std::cout << "Calling Plane: " << i + 1 << " Bar " << j + 1 << "\n"; + Int_t index = i * fPaddlesPerPlane + j; + //int a = zcorr(fhQvsPos[i][j], min, max, pars, index); + int a = zcorr((TH2F*)histofilename->Get(Form("QX_vs_Pos_Plane_%i_Bar_%i", i + 1, j + 1)), min, max, pars, i, j, index); + std::cout << "Write parameter: " << pars[0] << " " << pars[1] << " " << pars[2] << "\n"; + cout<<"5"< 0) j = j - a; + + Double_t offset1 = par->GetOffset1(); + Double_t offset2 = par->GetOffset2(); + Double_t veff = par->GetVeff(); + Double_t sync = par->GetSync(); + par->SetPar1za(pars[0]); + par->SetPar1zb(pars[1]); + par->SetPar1zc(pars[2]); + //fhQvsPos[i][j]->Write(); + cout<<"6"< 0) continue; + cout<<"2"<Get(Form("QX_vs_Pos_Plane_%i_Bar_%i", i + 1, j + 1))) + { + cout<<"4"<GetModuleParAt(i + 1, j + 1); + std::cout << "Calling Plane: " << i + 1 << " Bar " << j + 1 << "\n"; + Int_t index = i * fPaddlesPerPlane + j; + int a = zcorr(fhQXvsPos[i][j], min, max, pars, i, j, index); + //int a = zcorr((TH2F*)histofilename->Get(Form("QX_vs_Pos_Plane_%i_Bar_%i", i + 1, j + 1)), min, max, pars, i, j, index); + cout<<"5"< 0) j = j - a; + + Double_t offset1 = par->GetOffset1(); + Double_t offset2 = par->GetOffset2(); + Double_t veff = par->GetVeff(); + Double_t sync = par->GetSync(); + para[0] = par->GetPar1za(); + para[1] = par->GetPar1zb(); + + std::cout << "Old parameter: " << para[0] << " " << para[1] << " " << pars[2] << "\n"; + std::cout << "New parameter: " << pars[0] << " " << pars[1] << " " << pars[2] << "\n"; + std::cout << "Write parameter: " << pars[0]*para[1]+para[0] << " " << pars[1]*para[1] << " " << pars[2] << "\n"; + + par->SetPar1za(pars[0]*para[1]+para[0]); + par->SetPar1zb(pars[1]*para[1]); + par->SetPar1zc(pars[2]); + //fhQvsPos[i][j]->Write(); + cout<<"6"< 0) continue; + + + + //if (fhQvsPos[i][j]) + if (histofilename->Get(Form("tofd_hit_Q_plane_%d", i + 1))) + { + pars[0] = 0.; + pars[1] = 1.; + pars[2] = 0.; + para[0] = 0.; + para[1] = 0.; + para[2] = 0.; + TH2F* histo = (TH2F*)histofilename->Get(Form("tofd_hit_Q_plane_%d", i + 1)); + for (Int_t j = 0; j < fPaddlesPerPlane; j++) + { + if(j>41) continue; + TCanvas* c1 = new TCanvas("c1", "c1", 150, 150, 1200, 1200); + c1->Clear(); + c1->Divide(1, 3); + c1->cd(1); + gPad->SetLogz(); + TH2F* histo1 = (TH2F*)histo->Clone(); + histo1->SetAxisRange(7.69, 8.15, "Y"); + TH2F* histo2 = (TH2F*)histo->Clone(); + histo2->SetAxisRange(1.8, 2.1, "Y"); + TH1F* histo_py = (TH1F*)histo->ProjectionY("histo_py", j+1, j+1, ""); + TH1F* histo_py8 = (TH1F*)histo1->ProjectionY("histo_py8", j+1, j+1, ""); + //histo_py8->GetXaxis()->SetRange(7.69, 8.15); + TH1F* histo_py2 = (TH1F*)histo2->ProjectionY("histo_py2", j+1, j+1, ""); + //histo_py2->GetXaxis()->SetRange(1.8, 2.1); + auto par = fHitPar->GetModuleParAt(i + 1, j + 1); + std::cout << "Calling Plane: " << i + 1 << " Bar " << j + 1 << "\n"; + + histo->Draw("colz"); + c1->cd(2); + histo_py8->Draw(); + c1->cd(3); + histo_py2->Draw(); + //histo_py8->SetAxisRange(7.69, 8.15, "X"); + //dhisto_py2->SetAxisRange(1.8, 2.1,"X"); + + //Int_t binmax8 = histo_py8->GetMaximum(); + //Double_t Max8 = histo_py8->GetXaxis()->GetBinCenter(binmax8); + //Int_t binmax2 = histo_py2->GetMaximum(); + //Double_t Max2 = histo_py2->GetXaxis()->GetBinCenter(binmax2); + + + Int_t binmax8 = histo_py8->GetMaximumBin(); + + //histo_py8->GetXaxis()->SetRange(7.69, 8.15); + //Int_t binmax8 = histo->GetYaxis()->GetMaximum(); + //histo_py->SetAxisRange(1.8, 2.1, "X"); + Int_t binmax2 = histo_py2->GetMaximumBin(); + + + + //histo_py->SetAxisRange(0., 12., "X"); + Double_t Max8 = histo_py8->GetBinCenter(binmax8); + Double_t Max2 = histo_py2->GetBinCenter(binmax2); + + TF1* fgaus1 = new TF1("fgaus1", "gaus(0)", Max8 - 2., Max8 + 2.); + histo_py8->Fit("fgaus1", "QR0"); + Double_t m8 = fgaus1->GetParameter(1); + TF1* fgaus2 = new TF1("fgaus2", "gaus(0)", Max2 - 2., Max2 + 2.); + histo_py2->Fit("fgaus2", "QR0"); + Double_t m2 = fgaus2->GetParameter(1); + + + cout<<"binmax8 "<GetPar1za(); + para[1] = par->GetPar1zb(); + + pars[0] = para[0] * f - d; + pars[1] = para[1] * f; + + cout<<"Old pars: "<SetPar1za(pars[0]); + par->SetPar1zb(pars[1]); + + char check = 'x'; + + if(j==-1){ + cout<GetBinContent(0); + gPad->WaitPrimitive(); + gSystem->Sleep(3000); + } + cout<<"Checking"; + //cin>>check; + delete histo1; + delete histo2; + delete histo_py8; + delete histo_py2; + } + + } + + } + fHitPar->setChanged();*/ + } + if(fParameter==7){ + TFile* histofilename = TFile::Open(fHistoFile); + Double_t para[8]; + Double_t pars[3]; + pars[0] = 0.; + pars[1] = 0.; + pars[2] = 0.; + //Int_t min = fMinQ, max = fMaxQ; // select range for peak search + Double_t min = 0., max = 12.; // select range for peak search + cout<<"1"< 0) continue; + + + for (Int_t j = 0; j < fPaddlesPerPlane; j++) + { + + //if (fhQvsPos[i][j]) + if (histofilename->Get(Form("tofd_hit_Q_plane_%d", i + 1))) + { + + auto parY = fHitPar->GetModuleParAt(i + 1, j + 1); + std::cout << "Calling Plane: " << i + 1 << " Bar " << j + 1 << "\n"; + Int_t index = i * fPaddlesPerPlane + j; + //int a = zcorr(fhQvsPos[i][j], min, max, pars, index); + auto* hist = (TH2F*)histofilename->Get(Form("tofd_hit_Q_plane_%d", i + 1)); + auto* histo = hist->ProjectionY("p_y", j+2, j+2,""); + + //begin zcorr + if (histo->GetEntries() < 10) + { + R3BLOG(warning, "Nb of events below 10 with "<< histo->GetEntries() << " entries for histo with index" << index); + + continue; + } + + Double_t par[3000] = { 0 }; + Double_t x[3000] = { 0 }; + char strName[255]; + int a = 0; + sprintf(strName, "canvas_%d", index); + TCanvas* c1 = new TCanvas(strName, "", 100, 100, 800, 800); + c1->Divide(1, 3); + c1->cd(1); + gPad->SetLogz(); + hist->Draw("colz"); + // Projection of charge axis + auto* h1 = (TH2F*)histo->Clone(); + + c1->cd(2); + gPad->SetLogy(); + h1->Draw(); + // Use TSpectrum to find the peak candidates + Double_t threshhold = 0.001; + Int_t sigma = 10; + Double_t sqrtf = 1.; + Int_t nPeaks = fNbZPeaks; + doagainfit: + TSpectrum* s = new TSpectrum(nPeaks); + Int_t nfound = s->Search(h1, sigma, "nobackground", threshhold); + std::cout << "Found " << nfound << " candidate peaks to fit\n"; + + if (nfound == 0) + { + delete s; + delete c1; + delete histo; + delete h1; + continue; + } + + c1->Update(); + // Eliminate background peaks + Int_t mPeaks = 0; + Double_t* xpeaks = s->GetPositionX(); + for (Int_t p = 0; p <= nfound; p++) + { + Float_t xp = xpeaks[p]; + Int_t bin = h1->GetXaxis()->FindBin(xp); + Float_t yp = h1->GetBinContent(bin); + if (yp - sqrtf * TMath::Sqrt(yp) < 1.){ + cout<<"jump "<Update(); + + if (mPeaks < 2 && 0) + { + pars[0] = 0.; + //pars[1] = fMaxQ / peaks[0]; + pars[1] = 8. / peaks[0]; + pars[2] = 0.; + + delete s; + delete c1; + delete histo; + delete h1; + continue; + } + cout<<"mPeaks "<> peakscheck; + if(peakscheck == "n"){ + delete s; + std::cout << "New sigma: "; + std::cin >> sigma; + std::cout << "New threshhold: "; + std::cin >> threshhold; + std::cout << "New Sqrt factor: "; + std::cin >> sqrtf; + std::cout << "New # of peaks: "; + std::cin >> nPeaks; + goto doagainfit; + } + if(peakscheck == "goback"){ + a=0; + std::cout << "How many bars? "; + std::cin >> a; + } + while(peakscheck == 'c'){ + cout<< "We see "<> peakscheck; + } + + nfp = 0; + for (Int_t ii = 0; ii < mPeaks; ii++) + { + std::cout << "Peak @ " << peaks[ii]; + while ((std::cout << " corresponds to Z=") && !(std::cin >> z)) + { + std::cout << "That's not a number;"; + std::cin.clear(); + std::cin.ignore(std::numeric_limits::max(), '\n'); + } + if (z == 0) + continue; + x[nfp] = (Double_t)z; + zpeaks[nfp] = peaks[ii]; + // std::cout<<"z real " << x[nfp] << " z found " << zpeaks[nfp] <<"\n"; + nfp++; + } + doagain = 'n'; +// std::cout << "Do again? (y/n) "; +// std::cin >> doagain; + } while (doagain != "n"); + + if (mPeaks < 2) + { + pars[0] = 0.; + //pars[1] = fMaxQ / peaks[0]; + pars[1] = z / peaks[0]; + pars[2] = 0.; + + delete s; + delete c1; + delete histo; + delete h1; + continue; + } + + // fit charge axis + std::cout << "Selected " << nfp << " useful peaks to fit\nStart fitting...\n"; + c1->cd(3)->Clear(); + c1->Update(); + c1->cd(3); + auto gr1 = new TGraph(nfp, zpeaks, x); + gr1->SetMarkerColor(4); + gr1->SetMarkerSize(1.5); + gr1->SetMarkerStyle(20); + gr1->Draw("AP"); + //if(fitz) delete fitz; + //TF1* fitz = new TF1("fitz", "[0] + [1]*x + [2]*x*x", min, max); + if (fZfitType != "pol1" && fZfitType != "pol2") + { + R3BLOG(error, "Fit " << fZfitType << " is not allowed, use pol1 or pol2 "); + continue; + } + auto fitz = new TF1("fitz", /*fZfitType*/"pol1", 0., 12.); + fitz->SetLineColor(2); + fitz->SetLineWidth(2); + fitz->SetLineStyle(1); + // fitz->SetParameters(1.5, 2., .1); + gr1->Fit(fitz); + fitz->Draw("lsame"); + c1->Update(); + + // write parameters + std::cout << "Is OK? (y/n) "; + std::cin >> doagain; + if (doagain == "n"){ + std::cout << "New sigma: "; + std::cin >> sigma; + std::cout << "New threshhold: "; + std::cin >> threshhold; + delete s; + goto doagainfit; + } + + int nbpars = 2; + if (fZfitType == "pol2") + nbpars = 3; + + pars[0] = fitz->GetParameter(0); + pars[1] = fitz->GetParameter(1); + + for (Int_t k = 0; k < nbpars; k++) + { + pars[k] = fitz->GetParameter(k); + std::cout<GetParameter(k)<<"\n"; + //std::cout<Write(); + delete s; + delete histo; + delete h1; + delete gr1; + delete c1; + delete fitz; + + //end zcorr + + + std::cout << "Write parameter: " << pars[0] << " " << pars[1] << " " << pars[2] << "\n"; + + + if(a > 0){ + j = j - a; + continue; + } + + + Double_t offset1 = parY->GetOffset1(); + Double_t offset2 = parY->GetOffset2(); + Double_t veff = parY->GetVeff(); + Double_t sync = parY->GetSync(); + para[0] = parY->GetPar1za(); + para[1] = parY->GetPar1zb(); + para[2] = parY->GetPar1zc(); + + std::cout << "Old parameter: " << para[0] << " " << para[1] << " " << para[2] << "\n"; + std::cout << "New parameter: " << pars[0] << " " << pars[1] << " " << pars[2] << "\n"; + std::cout << "Write parameter: " << para[0]*pars[1]+pars[0] << " " << pars[1]*para[1] << " " << para[2]*pars[1] << "\n"; + + parY->SetPar1za(para[0]*pars[1]+pars[0]); + parY->SetPar1zb(pars[1]*para[1]); + parY->SetPar1zc(para[2]*pars[1]); + } + + } + + } + + } + + /*if (fParameter == 5) + { + // Determine effective speed of light in [cm/s] for each paddle + LOG(info) << "Calling function"; + calcZoffset(); + // Determine light attenuation lambda for each paddle + //LOG(info) << "Calling function calcLambda"; + //calcLambda(fTofdTotLow, fTofdTotHigh); + } +*/ +// for (Int_t i = 0; i < N_TOFD_HIT_PLANE_MAX; i++) +// { + //if(fhlos[i]) + //fhlos[i]->Write(); + // if(fhTofsync[i]) + // fhTofsync[i]->Write(); + // if (fh_tofd_TotPm[i]) + // fh_tofd_TotPm[i]->Write(); // control histogram for ToT + // if (fhTsync[i]) + // fhTsync[i]->Write(); // histogram for sync calculation + // if (fhTdiff[i]) + // fhTdiff[i]->Write(); // histogram for offset and veff calculation + // for (Int_t j = 0; j < N_TOFD_HIT_PADDLE_MAX; j++) + // { + // if(fhlost[1][11]) + // fhlost[1][11]->Write(); + // if(fhtott[i][j]) + // fhtott[i][j]->Write(); + // if (fh1_tofsync[i][j]) + // fh1_tofsync[i][j]->Write(); // histogram for ToF sync calculation + // if (fhLogTot1vsLogTot2[i][j]) + // fhLogTot1vsLogTot2[i][j]->Write(); // control histogram Log(ToT) Pm1 vs Log(ToT) Pm2 + // if (fhSqrtQvsPosToT[i][j]) + // fhSqrtQvsPosToT[i][j]->Write(); // histogram for ToT offset calculation + // if (fhQvsPos[i][j]) + // fhQvsPos[i][j]->Write(); // histogram for charge fit + // if (fh1_walk_b[i][j]) + // fh1_walk_b[i][j]->Write(); + // if (fh1_walk_t[i][j]) + // fh1_walk_t[i][j]->Write(); + //if (fZoffset[i][j]) + // fZoffset[i][j]->Write(); + + + + // if (!fTofdSmiley) + // { + // if (fhTot1vsPos[i][j]) + // fhTot1vsPos[i][j]->Write(); // histogram for position dependence of charge 1 + // if (fhTot2vsPos[i][j]) + // fhTot2vsPos[i][j]->Write(); // histogram for position dependence of charge 2 + // } + /* + if (fhTot1vsTot2[i][j]) + fhTot1vsTot2[i][j]->Write(); // control histogram ToT Pm1 vs ToT Pm2 + */ + //} + //} } Double_t R3BTofDCal2HitPar::walk(Double_t Q, Double_t par1, Double_t par2, Double_t par3, Double_t par4, Double_t par5) { - Double_t y = 0; - y = -30.2 + par1 * TMath::Power(Q, par2) + par3 / Q + par4 * Q + par5 * Q * Q; - return y; + Double_t y = 0; + y = -30.2 + par1 * TMath::Power(Q, par2) + par3 / Q + par4 * Q + par5 * Q * Q; + return y; } Double_t R3BTofDCal2HitPar::saturation(Double_t x) { - Double_t kor; - Int_t voltage = 600; - if (voltage == 600) - { - if (x < 173) - { - kor = 0.; - } - else if (x > 208) - { - kor = -1.73665e+03 + 2.82009e+01 * 208. - 1.53846e-01 * (208. * 208.) + 2.82425e-04 * (208. * 208. * 208.); - } - else - { - kor = -1.73665e+03 + 2.82009e+01 * x - 1.53846e-01 * (x * x) + 2.82425e-04 * (x * x * x); - } - } - if (voltage == 500) - { - if (x < 95.5) - { - kor = 0.; - } - else if (x > 124) - { - kor = 1.08 * x - 112.44; - } - else - { - kor = 643.257 - 16.7823 * x + 0.139822 * (x * x) - 0.000362154 * (x * x * x); - } - } - if (voltage == 700) - { - if (x < 198) - { - kor = 0.; - } - else if (x > 298) - { - kor = 0.21 * x - 45.54; - } - else - { - kor = -19067 + 383.93 * x - 3.05794 * (x * x) + 0.0120429 * (x * x * x) - 2.34619e-05 * (x * x * x * x) + - 1.81076e-08 * (x * x * x * x * x); - } - } - return kor; + Double_t kor; + Int_t voltage = 600; + if (voltage == 600) + { + if (x < 173) + { + kor = 0.; + } + else if (x > 208) + { + kor = -1.73665e+03 + 2.82009e+01 * 208. - 1.53846e-01 * (208. * 208.) + 2.82425e-04 * (208. * 208. * 208.); + } + else + { + kor = -1.73665e+03 + 2.82009e+01 * x - 1.53846e-01 * (x * x) + 2.82425e-04 * (x * x * x); + } + } + if (voltage == 500) + { + if (x < 95.5) + { + kor = 0.; + } + else if (x > 124) + { + kor = 1.08 * x - 112.44; + } + else + { + kor = 643.257 - 16.7823 * x + 0.139822 * (x * x) - 0.000362154 * (x * x * x); + } + } + if (voltage == 700) + { + if (x < 198) + { + kor = 0.; + } + else if (x > 298) + { + kor = 0.21 * x - 45.54; + } + else + { + kor = -19067 + 383.93 * x - 3.05794 * (x * x) + 0.0120429 * (x * x * x) - 2.34619e-05 * (x * x * x * x) + + 1.81076e-08 * (x * x * x * x * x); + } + } + return kor; } ClassImp(R3BTofDCal2HitPar); diff --git a/tofd/pars/R3BTofDCal2HitPar.h b/tofd/pars/R3BTofDCal2HitPar.h index 99aa398d3..41cedbac0 100644 --- a/tofd/pars/R3BTofDCal2HitPar.h +++ b/tofd/pars/R3BTofDCal2HitPar.h @@ -1,6 +1,6 @@ /****************************************************************************** * Copyright (C) 2019 GSI Helmholtzzentrum für Schwerionenforschung GmbH * - * Copyright (C) 2019-2023 Members of R3B Collaboration * + * Copyright (C) 2019 Members of R3B Collaboration * * * * This software is distributed under the terms of the * * GNU General Public Licence (GPL) version 3, * @@ -161,6 +161,10 @@ class R3BTofDCal2HitPar : public FairTask * Method for setting the upper range of ToT for offset calibration */ void SetTofdTotHigh(Double_t TotHigh) { fTofdTotHigh = TotHigh; } + + void SetTotRange(Double_t totmin, Double_t totmax) {fToTMin = totmin; fToTMax = totmax;} + + inline void ReadHistoFile(TString file) { fHistoFile = file; } private: R3BCoarseTimeStitch* fTimeStitch; @@ -168,7 +172,8 @@ class R3BTofDCal2HitPar : public FairTask * Method for creating histograms. */ void CreateHistograms(Int_t, Int_t); - + void SetHistogramsNULL(Int_t, Int_t); +void calcZoffset(); void calcOffset(); void calcSync(); void calcVeff(); @@ -181,7 +186,7 @@ class R3BTofDCal2HitPar : public FairTask /** * Method for calculation of z correction. */ - void zcorr(TH2F* histo, Int_t min, Int_t max, Double_t*, Int_t index); + int zcorr(TH2F* histo, Double_t min, Double_t max, Double_t*, Int_t, Int_t, Int_t index); /** * Method for calculation of ToT offset. @@ -193,6 +198,7 @@ class R3BTofDCal2HitPar : public FairTask */ Double_t walk(Double_t Q, Double_t par1, Double_t par2, Double_t par3, Double_t par4, Double_t par5); + /** * Method for calculation of saturation. */ @@ -216,7 +222,10 @@ class R3BTofDCal2HitPar : public FairTask R3BEventHeader* fHeader; /* Event header */ Double_t fTofdY; Double_t fTofdQ; + Double_t fMinQ; Double_t fMaxQ; + Double_t fToTMin; + Double_t fToTMax; Int_t fNbZPeaks; Bool_t fTofdSmiley; Bool_t fTofdZ; @@ -224,19 +233,29 @@ class R3BTofDCal2HitPar : public FairTask Double_t fTofdTotLow; Double_t fTofdTotHigh; Double_t fMeanTof; + TString fHistoFile; // Arrays of control histograms - TH2F* fh_tofd_TotPm[N_TOFD_HIT_PLANE_MAX]; + TH2F* fh_tofd_TotPm[N_TOFD_HIT_PLANE_MAX][3]; + TH2F* fh_tofd_TotPm1[N_TOFD_HIT_PADDLE_MAX]; + TH2F* fh_tofd_TotPm2[N_TOFD_HIT_PADDLE_MAX]; TH2F* fhTdiff[N_TOFD_HIT_PLANE_MAX]; TH2F* fhTsync[N_TOFD_HIT_PLANE_MAX]; TH1F* fh1_tofsync[N_TOFD_HIT_PLANE_MAX][N_TOFD_HIT_PADDLE_MAX]; TH2F* fhLogTot1vsLogTot2[N_TOFD_HIT_PLANE_MAX][N_TOFD_HIT_PADDLE_MAX]; TH2F* fhSqrtQvsPosToT[N_TOFD_HIT_PLANE_MAX][N_TOFD_HIT_PADDLE_MAX]; TH2F* fhQvsPos[N_TOFD_HIT_PLANE_MAX][N_TOFD_HIT_PADDLE_MAX]; - // TH2F* fhTot1vsTot2[N_TOFD_HIT_PLANE_MAX][N_TOFD_HIT_PADDLE_MAX]; + TH2F* fhQXvsPos[N_TOFD_HIT_PLANE_MAX][N_TOFD_HIT_PADDLE_MAX]; + TH2F* fhTot1vsTot2[N_TOFD_HIT_PLANE_MAX][N_TOFD_HIT_PADDLE_MAX]; TH2F* fhTot1vsPos[N_TOFD_HIT_PLANE_MAX][N_TOFD_HIT_PADDLE_MAX]; TH2F* fhTot2vsPos[N_TOFD_HIT_PLANE_MAX][N_TOFD_HIT_PADDLE_MAX]; - + TH2F* fhTofsync[N_TOFD_HIT_PLANE_MAX]; + TH2F* fhlos[N_TOFD_HIT_PLANE_MAX]; + TH2F* fhlost[N_TOFD_HIT_PLANE_MAX][N_TOFD_HIT_PADDLE_MAX]; + TH2F* fhtott[N_TOFD_HIT_PLANE_MAX][N_TOFD_HIT_PADDLE_MAX]; + TH2F* fh1_walk_t[N_TOFD_HIT_PLANE_MAX][N_TOFD_HIT_PADDLE_MAX]; + TH2F* fh1_walk_b[N_TOFD_HIT_PLANE_MAX][N_TOFD_HIT_PADDLE_MAX]; + TH2F* fZoffset[N_TOFD_HIT_PLANE_MAX][N_TOFD_HIT_PADDLE_MAX]; public: ClassDef(R3BTofDCal2HitPar, 1) };