BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
Tof400DigitAnalysis_period6.C
Go to the documentation of this file.
1/*
2 * To change this license header, choose License Headers in Project Properties.
3 * To change this template file, choose Tools | Templates
4 * and open the template in the editor.
5 */
6
7#include "TCanvas.h"
8#include "TChain.h"
9#include "TClonesArray.h"
10#include "TDirectory.h"
11#include "TFile.h"
12#include "TH1.h"
13#include "TH1F.h"
14#include "TH1S.h"
15#include "TH2F.h"
16#include "TH2S.h"
17#include "TList.h"
18#include "TPad.h"
19#include "TVector3.h"
20
21#include <TMath.h>
22#include <ctime>
23#include <fstream>
24#include <iostream>
25#include <istream>
26#include <stdio.h>
27
29
30Double_t CorrLR[10][48] = {0.};
31
34
35class BmnTrigDigit;
36class BmnTOF1Detector;
37
38void Tof400DigitAnalysis_period6(TString file = "", Int_t nEvForRead = 0, Int_t Periud = 5)
39{
40 TStopwatch timer;
41 timer.Start();
42
43 gROOT->LoadMacro("$VMCWORKDIR/macro/run/bmnloadlibs.C");
44 bmnloadlibs(); // load BmnRoot libraries
45
46 ifstream f_call;
47 char line[256];
48 Int_t Pl, St;
49 Double_t Temp;
50 TString NameCallFile = "Tof400LRcorr.dat";
51 f_call.open(NameCallFile.Data());
52 f_call.getline(line, 256);
53 f_call.getline(line, 256);
54 while (!f_call.eof()) {
55 f_call >> Pl >> St >> Temp >> CorrLR[Pl][St];
56 if (TMath::Abs(Temp - CorrLR[Pl][St]) > 2.)
57 CorrLR[Pl][St] = -11.9766;
58 }
59
60 BmnTOF1Detector* Plane[10];
61
62 TString name;
63 for (Int_t i = 0; i < 10; i++) {
64 name = Form("Plane%d", i);
65 Plane[i] = new BmnTOF1Detector(i, 2);
66 // Plane[i]->SetCorrLR(CorrLR[i]);
67 Plane[i]->SetCorrLR("Tof400LRcorr.dat");
68 Plane[i]->SetCorrSlewing("Tof400SlewingCorr_period6.root");
69 Plane[i]->SetGeoFile("geofile_full.root");
70 // Plane[i]->SetGeoFile("geometry_run6.root");
71 }
72
73 TList* fList = new TList();
74 TList* fListTof = new TList();
75 TH1I* h_dtBDT0 = new TH1I("h_dtBDT0", "h_dtBDT0", 1024, -24, 24);
76 fList->Add(h_dtBDT0);
77 TH2F* h_dtBDT0_vs_AmpT0 = new TH2F("h_dtBDT0_vs_AmpT0", "h_dtBDT0_vs_AmpT0", 1024, 0, 48, 1024, -24, 24);
78 fList->Add(h_dtBDT0_vs_AmpT0);
79 TH2F* h_dtBDT0_vs_AmpBD = new TH2F("h_dtBDT0_vs_AmpBD", "h_dtBDT0_vs_AmpBD", 1024, 0, 48, 1024, -24, 24);
80 fList->Add(h_dtBDT0_vs_AmpBD);
81 TH2F* h_TimeBD_vs_TimeT0 = new TH2F("h_TimeBD_vs_TimeT0", "h_TimeBD_vs_TimeT0", 680, 250, 420, 680, 250, 420);
82 fList->Add(h_TimeBD_vs_TimeT0);
83 TH1I* h_dtBDT0_Mod = new TH1I("h_dtBDT0_Mod", "h_dtBDT0_Mod", 1024, -24, 24);
84 fList->Add(h_dtBDT0_Mod);
85 TH1I* h_NHitBD = new TH1I("h_NHitBD", "h_NHitBD", 40, 0, 40);
86 fList->Add(h_NHitBD);
87
88 TH1I* h_dtBC2T0 = new TH1I("h_dtBC2T0", "h_dtBC2T0", 1024, -6., 6.);
89 fList->Add(h_dtBC2T0);
90 TH2F* h_dtBC2T0_vs_AmpT0 = new TH2F("h_dtBC2T0_vs_AmpT0", "h_dtBC2T0_vs_AmpT0", 1024, 0, 48, 1024, -12., 12.);
91 fList->Add(h_dtBC2T0_vs_AmpT0);
92 TH2F* h_dtBC2T0_vs_AmpBC2 = new TH2F("h_dtBC2T0_vs_AmpBC2", "h_dtBC2T0_vs_AmpBC2", 1024, 0, 48, 1024, -12., 12.);
93 fList->Add(h_dtBC2T0_vs_AmpBC2);
94 TH2F* h_TimeBC2_vs_TimeT0 = new TH2F("h_TimeBC2_vs_TimeT0", "h_TimeBC2_vs_TimeT0", 680, 250, 420, 680, 250, 420);
95 fList->Add(h_TimeBC2_vs_TimeT0);
96 TH1I* h_dtBC2T0_Corr = new TH1I("h_dtBC2T0_Corr", "h_dtBC2T0_Corr", 1024, -6., 6.);
97 fList->Add(h_dtBC2T0_Corr);
98
99 TH2I* hCrosing = new TH2I("hCrosing", "hCrosing", 48, 0, 48, 48, 0, 48);
100 fList->Add(hCrosing);
101 TH1I* hCrosingDx = new TH1I("hCrosingDx", "hCrosingDx", 48, 0., 60.);
102 fList->Add(hCrosingDx);
103 TH1I* hCrosingDy = new TH1I("hCrosingDy", "hCrosingDy", 120, -30., 30.);
104 fList->Add(hCrosingDy);
105 TH1I* hCrosingDr = new TH1I("hCrosingDr", "hCrosingDr", 600, 0., 60.);
106 fList->Add(hCrosingDr);
107 TH1I* hCrosingDt = new TH1I("hCrosingDt", "hCrosingDt", 1024, -12., 12.);
108 fList->Add(hCrosingDt);
109 TH2I* hCrosingWidth_Width = new TH2I("hCrosingWidth_Width", "hCrosingWidth_Width", 1024, 12, 60, 1024, 12, 60);
110 fList->Add(hCrosingWidth_Width);
111 TH1I* hCrosingDwidth = new TH1I("hCrosingDwidth", "hCrosingDwidth", 512, -12., 12.);
112 fList->Add(hCrosingDwidth);
113
114 TH2F* h_HitRPC = new TH2F("h_HitRPC", "h_HitRPC", 96, -96, 0, 5, 0, 5);
115 fListTof->Add(h_HitRPC);
116
117 TH2F* h_XYRPC = new TH2F("h_XYRPC", "h_XYRPC", 240, -150, 150, 120, -75, 75);
118 fListTof->Add(h_XYRPC);
119
120 TChain* eveTree = new TChain("cbmsim");
121 TString inName = Form("/home/storage/digi/%s", file.Data());
122 // inName = file;
123 cout << "Open file " << inName << endl << endl;
124 eveTree->Add(inName);
125
126 TClonesArray* ToF400Digits;
127 eveTree->SetBranchAddress("TOF400", &ToF400Digits);
128
129 TClonesArray* T0Digits;
130 eveTree->SetBranchAddress("T0", &T0Digits);
131
132 TClonesArray* VetoDigits;
133 eveTree->SetBranchAddress("VETO", &VetoDigits);
134
135 TClonesArray* BDDigits;
136 eveTree->SetBranchAddress("BD", &BDDigits);
137
138 TClonesArray* BC2Digits;
139 eveTree->SetBranchAddress("BC2", &BC2Digits);
140
141 TClonesArray* EventHeader;
142 eveTree->SetBranchAddress("EventHeader", &EventHeader);
143
144 Long64_t nEvents = eveTree->GetEntries();
145 if (nEvForRead == 0 || nEvForRead > nEvents)
146 nEvForRead = nEvents;
147 cout << "Will be readed " << nEvForRead << " events from " << nEvents << endl;
148
149 for (Int_t iEv = 0; iEv < nEvForRead; iEv++) {
150
151 if (iEv % 50000 == 0)
152 cout << "EVENT: " << iEv << endl;
153 eveTree->GetEntry(iEv);
154 if ((T0Digits->GetEntriesFast()) == 1 && (BC2Digits->GetEntriesFast()) == 1
155 && (VetoDigits->GetEntriesFast()) == 0 && (BDDigits->GetEntriesFast()) >= 2
156 //&& (BDDigits->GetEntriesFast()) == 0
157 )
158 {
159 BmnTrigDigit* digT0 = (BmnTrigDigit*)T0Digits->At(0);
160 if (digT0->GetAmp() >= 17.3 && digT0->GetAmp() <= 19.2) {
161 AmpBD = 0.;
162 TimeBD = 0.;
163 NHitBD = 0.;
164
165 BmnTrigDigit* digT0 = (BmnTrigDigit*)T0Digits->At(0);
166
167 //--------------------------- EventHeader --------------------------------------------------
168 Int_t iEvDig = EventHeader->GetEntriesFast();
169 if (iEvDig != 1)
170 cout << "iEvDig == " << iEvDig << endl;
171 BmnEventHeader* digEvent = (BmnEventHeader*)EventHeader->At(0);
172
173 //--------------------------- RPC --------------------------------------------------
174 for (Int_t i = 0; i < 10; i++)
175 Plane[i]->Clear();
176
177 for (Int_t iDig = 0; iDig < ToF400Digits->GetEntriesFast(); ++iDig) {
178 FlagHit = kFALSE;
179 BmnTof1Digit* digTof = (BmnTof1Digit*)ToF400Digits->At(iDig);
180 FlagHit = Plane[digTof->GetPlane()]->SetDigit(digTof);
181 if (FlagHit == kTRUE) {
182 if (digTof->GetPlane() < 5) {
183 XforHit = digTof->GetStrip() * -1. - 48;
184 YforHit = 4 - digTof->GetPlane();
185 } else {
186 XforHit = digTof->GetStrip() * -1.;
187 YforHit = 4 - (digTof->GetPlane() - 5);
188 }
189 h_HitRPC->Fill(XforHit, YforHit);
190 }
191 }
192
193 TVector3 XYZ;
194 XYZ.SetXYZ(0., 0., 0.);
195 Double_t ToF = 0;
196 FlagHit = kFALSE;
197 for (Int_t i = 0; i < 10; i++) {
198 Plane[i]->FindHits(digT0);
199 for (Int_t s = 0; s < 47; s++) {
200 XYZ.SetXYZ(0., 0., 0.);
201 ToF = 0;
202 FlagHit = kFALSE;
203 FlagHit = Plane[i]->GetXYZTime(s, &XYZ, &ToF);
204 if (FlagHit == kTRUE)
205 h_XYRPC->Fill(XYZ.x(), XYZ.y());
206 }
207 }
208 CrossingDet(Plane[2], Plane[7], hCrosing, hCrosingDx, hCrosingDy, hCrosingDr, hCrosingDt);
209
210 //--------------------------- BD --------------------------------------------------
211 for (Int_t iDig = 0; iDig < BDDigits->GetEntriesFast(); ++iDig) {
212 BmnTrigDigit* digBD = (BmnTrigDigit*)BDDigits->At(iDig);
213 h_TimeBD_vs_TimeT0->Fill(digT0->GetTime(), digBD->GetTime());
214 if (TMath::Abs(digBD->GetTime() - digT0->GetTime() - 13.) < 1.8) {
215 AmpBD = AmpBD + digBD->GetAmp();
216 TimeBD = TimeBD + digBD->GetTime();
217 NHitBD++;
218 }
219 h_dtBDT0_Mod->Fill(CalculateDtT0BD(digBD, digT0));
220 }
221 h_NHitBD->Fill(NHitBD);
222 AmpBD = AmpBD / NHitBD;
223 TimeBD = TimeBD / NHitBD;
224
225 Double_t dtBDT0 = CalculateDtT0BD(TimeBD, AmpBD, digT0);
226 h_dtBDT0->Fill(dtBDT0);
227 h_dtBDT0_vs_AmpT0->Fill(digT0->GetAmp(), dtBDT0);
228 h_dtBDT0_vs_AmpBD->Fill(AmpBD, dtBDT0);
229
230 //--------------------------- BC2 --------------------------------------------------
231 BmnTrigDigit* digBC2 = (BmnTrigDigit*)BC2Digits->At(0);
232 if (digBC2->GetAmp() >= 16.2 && digBC2->GetAmp() <= 19.3) {
233 Double_t dtBC2T0 = CalculateDtT0BC2(digBC2, digT0);
234 h_dtBC2T0->Fill(digBC2->GetTime() - digT0->GetTime());
235 h_dtBC2T0_Corr->Fill(dtBC2T0);
236 h_dtBC2T0_vs_AmpT0->Fill(digT0->GetAmp(), dtBC2T0);
237 h_dtBC2T0_vs_AmpBC2->Fill(digBC2->GetAmp(), dtBC2T0);
238 h_TimeBC2_vs_TimeT0->Fill(digBC2->GetTime(), digT0->GetTime());
239 } //*/
240
241 } // end if (digT0->GetAmp() >= 17.3 && digT0->GetAmp() <= 19.2)
242 } // end if ((T0Digits->GetEntriesFast()) == 1 && (VetoDigits->GetEntriesFast()) == 0 &&
243 // (BDDigits->GetEntriesFast()) >= 2)
244 } // end for (Int_t iEv = 0; iEv < nEvForRead; iEv++)
245
246 cout << "End Data Analysis" << endl;
247
248 timer.Stop();
249 cout << "Readed " << nEvForRead << " events" << endl;
250 cout << "Time = " << timer.RealTime() << " s" << endl;
251 cout << "Time/Event = " << timer.RealTime() / (Double_t)nEvForRead * 1000. << " ms/Event" << endl;
252 cout << "TimeCPU = " << timer.CpuTime() << " s" << endl;
253 cout << "TimeCPU/Event = " << timer.CpuTime() / (Double_t)nEvForRead * 1000. << " ms/Event" << endl;
254
255 ofstream f_time;
256 TString NameTimeFile = "Tof400An_Time.dat";
257 // TString NameCallFile = file;
258 // Point = NameCallFile.First('.');
259 // NameCallFile.Replace(Point, 15, "_LRcorr.dat");
260 cout << "Write time to " << NameTimeFile.Data() << endl;
261 f_time.open(NameTimeFile.Data(), std::ofstream::out | std::ofstream::app);
262 f_time << file << "\t" << timer.RealTime() << "\t" << timer.RealTime() / (Double_t)nEvForRead * 1000. << "\t"
263 << timer.CpuTime() << "\t" << timer.CpuTime() / (Double_t)nEvForRead * 1000. << endl;
264 f_time.close();
265
266 timer.Reset();
267 timer.Start();
268 TString outName = file;
269 Int_t Point = outName.First('.');
270 outName.Replace(Point, 15, "_TofAn6.root");
271 outName = Form("/home/storage/analysis/%s", outName.Data());
272 cout << "Save data to " << outName.Data() << endl;
273 TFile* fileout = new TFile(outName.Data(), "RECREATE");
274
275 TDirectory* Dir;
276 Dir = fileout->mkdir("ToF");
277 Dir->cd();
278 fListTof->Write();
279
280 TDirectory* DirPlane[10];
281 TDirectory* Dir1Plane[10];
282 TDirectory* Dir2Plane[10];
283 for (Int_t i = 0; i < 10; i++) {
284 DirPlane[i] = Dir->mkdir(Plane[i]->GetName());
285 DirPlane[i]->cd();
286 Plane[i]->GetList(0)->Write();
287 Dir1Plane[i] = DirPlane[i]->mkdir("Detector");
288 Dir1Plane[i]->cd();
289 Plane[i]->GetList(1)->Write();
290 Dir2Plane[i] = DirPlane[i]->mkdir("ToF");
291 Dir2Plane[i]->cd();
292 Plane[i]->GetList(2)->Write();
293 }
294
295 fileout->cd();
296 fList->Write();
297 fileout->Close();
298
299 timer.Stop();
300
301 cout << "Time for write root = " << timer.RealTime() << " s" << endl;
302
303} // end of macros
304
305//----Additional functions ----------------------------------------------------------------
306
308{
309 Double_t dt = bd->GetTime() - t0->GetTime();
310 dt = dt
311 - (14.55 + 0.1218 * bd->GetAmp() - 0.02473 * bd->GetAmp() * bd->GetAmp()
312 + 0.0006237 * bd->GetAmp() * bd->GetAmp() * bd->GetAmp());
313 dt = dt
314 - (1.947 - 0.5363 * t0->GetAmp() + 0.03428 * t0->GetAmp() * t0->GetAmp()
315 - 0.0005853 * t0->GetAmp() * t0->GetAmp() * t0->GetAmp());
316 /* dt = dt - (-3.412 + 0.1866 * t0->GetAmp());//*/
317 return dt;
318}
319
321{
322 Double_t dt = timeBD - t0->GetTime();
323 dt = dt - (14.55 + 0.1218 * ampBD - 0.02473 * ampBD * ampBD + 0.0006237 * ampBD * ampBD * ampBD);
324 dt = dt
325 - (1.947 - 0.5363 * t0->GetAmp() + 0.03428 * t0->GetAmp() * t0->GetAmp()
326 - 0.0005853 * t0->GetAmp() * t0->GetAmp() * t0->GetAmp());
327 // dt = dt - (-3.412 + 0.1866 * t0->GetAmp());//*/
328 return dt;
329}
330
332{
333 Double_t dt = bc2->GetTime() - t0->GetTime();
334 dt = dt
335 - (40.47 - 6.338 * bc2->GetAmp() + 0.35 * bc2->GetAmp() * bc2->GetAmp()
336 - 0.006571 * bc2->GetAmp() * bc2->GetAmp() * bc2->GetAmp());
337 dt = dt
338 - (1.947 - 0.5363 * t0->GetAmp() + 0.03428 * t0->GetAmp() * t0->GetAmp()
339 - 0.0005853 * t0->GetAmp() * t0->GetAmp() * t0->GetAmp());
340 /* dt = dt - (0.4034 - 0.01599 * t0->GetAmp());
341 dt = dt + 0.001108; //zero for C 3.5 Gec/nuclon */
342 return dt;
343}
344
345void CrossingDet(BmnTOF1Detector* Det1, BmnTOF1Detector* Det2, TH2* hD1D2, TH1* hDx, TH1* hDy, TH1* hDr, TH1* hDt)
346{
347 TVector3 XyzD1, XyzD2;
348 Double_t ToF1 = 0, ToF2 = 0, dT = 0.;
349 for (Int_t iD1 = 0; iD1 < 12; iD1++) {
350 XyzD1.SetXYZ(0., 0., 0.);
351 ToF1 = 0;
352 if (Det1->GetXYZTime(iD1, &XyzD1, &ToF1))
353 for (Int_t iD2 = 36; iD2 < 48; iD2++) {
354 XyzD2.SetXYZ(0., 0., 0.);
355 ToF2 = 0.;
356 if (Det2->GetXYZTime(iD2, &XyzD2, &ToF2)) {
357 hD1D2->Fill(iD1, iD2);
358 Double_t dX = XyzD1.x() - XyzD2.x();
359 Double_t dY = XyzD1.y() - XyzD2.y();
360 Double_t dR = TMath::Sqrt(dX * dX + dY * dY);
361 hDx->Fill(dX);
362 hDy->Fill(dY);
363 hDr->Fill(dR);
364 hDt->Fill(ToF1 - ToF2);
365 hCrosingWidth_Width->Fill(Det1->GetWidth(iD1), Det2->GetWidth(iD2));
366 hCrosingDwidth->Fill(Det1->GetWidth(iD1) - Det2->GetWidth(iD2));
367 // cout << 3 << endl;
368 }
369 }
370 }
371}
int i
Definition P4_F32vec4.h:22
bool Bool_t
int Int_t
constexpr Bool_t kFALSE
double Double_t
long long Long64_t
constexpr Bool_t kTRUE
#define gROOT
char * Form(const char *fmt,...)
Bool_t GetXYZTime(Int_t Str, TVector3 *XYZ, Double_t *ToF)
Bool_t SetGeoFile(TString NameFile)
Bool_t SetCorrLR(TString NameFile)
Double_t GetWidth(Int_t Str)
TList * GetList(Int_t n)
Bool_t SetCorrSlewing(TString NameFile)
Short_t GetPlane() const
Short_t GetStrip() const
Double_t GetAmp() const
Double_t GetTime() const
Int_t SetBranchAddress(const char *bname, T **add, TBranch **ptr=nullptr)
virtual Int_t Add(const char *name, Long64_t nentries=TTree::kMaxEntries)
Long64_t GetEntries() const override
Int_t GetEntry(Long64_t entry=0, Int_t getall=0) override
Int_t Write(const char *name=nullptr, Int_t option=0, Int_t bufsize=0) const override
Bool_t cd() override
TDirectory * mkdir(const char *name, const char *title="", Bool_t returnExistingDirectory=kFALSE) override
virtual Bool_t cd()
virtual TDirectory * mkdir(const char *name, const char *title="", Bool_t returnExistingDirectory=kFALSE)
void Close(Option_t *option="") override
virtual Int_t Fill(const char *name, Double_t w)
virtual Int_t Fill(const char *namex, const char *namey, Double_t w)
void Add(TObject *obj) override
Int_t GetEntriesFast() const
TObject * At(Int_t idx) const override
Double_t RealTime()
void Start(Bool_t reset=kTRUE)
Double_t CpuTime()
void Stop()
void Reset()
TString & Replace(Ssiz_t pos, Ssiz_t n, const char *s)
Ssiz_t First(char c) const
const char * Data() const
void SetXYZ(Double_t x, Double_t y, Double_t z)
Double_t x() const
Double_t y() const
Double_t NHitBD
Double_t TimeBD
Double_t AmpBD
Double_t CorrLR[10][48]
void CrossingDet(BmnTOF1Detector *Det1, BmnTOF1Detector *Det2, TH2 *hD1D2, TH1 *hDx, TH1 *hDy, TH1 *hDr, TH1 *hDt)
void Tof400DigitAnalysis_period6(TString file="", Int_t nEvForRead=0, Int_t Periud=5)
Double_t CalculateDtT0BC2(BmnTrigDigit *bc2, BmnTrigDigit *t0)
Double_t CalculateDtT0BD(BmnTrigDigit *bd, BmnTrigDigit *t0)
Double_t Sqrt(Double_t x)
Double_t Abs(Double_t d)