BmnRoot
Loading...
Searching...
No Matches
BmnToCbmHitConverter.cxx
Go to the documentation of this file.
1
2#include "BmnToCbmHitConverter.h"
3
4#include "BmnVSPHit.h" // AZ-060625
5#include "CbmMCTrack.h" //AZ-291223
6#include "DstEventHeader.h" //AZ-010723
7#include "FairField.h" //AZ-030323
8#include "FairRunAna.h" //AZ-030323
9
10#include <TDatabasePDG.h> //AZ-291223
11#include <TFile.h>
12#include <TH2.h>
13#include <TObjString.h>
14#include <TRandom.h> //AZ-191223
15#include <TStopwatch.h>
16
17static Double_t workTime = 0.0;
18
19// ------ Default constructor -----------------------------------------
20BmnToCbmHitConverter::BmnToCbmHitConverter()
21 : BmnTask("BMN to CBM Hits Converter", 1)
22 , fBmnGemHitsBranchName("BmnGemStripHit")
23 , fBmnSilHitsBranchName("BmnSiliconHit")
24 , fBmnVspHitsBranchName("BmnVSPHit")
25 , fCbmHitsBranchName("StsHit")
26 , fBmnGemHitsArray(nullptr)
27 , fBmnGemUpperClusters(nullptr)
28 , fBmnGemLowerClusters(nullptr)
29 , fBmnSilHitsArray(nullptr)
30 , fBmnSilUpperClusters(nullptr)
31 , fBmnSilLowerClusters(nullptr)
32 , fBmnVspHitsArray(nullptr)
33 , fBmnVspUpperClusters(nullptr)
34 , fBmnVspLowerClusters(nullptr)
35 , fCbmHitsArray(nullptr)
36 , fGemConfigFile("GemRun8.xml")
37 , fSilConfigFile("SiliconRun8.xml")
38 , fVspConfigFile("")
39 , GemStationSet(nullptr)
40 , SilStationSet(nullptr)
41 , fVspStationSet(nullptr)
42 , fUseFixedErrors(kFALSE)
43 , fExp(kFALSE)
44 , fMCcorr(kTRUE)
45 , fieldFlag(kFALSE)
46 , fAlignmentEnabled(kTRUE)
47{}
48// -------------------------------------------------------------------------
49
50// ----- Standard constructor ------------------------------------------
51// AZ-170323 BmnToCbmHitConverter::BmnToCbmHitConverter(Int_t iVerbose)
52BmnToCbmHitConverter::BmnToCbmHitConverter(Int_t iVerbose, Bool_t isExp)
53 : BmnTask("BMN to CBM Hits Converter", iVerbose)
54 , fBmnGemHitsBranchName("BmnGemStripHit")
55 , fBmnSilHitsBranchName("BmnSiliconHit")
56 , fBmnVspHitsBranchName("BmnVSPHit")
57 , fCbmHitsBranchName("StsHit")
58 , fBmnGemHitsArray(nullptr)
59 , fBmnGemUpperClusters(nullptr)
60 , fBmnGemLowerClusters(nullptr)
61 , fBmnSilHitsArray(nullptr)
62 , fBmnSilUpperClusters(nullptr)
63 , fBmnSilLowerClusters(nullptr)
64 , fBmnVspHitsArray(nullptr)
65 , fBmnVspUpperClusters(nullptr)
66 , fBmnVspLowerClusters(nullptr)
67 , fCbmHitsArray(nullptr)
68 , fGemConfigFile("GemRun8.xml")
69 , fSilConfigFile("SiliconRun8.xml")
70 , fVspConfigFile("")
71 , GemStationSet(nullptr)
72 , SilStationSet(nullptr)
73 , fVspStationSet(nullptr)
74 , fUseFixedErrors(kFALSE)
75 , fExp(isExp)
76 , fMCcorr(kTRUE)
77 , fieldFlag(kFALSE)
78 , fAlignmentEnabled(kTRUE)
79{}
80// -------------------------------------------------------------------------
81
82// ----- Destructor ----------------------------------------------------
83BmnToCbmHitConverter::~BmnToCbmHitConverter() {}
84// -------------------------------------------------------------------------
85
86// ----- Public method Exec --------------------------------------------
87void BmnToCbmHitConverter::Exec(Option_t* opt)
88{
89
90 TStopwatch sw;
91 sw.Start();
92
93 if (!IsActive())
94 return;
95
96 fCbmHitsArray->Delete();
97
98 int nVspStat = (fVspStationSet) ? fVspStationSet->GetNStations() : 0; // AZ-070625
99
100 if (fBmnGemHitsArray) {
101
102 for (Int_t iHit = 0; iHit < fBmnGemHitsArray->GetEntriesFast(); ++iHit) {
103 BmnGemStripHit* bmnHit = (BmnGemStripHit*)fBmnGemHitsArray->At(iHit);
104 if (bmnHit->IsPseudo())
105 continue; // AZ-180723 - exclude 1-dim hits (one-sided pseudo hit)
106
107 // if (bmnHit->GetRefIndex() == -1) continue;
108
109 // Section for hit filtration by signal asymmetry
110 // StripCluster* uc = (StripCluster*)fBmnGemUpperClusters->At(bmnHit->GetUpperClusterIndex());
111 // StripCluster* lc = (StripCluster*)fBmnGemLowerClusters->At(bmnHit->GetLowerClusterIndex());
112 // Float_t ls = lc->TotalSignal;
113 // Float_t us = uc->TotalSignal;
114 // if (us < ls - 1000 || us > ls + 1000) continue;
115 // if (Abs((us - ls) / (us + ls)) > 0.75) //asymmetry
116 // continue;
117
118 TVector3 pos;
119 bmnHit->Position(pos);
120 TVector3 dpos;
121 if (fUseFixedErrors) {
122 // dpos[0] = 0.08/TMath::Sqrt(12); //AZ
123 // dpos[1] = 0.1234; //AZ
124 dpos[0] = fDXgem; // 0.015; //AZ - as in cbmroot
125 dpos[1] = fDYgem; // 0.058; //AZ - as in cbmroot
126 } else {
127 bmnHit->PositionError(dpos);
128 }
129
130 Int_t stat = bmnHit->GetStation();
131 Int_t mod = bmnHit->GetModule();
132 // ElectronDriftDirectionInModule driftDir =
133 // GemStationSet->GetStation(stat)->GetModule(mod)->GetElectronDriftDirection(); //AZ-200322
134 // if (driftDir == BackwardZAxisEDrift) pos[2] -= 0.9; //AZ-200322
135 // if (driftDir == ForwardZAxisEDrift) pos[2] -= 0.9; //AZ-200322
136
137 Int_t lay = 0;
138 for (lay = 0; lay < GemStationSet->GetStation(stat)->GetModule(mod)->GetNStripLayers(); lay++) {
139 BmnGemStripLayer* layer = &(GemStationSet->GetStation(stat)->GetModule(mod)->GetStripLayer(lay));
140 if (layer->IsPointInsideStripLayer(-bmnHit->GetX(), bmnHit->GetY()))
141 break;
142 }
143 if (lay == GemStationSet->GetStation(stat)->GetModule(mod)->GetNStripLayers())
144 continue; // AZ-230322 - strange case //CHECK IT!!!
145
146 // formula for converting from the bm@n system of modules and layers to the cbm one
147 Int_t sect = 2 * mod + 1 + lay / 2;
148
149 Int_t sens = 1;
150
151 Int_t nSilStations = SilStationSet->GetNStations();
152 Int_t detId = kGEM << 24 | (stat + 1 + nVspStat + nSilStations) << 16 | sect << 4 | sens << 1;
153
154 // AZ-140226 - remove beam spot in GEMS for mag. field == 0
155 if (!fieldFlag) {
156 if ((stat == 1 && pos[0] >= -16 && pos[0] < -2) || (stat == 2 && pos[0] >= -18 && pos[0] < -2)
157 || (stat == 3 && pos[0] >= -12 && pos[0] < -2) || (stat == 4 && pos[0] >= -10 && pos[0] < 3)
158 || (stat == 5 && pos[0] >= -10 && pos[0] < 3) || (stat == 6 && pos[0] >= -6 && pos[0] < 4))
159 {
160 if (TMath::Abs(pos[1]) < 15)
161 continue;
162 }
163 }
164
165 new ((*fCbmHitsArray)[fCbmHitsArray->GetEntriesFast()]) CbmStsHit(detId, pos, dpos, 0.0, 0, 0);
166 CbmStsHit* hit = (CbmStsHit*)fCbmHitsArray->At(fCbmHitsArray->GetEntriesFast() - 1);
167
168 // AZ-120223 - apply alignment
169 if (fExp && fAlignmentEnabled)
170 ApplyAlignment(hit);
171 else if (fieldFlag && !fExp) {
172 // AZ-140824 - shift coordinates
173 hit->SetY(hit->GetY() + 0.30);
174 if (fMCcorr) {
175 int iq = 0;
176 Double_t pmom = 0;
177 if (bmnHit->GetRefIndex() >= 0) {
178 // Different smearing for negative tracks
179 FairMCPoint* mcp = (FairMCPoint*)fGemPoints->UncheckedAt(bmnHit->GetRefIndex());
180 CbmMCTrack* mctr = (CbmMCTrack*)fMCTracks->UncheckedAt(mcp->GetTrackID());
181 pmom = mctr->GetP();
182 if (TDatabasePDG::Instance()->GetParticle(mctr->GetPdgCode()))
183 iq = TMath::Nint(TDatabasePDG::Instance()->GetParticle(mctr->GetPdgCode())->Charge() / 3);
184 }
185 // Correct MC hit
186 CorrectHitMC(hit, pmom, iq);
187 if (iq == 0)
188 iq = 1; // for backward compatibility
189
190 Bool_t ok = CorrectEffic(hit, iq); // efficiency correction
191 if (!ok) {
192 // Remove hit
193 fCbmHitsArray->RemoveLast();
194 continue;
195 }
196 } // if (fMCcorr)
197 }
198
199 FairRootManager::Instance()->SetUseFairLinks(kTRUE);
200 hit->ResetLinks();
201 hit->SetLinks(bmnHit->GetLinks());
202 FairRootManager::Instance()->SetUseFairLinks(kFALSE);
203 hit->SetRefIndex(bmnHit->GetRefIndex());
204 if (fExp)
205 hit->SetRefIndex(iHit); // AZ-180226 - original hit index (for original coordinates)
206 hit->fDigiF = bmnHit->GetLowerClusterIndex() + 0; // AZ-250322
207 hit->fDigiB = bmnHit->GetUpperClusterIndex() + 1000000; // AZ-250322
208 // AZ-221022
209 StripCluster* upp = (StripCluster*)fBmnGemUpperClusters->UncheckedAt(hit->fDigiB % 1000000);
210 StripCluster* low = (StripCluster*)fBmnGemLowerClusters->UncheckedAt(hit->fDigiF % 1000000);
211 Double_t uppq = upp->TotalSignal;
212 Double_t lowq = low->TotalSignal;
213 Double_t corrq = (uppq - lowq) / (uppq + lowq);
214 hit->SetSignalDiv(corrq);
215 uppq = upp->Width; // N of strips
216 lowq = low->Width;
217 corrq = (uppq - lowq) / (uppq + lowq);
218 hit->SetTimeStamp(corrq); // just using available space
219 hit->SetStrips(low->MeanPosition, upp->MeanPosition); // AZ-190823
220 // hit->SetWidths (low->Width, upp->Width); //AZ-190823
221 } // for (Int_t iHit = 0; iHit < fBmnGemHitsArray->GetEntriesFast();
222 }
223
224 if (fBmnSilHitsArray) {
225 for (Int_t iHit = 0; iHit < fBmnSilHitsArray->GetEntriesFast(); ++iHit) {
226 BmnSiliconHit* bmnHit = (BmnSiliconHit*)fBmnSilHitsArray->At(iHit);
227 if (bmnHit->IsPseudo())
228 continue; // AZ-180723 - exclude 1-dim hits (one-sided pseudo hit)
229 // if (bmnHit->GetRefIndex() == -1) continue;
230 // Section for hit filtration by signal asymmetry
231 // StripCluster* uc = (StripCluster*)fBmnSilUpperClusters->At(bmnHit->GetUpperClusterIndex());
232 // StripCluster* lc = (StripCluster*)fBmnSilLowerClusters->At(bmnHit->GetLowerClusterIndex());
233 // Float_t ls = lc->TotalSignal;
234 // Float_t us = uc->TotalSignal;
235 // if (us < ls - 1000 || us > ls + 1000) continue;
236 // if (Abs((us - ls) / (us + ls)) > 0.75) //asymmetry
237 // continue;
238
239 TVector3 pos;
240 bmnHit->Position(pos);
241 pos[2] -= 0.0150; // AZ - shift to the entrance
242 TVector3 dpos;
243 if (fUseFixedErrors) {
244 dpos[0] = fDXsil; // 0.01 / TMath::Sqrt(12); //AZ
245 // dpos[1] = 0.1234; //AZ
246 dpos[1] = fDYsil; // 0.021; //AZ - as in cbmroot
247 } else {
248 bmnHit->PositionError(dpos);
249 }
250
251 Int_t sens = 1;
252 Int_t detId = kSILICON << 24 | (bmnHit->GetStation() + nVspStat + 1) << 16 | (bmnHit->GetModule() + 1) << 4
253 | sens << 1;
254 new ((*fCbmHitsArray)[fCbmHitsArray->GetEntriesFast()]) CbmStsHit(detId, pos, dpos, 0.0, 0, 0);
255 CbmStsHit* hit = (CbmStsHit*)fCbmHitsArray->At(fCbmHitsArray->GetEntriesFast() - 1);
256 if (fExp && fAlignmentEnabled)
257 ApplyAlignment(hit); // AZ-120223 - alignment
258 else if (!fExp) {
259 // AZ-140824 - shift coordinates
260 hit->SetY(hit->GetY() + 0.30);
261 if (fMCcorr) {
262 int iq = 0;
263 Double_t pmom = 0;
264 if (bmnHit->GetRefIndex() >= 0) {
265 FairMCPoint* mcp = (FairMCPoint*)fSilPoints->UncheckedAt(bmnHit->GetRefIndex());
266 CbmMCTrack* mctr = (CbmMCTrack*)fMCTracks->UncheckedAt(mcp->GetTrackID());
267 pmom = mctr->GetP();
268 if (TDatabasePDG::Instance()->GetParticle(mctr->GetPdgCode()))
269 iq = TMath::Nint(TDatabasePDG::Instance()->GetParticle(mctr->GetPdgCode())->Charge() / 3);
270 }
271 // Correct MC hit
272 CorrectHitMC(hit, pmom, iq);
273 if (iq == 0)
274 iq = 1; // for backward compatibility
275
276 Bool_t ok = CorrectEffic(hit, iq); // efficiency correction
277 if (!ok) {
278 // Remove hit
279 fCbmHitsArray->RemoveLast();
280 continue;
281 }
282 } // if (fMCcorr)
283 }
284 FairRootManager::Instance()->SetUseFairLinks(kTRUE);
285 hit->ResetLinks();
286 hit->SetLinks(bmnHit->GetLinks());
287 FairRootManager::Instance()->SetUseFairLinks(kFALSE);
288 hit->SetRefIndex(bmnHit->GetRefIndex());
289 if (fExp)
290 hit->SetRefIndex(iHit); // AZ-180226 - original hit index (for original coordinates)
291 hit->fDigiF =
292 bmnHit->GetLowerClusterIndex() + 2000000;
293 hit->fDigiB =
294 bmnHit->GetUpperClusterIndex() + 3000000;
295 // AZ-221022
296 StripCluster* upp = (StripCluster*)fBmnSilUpperClusters->UncheckedAt(hit->fDigiB % 1000000);
297 StripCluster* low = (StripCluster*)fBmnSilLowerClusters->UncheckedAt(hit->fDigiF % 1000000);
298 Double_t uppq = upp->TotalSignal;
299 Double_t lowq = low->TotalSignal;
300 Double_t corrq = (uppq - lowq) / (uppq + lowq);
301 hit->SetSignalDiv(corrq);
302 uppq = upp->Width; // N of strips
303 lowq = low->Width;
304 corrq = (uppq - lowq) / (uppq + lowq);
305 hit->SetTimeStamp(corrq); // just using available space
306 hit->SetStrips(low->MeanPosition, upp->MeanPosition); // AZ-190823
307 // hit->SetWidths (low->Width, upp->Width); //AZ-190823
308 } // for (Int_t iHit = 0; iHit < fBmnSilHitsArray->GetEntriesFast();
309 } // if (fBmnSilHitsArray)
310
311 if (fBmnVspHitsArray)
312 AddVspHits(); // AZ-060625
313
314 sw.Stop();
315 workTime += sw.RealTime();
316}
317// -------------------------------------------------------------------------
318
319// ----- Private method Init -------------------------------------------
320InitStatus BmnToCbmHitConverter::Init()
321{
322 // Get input array
323 FairRootManager* ioman = FairRootManager::Instance();
324 if (!ioman)
325 Fatal("Init", "No FairRootManager");
326 fBmnGemHitsArray = (TClonesArray*)ioman->GetObject(fBmnGemHitsBranchName);
327 fBmnSilHitsArray = (TClonesArray*)ioman->GetObject(fBmnSilHitsBranchName);
328 fBmnVspHitsArray = (TClonesArray*)ioman->GetObject(fBmnVspHitsBranchName);
329
330 fBmnGemLowerClusters = (TClonesArray*)ioman->GetObject("BmnGemLowerCluster");
331 fBmnGemUpperClusters = (TClonesArray*)ioman->GetObject("BmnGemUpperCluster");
332 fBmnSilLowerClusters = (TClonesArray*)ioman->GetObject("BmnSiliconLowerCluster");
333 fBmnSilUpperClusters = (TClonesArray*)ioman->GetObject("BmnSiliconUpperCluster");
334 fBmnVspLowerClusters = (TClonesArray*)ioman->GetObject("BmnVSPLowerCluster");
335 fBmnVspUpperClusters = (TClonesArray*)ioman->GetObject("BmnVSPUpperCluster");
336 if (!fBmnGemHitsArray && !fBmnSilHitsArray) {
337 cout << "BmnToCbmHitConverter::Init(): branches with inner hits not found! Task will be deactivated" << endl;
338 SetActive(kFALSE);
339 return kERROR;
340 }
341 fGemPoints = (TClonesArray*)ioman->GetObject("StsPoint"); // AZ-291223
342 fSilPoints = (TClonesArray*)ioman->GetObject("SiliconPoint"); // AZ-291223
343 fMCTracks = (TClonesArray*)ioman->GetObject("MCTrack"); // AZ-291223
344
345 // Register output array
346 fCbmHitsArray = new TClonesArray("CbmStsHit");
347 ioman->Register("StsHit", "STSHIT", fCbmHitsArray, kTRUE);
348
349 TString gPathConfig = gSystem->Getenv("VMCWORKDIR");
350
351 // Create GEM detector ------------------------------------------------------
352 GemStationSet = new BmnGemStripStationSet(gPathConfig + "/parameters/gem/XMLConfigs/" + fGemConfigFile);
353 SilStationSet = new BmnSiliconStationSet(gPathConfig + "/parameters/silicon/XMLConfigs/" + fSilConfigFile);
354 if (fVspConfigFile != "")
355 fVspStationSet = new BmnVSPStationSet(gPathConfig + "/parameters/vsp/XMLConfigs/" + fVspConfigFile);
356
357 FairField* mField = FairRunAna::Instance()->GetField();
358 if (mField == nullptr) {
359 LOG(error) << "BmnToCbmHitConverter::Init(): mField has null pointer! Task will be deactivated";
360 SetActive(kFALSE);
361 return kERROR;
362 }
363 fieldFlag = (TMath::Abs(mField->GetBy(0, 0, 135)) > 0.5) ? kTRUE : kFALSE;
364
365 return kSUCCESS;
366}
367// -------------------------------------------------------------------------
368
369void BmnToCbmHitConverter::OnlineWrite(const std::unique_ptr<TTree>& resultTree)
370{
371 if (!IsActive())
372 return;
373
374 resultTree->Branch("StsHit", &fCbmHitsArray);
375 resultTree->Fill();
376}
377
378void BmnToCbmHitConverter::SetFixedErrors(Float_t dXgem, Float_t dYgem, Float_t dXsil, Float_t dYsil)
379{
380 fUseFixedErrors = kTRUE;
381 fDXgem = dXgem;
382 fDYgem = dYgem;
383 fDXsil = dXsil;
384 fDYsil = dYsil;
385}
386
387void BmnToCbmHitConverter::Finish()
388{
389 printf("Work time of BmnToCbmHitConverter: %4.2f sec.\n", workTime);
390}
391
392// -------------------------------------------------------------------------
393
394void BmnToCbmHitConverter::ApplyAlignment(CbmStsHit* hit)
395{
396 // AZ-120223 - apply alignment
397
398 // AZ-070223 - alignment parameters from Igor
399 // /home/ruf/NICA/bmnroot230107
400 // /data4/ruf/Run8raw/7759/00/runAll.sh
401 // GEM: xhit-=StsDX[station][sector], yhit-=StsDY[station][sector], zhit+=GemDZ10[station-4][sector/2]
402 // Si: xhit-= SiDX[stn][sct]; yhit-= SiDY[stn][sct]; zhit+=SiDZ[stn][sct];
403
404 const Int_t MStn = 11 /*, NSct[MStn] = { 6,10,14,18, 8,8,8,8,8,8,8 }*/;
405 /*Double_t StsZ[MStn][18] = {
406 {19.3535, 20.5035, 19.2535, 16.2435, 15.0335, 16.4235},
407 {29.406, 28.231, 29.401, 28.086, 29.256, 23.807, 24.977, 23.827, 25.112, 23.892},
408 {36.9615, 38.1615, 36.8915, 38.1315, 36.8215, 37.9965, 36.8515, 33.534, 32.364, 33.519, 32.344, 33.709, 32.524, 33.719},
409 {46.788, 45.583, 46.733, 45.588, 46.773, 45.438, 46.613, 45.433, 46.618, 41.217, 42.342, 41.212, 42.367, 41.222, 42.512,
410 41.342, 42.467, 41.312}, {65.61, 65.61, 65.61, 65.61, 60.91, 60.91, 60.91, 60.91},
411 {92.392, 92.392, 92.392, 92.392, 97.092, 97.092, 97.092, 97.092},
412 {128.362, 128.362, 128.362, 128.362, 123.662, 123.662, 123.662, 123.662},
413 {154.842, 154.842, 154.842, 154.842, 159.542, 159.542, 159.542, 159.542},
414 {190.839, 190.839, 190.839, 190.839, 186.139, 186.139, 186.139, 186.139},
415 {217.437, 217.437, 217.437, 217.437, 222.137, 222.137, 222.137, 222.137},
416 {253.251, 253.251, 253.251, 253.251, 248.551, 248.551, 248.551, 248.551} };
417 */
418 double GemDZ10[7][4] = {{0.243654, 0.250775, 0.173749, 0.161366}, // for [station][Quater]
419 {-0.143876, -0.150758, 0.109959, 0.156433}, // same values from GEM-6 base
420 {-0.110145, -0.165981, -0.213889, -0.229543},
421 {0.103517, 0.036342, 0.086135, 0.0941369},
422 {-0.037648, -0.0941807, -0.180758, -0.176327},
423 {-0.0674491, -0.0367844, -0.0972473, -0.0584354},
424 {0, 0, 0, 0}};
425
426 // CbmStsHit: hit-track. Shifts of intermediate GEM stations with respect of first and last GEMs
427 Double_t StsDX[MStn][8] = {
428 {0., 0., 0., 0., 0., 0., 0., 0.},
429 {0., 0., 0., 0., 0., 0., 0., 0.},
430 {0., 0., 0., 0., 0., 0., 0., 0.},
431 {0., 0., 0., 0., 0., 0., 0., 0.},
432 {2.28613e-01, 2.34606e-01, 3.70449e-01, 4.39267e-01, 2.67664e-01, 2.84859e-01, 3.41005e-01, 4.16753e-01},
433 {3.17030e-01, 2.45618e-01, 4.29878e-01, 4.29169e-01, 2.60443e-01, 1.79925e-01, 3.76252e-01, 3.71674e-01},
434 {1.25543e-01, 1.37935e-01, 2.77096e-01, 3.51704e-01, 2.16012e-01, 2.12139e-01, 2.88912e-01, 3.57699e-01},
435 {3.17647e-01, 2.36295e-01, 3.75431e-01, 3.80575e-01, 1.68884e-01, 8.80524e-02, 2.51805e-01, 2.56513e-01},
436 {-3.01883e-02, -2.87880e-02, 1.09744e-01, 1.90429e-01, 1.78162e-02, 1.32671e-02, 9.23760e-02, 1.67779e-01},
437 {1.43966e-01, 5.53303e-02, 2.48555e-01, 2.49193e-01, 3.88544e-02, -3.31393e-02, 7.96085e-02, 7.41273e-02},
438 {0., -5.63645e-04, 0.05, 8.74027e-04, 0., 3.60763e-03, 0.074, 1.49797e-01}};
439 Double_t StsDY[MStn][8] = {
440 {0., 0., 0., 0., 0., 0., 0., 0.},
441 {0., 0., 0., 0., 0., 0., 0., 0.},
442 {0., 0., 0., 0., 0., 0., 0., 0.},
443 {0., 0., 0., 0., 0., 0., 0., 0.},
444 {-2.40830e-01, -2.27117e-01, -4.41393e-01, -2.98662e-01, -1.47376e-01, -5.64919e-01, 3.60250e-02, -2.15290e-01},
445 {-4.83443e-01, -3.63780e-01, -3.52015e-01, -3.36903e-01, -1.14008e-01, -3.49136e-01, -2.14579e-01,
446 -6.23524e-01},
447 {-1.91062e-01, -2.04515e-01, -3.79972e-01, -2.57092e-01, -1.01400e-01, -5.18137e-01, 5.39497e-02, -1.85679e-01},
448 {-4.55871e-01, -3.48696e-01, -2.33191e-01, -2.40044e-01, -5.37329e-02, -3.00596e-01, -1.98740e-01,
449 -6.04983e-01},
450 {-2.31062e-01, -2.58276e-01, -4.03759e-01, -2.70196e-01, -8.53213e-02, -5.14891e-01, 3.15495e-02, -2.08464e-01},
451 {-2.58821e-01, -1.58621e-01, -9.84743e-02, -1.05756e-01, 1.54258e-01, -1.04935e-01, 1.83280e-02, -3.86622e-01},
452 {0., -5.86937e-03, -0.23, -1.45342e-01, 0., -3.98117e-01, 0.15, -1.04839e-01}};
453
454 Double_t SiDX[4][18] = {{2.68835e-01, 2.39831e-01, 2.52506e-01, 2.05206e-01, 2.11052e-01, 0.20737190, 0, 0, 0, 0, 0,
455 0, 0, 0, 0, 0, 0, 0},
456 {3.15120e-01, 2.74671e-01, 2.55372e-01, 2.83446e-01, 2.79317e-01, 3.48623e-01, 3.18654e-01,
457 3.00181e-01, 3.16157e-01, 3.80363e-01, 0, 0, 0, 0, 0, 0, 0, 0},
458 {1.65037e-01, 2.71806e-01, 2.66914e-01, 2.40166e-01, 2.59614e-01, 2.66288e-01, 2.83516e-01,
459 2.88311e-01, 2.61280e-01, 3.02609e-01, 2.96469e-01, 2.98769e-01, 2.66650e-01, 3.72627e-01,
460 0, 0, 0, 0},
461 {3.01013e-01, 3.62381e-01, 3.03249e-01, 3.06502e-01, 2.75747e-01, 2.95530e-01, 3.05152e-01,
462 3.42706e-01, 3.09787e-01, 2.24795e-01, 2.31845e-01, 2.59297e-01, 2.17342e-01, 2.20424e-01,
463 2.06861e-01, 2.49563e-01, 2.64846e-01, 2.02358e-01}};
464 Double_t SiDY[4][18] = {
465 {-0.272688, -0.300649, -0.343639, -0.313878, -0.331129, -0.299529, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
466 {-0.342847, -0.343267, -0.328436, -0.354188, -0.332445, -0.349645, -0.332025, -0.294098, -0.249696, -0.217585,
467 0, 0, 0, 0, 0, 0, 0, 0},
468 {-0.398784, -0.342968, -0.325462, -0.325316, -0.338778, -0.338966, -0.324549, -0.339105, -0.314266, -0.301756,
469 -0.288295, -0.272864, -0.272913, -0.21294, 0, 0, 0, 0},
470 {-0.269827, -0.227419, -0.295412, -0.274608, -0.322856, -0.33417, -0.361013, -0.397192, -0.368905, -0.217688,
471 -0.240085, -0.266281, -0.285738, -0.269306, -0.292405, -0.297714, -0.323987, -0.379705}};
472 Double_t SiDZ[4][18] = {
473 {-0.4, -0.3, -0.23, -0.27, -0.3, -0.23, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
474 {-0.5, -0.4, -0.4, -0.28, -0.25, -0.45, -0.45, -0.3, -0.23, -0.2, 0, 0, 0, 0, 0, 0, 0, 0},
475 {-0.3, -0.5, -0.45, -0.4, -0.3, -0.2, -0.2, -0.25, -0.2, -0.3, -0.3, -0.4, -0.5, -0.27, 0, 0, 0, 0},
476 {-0.5, -0.6, -0.4, -0.5, -0.3, -0.35, -0.2, -0.1, -0.4, -0.6, -0.6, -0.6, -0.5, -0.35, -0.2, -0.1, 0.1, 0.}};
477
478 // AZ-090223 - Si module rotations around OZ (dx = p0 + p1*y)
479 const Double_t p0si[4][18] = {{3.865e-2, 6.497e-2, 3.848e-2, -3.535e-2, -2.711e-2, -4.438e-2}, // st1
480 {-3.011e-2, -6.724e-3, -2.784e-2, -1.761e-2, -1.431e-2, 2.650e-2, 1.818e-2, 1.052e-2,
481 1.510e-2, 1.977e-2}, // st2
482 {-1.639e-2, -1.628e-2, -4.100e-2, -2.524e-2, -2.082e-2, -3.096e-3, -1.535e-2,
483 -1.594e-2, 1.974e-2, 1.212e-2, 1.859e-2, 3.067e-2, 2.477e-2, -3.289e-3}, // st3
484 {5.774e-2, 5.280e-2, 4.679e-2, 6.866e-2, 6.183e-2, 4.256e-2, 4.247e-2, -5.477e-3,
485 2.597e-2, 1.403e-2, -2.016e-2, -3.445e-2, -3.345e-2, -3.685e-2, -4.381e-2, -6.655e-2,
486 -2.791e-2, -4.438e-2}};
487 const Double_t p1si[4][18] = {
488 {-8.887e-3, -1.128e-2, -6.696e-3, -1.042e-2, -7.374e-3, -1.303e-2}, // st1
489 {4.338e-3, 8.247e-4, 3.374e-3, 1.151e-3, 1.616e-3, 7.210e-3, 3.103e-3, 1.564e-3, 2.896e-3, 5.346e-3}, // st2
490 {3.678e-3, 1.777e-3, 4.869e-3, 2.708e-3, 2.610e-3, 4.911e-4, 2.753e-3, -2.473e-3, 3.427e-3, 1.598e-3, 2.388e-3,
491 3.975e-3, 3.602e-3, 1.143e-3}, // st3
492 {-4.152e-3, -6.035e-3, -4.377e-3, -6.320e-3, -5.058e-3, -3.510e-3, -3.763e-3, 4.285e-4, -5.000e-3, 1.801e-3,
493 -3.673e-3, -3.629e-3, -3.629e-3, -3.660e-3, -4.075e-3, -7.053e-3, -4.439e-3, -6.091e-3}};
494 // Second iteration
495 const Double_t p10si[4][18] = {{-2.385e-3, 6.372e-4, 8.679e-4, -2.758e-3, -5.495e-3, -4.838e-3}, // st1
496 {-3.930e-3, 2.765e-3, -3.730e-4, -5.525e-4, 1.433e-3, 1.042e-3, 2.216e-3, 3.619e-3,
497 1.800e-3, 1.076e-3}, // st2
498 {2.693e-3, 2.794e-3, -4.447e-3, -3.237e-4, 8.975e-4, 8.401e-4, -2.997e-3, -3.868e-3,
499 2.667e-3, -6.283e-4, 3.312e-3, 3.420e-3, 5.233e-3, -4.670e-3}, // st3
500 {7.950e-3, 2.818e-3, -1.178e-3, 3.522e-3, -1.457e-4, -2.013e-3, 3.139e-3, -1.060e-2,
501 -5.446e-3, 4.349e-3, 9.024e-4, -1.297e-3, -2.652e-3, -7.524e-3, -4.564e-3,
502 -9.118e-3, 8.307e-3, -5.750e-3}};
503 const Double_t p11si[4][18] = {
504 {3.525e-4, 4.803e-5, -4.898e-4, -5.709e-4, -1.111e-3, -1.295e-3}, // st1
505 {5.953e-4, -4.350e-4, 6.053e-5, 4.750e-5, 1.259e-4, 6.429e-4, 2.556e-4, 3.749e-4, 1.500e-4, 7.110e-4}, // st2
506 {1.719e-5, -3.435e-4, 6.088e-4, 1.112e-4, 2.683e-4, -6.623e-5, 4.340e-4, -6.263e-4, 4.717e-4, -1.268e-4,
507 4.721e-4, 5.582e-4, 6.571e-4, -4.956e-4}, // st3
508 {-4.333e-4, -3.287e-4, -7.432e-6, -4.067e-4, -6.684e-5, 3.667e-5, -4.441e-4, 8.977e-4, 7.495e-4, 5.340e-4,
509 -1.572e-4, -5.901e-5, -3.272e-4, -6.975e-4, -3.440e-4, -9.403e-4, 5.092e-4, -1.083e-3}};
510 // AZ-120223 - extra dZ-shifts of GEMs
511 /*Double_t dzAZ[9][4] = {{0.40,0.45,0.00,0.00}, {0.05,0.05,0.35,0.35}, {-0.15,-0.15,-0.30,-0.30},
512 {0.18,0.18,0.30,0.30}, {0.05,0.05,0.00,0.00}, {0.05,0.05,0.08,0.08},
513 {-0.15,-0.15,0.00,0.00}, {0}, {0}}; */
514
515 static TH2D *hLorCorX[2][17] = {{nullptr}, {nullptr}}, *hLorCorY[2][17] = {{nullptr}, {nullptr}};
516 static int irun = -1;
517 static Double_t zStat[19] = {0};
518
519 if (irun <= 0)
520 irun = ((DstEventHeader*)FairRootManager::Instance()->GetObject("DstEventHeader."))->GetRunId(); // AZ-010723
521
522 // Lorentz shifts
523 // AZ-010723 if (hLorCor[0][0] == nullptr) {
524 if (irun <= 0)
525 return;
526 else if (hLorCorX[0][0] == nullptr) {
527 // Station Z-positions
528 // AZ-160226 - VSP
529 int nsta0 = 0;
530 if (fVspStationSet)
531 nsta0 = fVspStationSet->GetNStations();
532 for (int ista = 0; ista < nsta0; ++ista)
533 zStat[ista] = fVspStationSet->GetZStationPosition(ista);
534 //
535 int nsta = SilStationSet->GetNStations();
536 for (int ista = 0; ista < nsta; ++ista)
537 zStat[ista + nsta0] = SilStationSet->GetZStationPosition(ista);
538 int nsta1 = GemStationSet->GetNStations();
539 for (int ista = 0; ista < nsta1; ++ista)
540 zStat[ista + nsta + nsta0] = GemStationSet->GetZStationPosition(ista);
541 if (fVerbose > 0)
542 for (int j = 0; j < nsta0 + nsta + nsta1; ++j)
543 cout << " Stat: " << j << " " << zStat[j] << endl;
544
545 map<int, vector<Double_t>> lorCorX[2], lorCorY[2];
546 int nx, ny;
547
548 for (int igem = 0; igem < 2; ++igem) {
549 ReadCorrections(irun, igem, lorCorX, lorCorY, nx, ny);
550 Double_t dx = 160.0 / nx, dy = 45.0 / ny;
551 // Double_t dx = 20.0, dy = 15.0;
552
553 for (int iupd = 0; iupd < 2; ++iupd) {
554 int istb = (igem == 0) ? 0 : nsta + nsta0;
555 int iste = (igem == 0) ? nsta + nsta0 : nsta + nsta1 + nsta0;
556 if (iste > (int)lorCorX[0].size())
557 --iste; // exclude last small GEM
558
559 for (int ist = istb; ist < iste; ++ist) {
560 TString hName = "hLorCorX_";
561 hName += iupd;
562 hName += "_";
563 hName += ist;
564 // hLorCor[iupd][ist] = new TH2D(hName, "", 8, -80, 80, 3, 0, 45);
565 int nxfsd = 0;
566 if (igem)
567 hLorCorX[iupd][ist] = new TH2D(hName, "", nx, -80, 80, ny, 0, 45);
568 else {
569 nxfsd = lorCorX[0][ist + 1].size();
570 hLorCorX[iupd][ist] = new TH2D(hName, "", nxfsd, 0, nxfsd, ny, 0, ny);
571 }
572
573 hName = "hLorCorY_";
574 hName += iupd;
575 hName += "_";
576 hName += ist;
577 if (igem)
578 hLorCorY[iupd][ist] = new TH2D(hName, "", nx, -80, 80, ny, 0, 45);
579 else
580 hLorCorY[iupd][ist] = new TH2D(hName, "", nxfsd, 0, nxfsd, ny, 0, ny);
581
582 if (igem) {
583 if (lorCorX[iupd].find(ist + 1) == lorCorX[iupd].end())
584 continue; // AZ-160226
585 // for (int j = 0; j < nx*ny; ++j) {
586 for (int j = 0; j < 24; ++j) {
587 // int ireg1 = j - 13;
588 // if (ireg1 < 0) ++ireg1;
589 int ireg1 = j - 12;
590 int ireg = TMath::Abs(ireg1);
591 if (ireg1 < 0)
592 --ireg;
593 int iy2 = ireg / 4;
594 int ix2 = ireg % 4;
595 Double_t x = (ix2 + 0.5) * dx;
596 if (ireg1 < 0)
597 x = -x;
598 Double_t y = (iy2 + 0.5) * dy;
599 hLorCorX[iupd][ist]->Fill(x, y, lorCorX[iupd][ist + 1][j]); // X-corrections
600 if (lorCorY[iupd].find(ist + 1) != lorCorY[iupd].end())
601 hLorCorY[iupd][ist]->Fill(x, y, lorCorY[iupd][ist + 1][j]); // Y-corrections
602 } // for (int j = 0;
603 } else {
604 for (int j = 0; j < nxfsd; ++j) {
605 hLorCorX[iupd][ist]->Fill(j + 0.5, 0.5, lorCorX[iupd][ist + 1][j]); // X-corrections
606 hLorCorY[iupd][ist]->Fill(j + 0.5, 0.5, lorCorY[iupd][ist + 1][j]); // Y-corrections
607 }
608 } // if (igem)
609 } // for (int ist = istb;
610 } // for (int iupd = 0;
611 } // for (int igem = 0;
612 } // else if (hLorCorX[0][0] == nullptr)
613
614 int ista = hit->GetStationNr() - 1;
615 int isec = hit->GetSectorNr() - 1;
616
617 if (hit->GetSystemId() != kGEM && irun < 9020) {
618 // if (0) {
619 // AZ-070223 - alignment for run 8
620 hit->SetX(hit->GetX() - SiDX[ista][isec]);
621 hit->SetY(hit->GetY() - SiDY[ista][isec]);
622 hit->SetZ(hit->GetZ() + SiDZ[ista][isec]);
623 // AZ-090223 - Si rotation around Z-axis (w/out changing Y-coordinate)
624 hit->SetX(hit->GetX()
625 - (p0si[ista][isec] / 2 + p1si[ista][isec] / 2 * hit->GetY() + p10si[ista][isec]
626 + p11si[ista][isec] * hit->GetY()));
627 }
628 // AZ-110223
629 if (hit->GetSystemId() == kGEM && irun < 9020) {
630 // if (0) {
631 // AZ-070223 - apply alignment (from Igor) for run 8
632 hit->SetX(hit->GetX() - StsDX[ista][isec]);
633 hit->SetY(hit->GetY() - StsDY[ista][isec]);
634 hit->SetZ(hit->GetZ() + GemDZ10[ista - 4][isec / 2]);
635
636 // hit->SetZ(hit->GetZ() + dzAZ[ista - 4][isec / 2]);
637 // 2 missed strips
638 if (ista == 4 && isec == 3) {
639 if (hit->GetX() >= -28.8 && hit->GetX() < -18.5) {
640 hit->SetX(hit->GetX() - 0.16);
641 hit->SetY(hit->GetY() - 0.16 / TMath::Tan(15.0 * TMath::DegToRad()));
642 }
643 }
644
645 // AZ-280223 Apply rotation around OY for lower GEMs
646 if (isec > 3) {
647 // Double_t xshift[9] = {0.32, 0.51, 0.63, 0.89, 1.02, 1.232, 1.368};
648 Double_t xshift[9] = {0.024, 0.058, 0.039, 0.050, 0.005, 0.028, 0.000};
649 Double_t alpha = 0.0055;
650 Double_t xnew = hit->GetX() * TMath::Cos(alpha) - (hit->GetZ() - zStat[ista]) * TMath::Sin(alpha);
651 Double_t znew = hit->GetX() * TMath::Sin(alpha) + (hit->GetZ() - zStat[ista]) * TMath::Cos(alpha);
652 hit->SetX(xnew + xshift[ista - 4]);
653 // if (ista == 6) znew -= 0.05;
654 // else if (ista == 7) znew -= 0.1;
655 // else if (ista == 8) znew -= 0.1;
656 // else if (ista == 10) znew += 0.2; //AZ-010323
657 hit->SetZ(znew + zStat[ista]);
658 } // else if (ista == 10) hit->SetZ (hit->GetZ() - 0.15); //AZ-010323
659
660 // AZ-030323 Lorentz shifts per station
661 /*
662 if (fieldFlag) {
663 Double_t lors[2][9] = {{ 0.017, -0.027, 0.021, -0.018, 0.014, 0.009, -0.038},
664 {-0.012, 0.023, -0.030, 0.039, -0.033, 0.013, -0.001}};
665 hit->SetX (hit->GetX() - lors[isec/4][ista-4] * 0.7);
666 }
667 */
668 } // if (hit->GetSystemId() == kGEM)
669
670 if (fieldFlag && hit->GetSystemId() == kGEM) {
671 // Lorentz corrections
672 int iupd = isec / 4;
673 Double_t xh = hit->GetX(), yh = TMath::Abs(hit->GetY()), xh1 = xh, yh1 = yh;
674 TH2D* h2 = hLorCorX[iupd][ista];
675 if (h2 == nullptr)
676 return;
677 if (xh1 >= h2->GetXaxis()->GetXmax())
678 xh1 = h2->GetXaxis()->GetXmax() - 0.1;
679 else if (xh1 < h2->GetXaxis()->GetXmin())
680 xh1 = h2->GetXaxis()->GetXmin() + 0.1;
681 if (yh1 >= h2->GetYaxis()->GetXmax())
682 yh1 = h2->GetYaxis()->GetXmax() - 0.1;
683 else if (yh1 < h2->GetYaxis()->GetXmin())
684 yh1 = h2->GetYaxis()->GetXmin() + 0.1;
685 // if (yh < h2->GetYaxis()->GetXmin() || yh >= h2->GetYaxis()->GetXmax()) cout << " yy " << isec << " " << yh <<
686 // endl;
687 Double_t dx = 0.0;
688 if (irun < 9020)
689 dx = h2->Interpolate(xh1, yh1); // AZ-170226 - for run8
690 else { // AZ-170226
691 int jx = h2->GetXaxis()->FindBin(xh1);
692 int jy = h2->GetYaxis()->FindBin(yh1);
693 dx = h2->GetBinContent(jx, jy);
694 }
695 hit->SetX(xh - dx * 0.65);
696 // hit->SetX (xh - dx * 0.85); // AZ-170226 - check
697
698 h2 = hLorCorY[iupd][ista];
699 Double_t dy = 0.0;
700 if (irun < 9020)
701 dy = h2->Interpolate(xh1, yh1); // AZ-170226 - for run8
702 else { // AZ-170226
703 int jx = h2->GetXaxis()->FindBin(xh1);
704 int jy = h2->GetYaxis()->FindBin(yh1);
705 dy = h2->GetBinContent(jx, jy);
706 }
707 if (iupd)
708 dy = -dy;
709 // cout << iupd << " " << ista << " " << dy << endl;
710 Double_t ynew = yh - dy * 0.65;
711 // Double_t ynew = yh - dy * 0.85; // AZ-170226 - check
712 if (iupd)
713 ynew = -ynew;
714 hit->SetY(ynew);
715 } else if (hit->GetSystemId() != kGEM) {
716 int iupd = isec / hLorCorX[0][ista]->GetNbinsX();
717 int ireg = isec % hLorCorX[0][ista]->GetNbinsX();
718 TH2D* h2 = hLorCorX[iupd][ista];
719 Double_t dx = h2->GetBinContent(ireg + 1, 1);
720 hit->SetX(hit->GetX() - dx * 0.65);
721 // hit->SetX(hit->GetX() - dx * 0.85); // AZ-170226 - check
722 h2 = hLorCorY[iupd][ista];
723 Double_t dy = h2->GetBinContent(ireg + 1, 1);
724 hit->SetY(hit->GetY() - dy * 0.65);
725 // hit->SetY(hit->GetY() - dy * 0.85); // AZ-170226 - check
726 }
727}
728
729// -------------------------------------------------------------------------
730
731void BmnToCbmHitConverter::ReadCorrections(int irun,
732 int igem,
733 std::map<int, std::vector<Double_t>>* lorCorX,
734 std::map<int, std::vector<Double_t>>* lorCorY,
735 int& nx,
736 int& ny)
737{
738 // Read file with "effective" Lorentz-corrections
739
740 ifstream fin;
741 TString fileName, chline;
742
743 fileName = gSystem->ExpandPathName("$VMCWORKDIR");
744 // AZ-030825 fileName += "/input/";
745 fileName += "/parameters/reco/"; // AZ-030825
746 if (igem)
747 fileName += "LorentzGemCorr.txt";
748 else
749 fileName += "LorentzFsdCorr.txt";
750 fin.open(fileName.Data());
751 if (!fin.is_open()) { // AZ-040825
752 cout << " !!! Lorentz correction file " << fileName << " not found !!!" << endl;
753 exit(1);
754 }
755
756 map<int, std::size_t> runsMap;
757 TObjArray* tokens = nullptr;
758
759 while (1) {
760 // while (!fin.eof()) {
761 size_t pos = fin.tellg();
762 chline.ReadLine(fin);
763 if (fin.eof())
764 break;
765 if (chline.Contains("Run")) {
766 // new run
767 // cout << chline << endl;
768 tokens = chline.Tokenize(" ");
769 TString run = ((TObjString*)tokens->UncheckedAt(1))->String();
770 run.Atoi();
771 // runsMap[run.Atoi()] = fin.tellg();
772 runsMap[run.Atoi()] = pos;
773 // cout << run.Atoi() << " " << pos << endl;
774 }
775 if (tokens) {
776 tokens->Delete();
777 delete tokens;
778 tokens = nullptr;
779 }
780 }
781 // if (chline[0] == '*') continue; // comment
782 // cout << runsMap.size() << endl;
783
784 auto iter = runsMap.upper_bound(irun);
785 if (iter != runsMap.begin())
786 --iter;
787 cout << " ---> Using corrections for run " << iter->first << " " << iter->second << " " << igem << endl;
788
789 fin.clear();
790 fin.seekg(iter->second);
791 chline.ReadLine(fin);
792 // cout << chline << endl;
793
794 // Read detectors
795 while (1) {
796 chline.ReadLine(fin);
797 if (fin.eof())
798 break;
799 if (chline.Contains("Station")) {
800 // cout << chline << endl;
801 tokens = chline.Tokenize(" ");
802 TString stat = ((TObjString*)tokens->UncheckedAt(1))->String();
803 int ista = stat.Atoi();
804 int iupd = 0; // half-station
805 if (chline.Contains("lower"))
806 iupd = 1;
807 nx = ((TObjString*)tokens->UncheckedAt(3))->String().Atoi();
808 ny = ((TObjString*)tokens->UncheckedAt(4))->String().Atoi();
809 int ixdata = 0;
810 if (chline.Contains("X-"))
811 ixdata = 1;
812 // cout << ista << " " << iupd << " " << nx << " " << ny << endl;
813 if (nx * ny < 1)
814 cout << " !!! Wrong input line !!! " << chline << " " << nx << " " << ny << endl;
815 if (tokens) {
816 tokens->Delete();
817 delete tokens;
818 tokens = nullptr;
819 }
820 chline.ReadLine(fin);
821 tokens = chline.Tokenize(", ");
822 int nreg = tokens->GetEntriesFast();
823 if (nreg != nx * ny) {
824 cout << " !!! Wrong input line !!! " << chline << " " << nreg << " " << nx * ny << endl;
825 exit(0);
826 }
827 auto lorCor = vector<Double_t>(nreg);
828 // cout << tokens->GetEntriesFast() << " " << ((TObjString*)tokens->UncheckedAt(0))->String() << endl;
829 for (int itok = 0; itok < nreg; ++itok) {
830 Double_t corr = ((TObjString*)tokens->UncheckedAt(itok))->String().Atof();
831 // lorCor[iupd][ista][itok] = corr;
832 lorCor[itok] = corr;
833 }
834 if (ixdata)
835 lorCorX[iupd][ista] = lorCor;
836 else
837 lorCorY[iupd][ista] = lorCor;
838 if (tokens) {
839 tokens->Delete();
840 delete tokens;
841 tokens = nullptr;
842 }
843 } else if (chline.Contains("*")) {
844 continue; // comment
845 } else
846 break;
847 } // while (1)
848 fin.close();
849}
850
851// -------------------------------------------------------------------------
852
853Bool_t BmnToCbmHitConverter::CorrectEffic(CbmStsHit* hit, int iq)
854{
855 // AZ-211223 Apply efficiency corrections for Monte Carlo
856
857 /*
858 if (hit->GetStationNr() == 3) {
859 // Apply efficiency for 1 station
860 if (gRandom->Rndm() > 0.90) return kFALSE; // 90%
861 }
862 */
863 // return kTRUE; // no inefficiency applied
864
865 static int first = 1;
866 const int nstat = 19;
867 static TH2D* hEff[2][nstat] = {{nullptr}, {nullptr}};
868
869 if (first) {
870 first = 0;
871 // Get histograms with station efficiencies
872 TString fileName;
873 fileName = gSystem->ExpandPathName("$VMCWORKDIR");
874 // AZ-220725 fileName += "/input/";
875 fileName += "/parameters/reco/"; // AZ-220725
876 TString fileEff = "SiGemEffic.root";
877 if (FairRun::Instance()->GetTask("STS Vector Finder V9")) {
878 fileEff = "SiGemEfficV8.root";
879 }
880 fileName += fileEff;
881 TFile* f = TFile::Open(fileName);
882 cout << " *** BmnToCbmHitConverter::CorrectEffic(): Reading Monte Carlo efficiency correction file " << endl;
883 cout << " *** " << fileName << endl;
884
885 for (int ipm = 0; ipm < 2; ++ipm) {
886
887 for (int i = 0; i < nstat; ++i) {
888 TString hName = (ipm == 0) ? "hEffpos_" : "hEffneg_";
889 hName += (1 + i);
890 TObject* obj = f->FindObjectAny(hName);
891 if (obj == nullptr)
892 break;
893 hEff[ipm][i] = (TH2D*)obj;
894 // cout << hName << endl;
895 }
896 }
897 }
898
899 int ipm = (iq > 0) ? 0 : 1;
900 int ista = hit->GetStationNr() - 1;
901 if (hEff[ipm][ista] == nullptr)
902 return kTRUE; // last small GEM
903 int ix = hEff[ipm][ista]->GetXaxis()->FindBin(hit->GetX());
904 int iy = hEff[ipm][ista]->GetYaxis()->FindBin(hit->GetY());
905 if (ix < 1 || ix > hEff[ipm][ista]->GetNbinsX() || iy < 1 || iy > hEff[ipm][ista]->GetNbinsY())
906 return kTRUE;
907 Double_t eff = hEff[ipm][ista]->GetBinContent(ix, iy);
908
909 // eff = 99; //AZ-170824 - no eff. correction
910
911 if (eff >= 0.999999)
912 return kTRUE;
913 if (gRandom->Rndm() <= eff)
914 return kTRUE;
915 return kFALSE;
916}
917
918// -------------------------------------------------------------------------
919
920void BmnToCbmHitConverter::CorrectHitMC(CbmStsHit* hit, Double_t pmom0, int iq)
921{
922 // AZ-270824 - apply coordinate corrections for MC hits
923 // AZ-040125 - use 2-component Gaussian description
924
925 Float_t dx = 0;
926
927 // AZ-310823 Extra shift for MC
928 if (hit->GetSystemId() == kGEM) {
929 dx = 0.0100;
930 if (hit->GetStationNr() % 2) {
931 if (hit->GetSectorNr() < 5)
932 dx = -dx;
933 } else if (hit->GetSectorNr() >= 5)
934 dx = -dx;
935 hit->SetX(hit->GetX() - dx * 0.65);
936 }
937
938 // return; // AZ-040125 - no corrections
939
940 static int first = 1;
941 // static map<int,Double_t> shiftsPos[2][19], sigmasPos[2][19], shiftsNeg[2][19], sigmasNeg[2][19];
942 static map<int, vector<Float_t>> paramsPos[2][19], paramsNeg[2][19];
943
944 // Read file with corrections
945 if (first) {
946 first = 0;
947 ifstream fin;
948 TString fileName, chline;
949
950 fileName = gSystem->ExpandPathName("$VMCWORKDIR");
951 // AZ-030825 fileName += "/input/SiGemHitCorrMC.txt";
952 fileName += "/parameters/reco/SiGemHitCorrMC.txt"; // AZ-030825
953 fin.open(fileName.Data());
954 if (!fin.is_open()) { // AZ-030825
955 cout << " !!! MC correction file " << fileName << " not found !!!" << endl;
956 exit(1);
957 }
958
959 TObjArray* tokens = nullptr;
960
961 while (1) {
962 chline.ReadLine(fin);
963 if (fin.eof())
964 break;
965 // cout << chline << endl;
966 if (chline.Contains("#"))
967 continue; // comment
968 tokens = chline.Tokenize(" ");
969 TString stations = ((TObjString*)tokens->UncheckedAt(0))->String();
970 int stat1 = stations.Atoi(), stat2 = stat1;
971 if (stations.Contains("-")) {
972 // Several stations
973 stat1 = stations.Atoi();
974 int ip = stations.First("-");
975 TString sstat = stations(ip + 1, stations.Length() - ip - 1);
976 stat2 = sstat.Atoi();
977 }
978 int pmin = (((TObjString*)tokens->UncheckedAt(1))->String()).Atof(); // in MeV
979 int pmax = (((TObjString*)tokens->UncheckedAt(2))->String()).Atoi(); // in MeV
980 Float_t ampl1 = (((TObjString*)tokens->UncheckedAt(3))->String()).Atof();
981 Float_t shift1 = (((TObjString*)tokens->UncheckedAt(4))->String()).Atof() * 1.1;
982 Float_t sigma1 = (((TObjString*)tokens->UncheckedAt(5))->String()).Atof() * 1.1;
983 Float_t ampl2 = (((TObjString*)tokens->UncheckedAt(6))->String()).Atof();
984 Float_t shift2 = (((TObjString*)tokens->UncheckedAt(7))->String()).Atof() * 1.1;
985 Float_t sigma2 = (((TObjString*)tokens->UncheckedAt(8))->String()).Atof() * 1.1;
986 int ixy = 0;
987 if ((((TObjString*)tokens->UncheckedAt(9))->String()).Contains("Y"))
988 ixy = 1;
989 Float_t ampl = ampl1 + ampl2;
990 // cout << stat1 << " " << stat2 << " " << pmax << " " << ampl1 << " " << shift1 << " " << sigma1 << " " <<
991 // ixy
992 // << endl;
993
994 for (int ista = stat1; ista <= stat2; ++ista) {
995 if (pmin < 0) {
996 paramsNeg[ixy][ista][pmax] = {ampl1 / ampl, shift1, sigma1, ampl2 / ampl, shift2, sigma2};
997 } else {
998 paramsPos[ixy][ista][pmax] = {ampl1 / ampl, shift1, sigma1, ampl2 / ampl, shift2, sigma2};
999 }
1000 }
1001 if (tokens) {
1002 tokens->Delete();
1003 delete tokens;
1004 tokens = nullptr;
1005 }
1006 }
1007 }
1008
1009 int ista = hit->GetStationNr(), pmom = 0;
1010 Float_t sigx = 0.03, sigy = 0.05, dy = 0; // for backward compatibility
1011 if (hit->GetSystemId() != kGEM) {
1012 sigx = 0.008;
1013 sigy = 0.04;
1014 }
1015 if (iq != 0) {
1016 int pmomev = int(pmom0 * 1000 / iq);
1017 auto par = (iq > 0) ? &paramsPos : &paramsNeg;
1018 auto params = *par;
1019 if (params[0][ista].size() != 0) {
1020 pmom = TMath::Min(pmomev, params[0][ista].rbegin()->first - 1);
1021 auto iter = params[0][ista].lower_bound(pmom);
1022 vector<Float_t>& pars = iter->second;
1023 dx = pars[1];
1024 sigx = pars[2];
1025 if (gRandom->Rndm() > pars[0]) {
1026 dx = pars[4];
1027 sigx = pars[5];
1028 }
1029 }
1030 if (params[1][ista].size() != 0) {
1031 pmom = TMath::Min(pmomev, params[1][ista].rbegin()->first - 1);
1032 auto iter = params[1][ista].lower_bound(pmom);
1033 vector<Float_t>& pars = iter->second;
1034 dy = pars[1];
1035 sigy = pars[2];
1036 if (gRandom->Rndm() > pars[0]) {
1037 dy = pars[4];
1038 sigy = pars[5];
1039 }
1040 }
1041 }
1042
1043 hit->SetX(hit->GetX() + gRandom->Gaus(dx, sigx));
1044 hit->SetY(hit->GetY() + gRandom->Gaus(dy, sigy));
1045
1046 // cout << " zzz " << iq << " " << ista << " " << pmom0 << " " << pmom << " " << dx << " " << sigx << " " << dy << "
1047 // " << sigy << endl;
1048}
1049
1050// AZ-060625 -------------------------------------------------------------------------
1051
1052void BmnToCbmHitConverter::AddVspHits()
1053{
1054 // Add VSP (Vertex Silicon Plane - CBM sensors) hits
1055
1056 int nhits = fBmnVspHitsArray->GetEntriesFast();
1057
1058 for (Int_t iHit = 0; iHit < nhits; ++iHit) {
1059 BmnVSPHit* bmnHit = (BmnVSPHit*)fBmnVspHitsArray->UncheckedAt(iHit);
1060 TVector3 pos;
1061 bmnHit->Position(pos);
1062 pos[2] -= 0.0150; // AZ - shift to the entrance
1063 TVector3 dpos;
1064 if (fUseFixedErrors) {
1065 dpos[0] = fDXsil; // 0.01 / TMath::Sqrt(12); //AZ
1066 // dpos[1] = 0.1234; //AZ
1067 dpos[1] = fDYsil; // 0.021; //AZ - as in cbmroot
1068 } else {
1069 bmnHit->PositionError(dpos);
1070 }
1071
1072 Int_t sens = 1;
1073 Int_t detId = kVSP << 24 | (bmnHit->GetStation() + 1) << 16 | (bmnHit->GetModule() + 1) << 4 | sens << 1;
1074 new ((*fCbmHitsArray)[fCbmHitsArray->GetEntriesFast()]) CbmStsHit(detId, pos, dpos, 0.0, 0, 0);
1075 CbmStsHit* hit = (CbmStsHit*)fCbmHitsArray->At(fCbmHitsArray->GetEntriesFast() - 1);
1076 if (fExp && fAlignmentEnabled)
1077 // if (0)
1078 ApplyAlignment(hit);
1079 else if (!fExp) {
1080 /*
1081 int iq = 0;
1082 Double_t pmom = 0;
1083 if (bmnHit->GetRefIndex() >= 0) {
1084 FairMCPoint* mcp = (FairMCPoint*)fSilPoints->UncheckedAt(bmnHit->GetRefIndex());
1085 CbmMCTrack* mctr = (CbmMCTrack*)fMCTracks->UncheckedAt(mcp->GetTrackID());
1086 pmom = mctr->GetP();
1087 if (TDatabasePDG::Instance()->GetParticle(mctr->GetPdgCode()))
1088 iq = TMath::Nint(TDatabasePDG::Instance()->GetParticle(mctr->GetPdgCode())->Charge() / 3);
1089 }
1090 // Correct MC hit
1091 CorrectHitMC(hit, pmom, iq);
1092 if (iq == 0)
1093 iq = 1; // for backward compatibility
1094
1095 Bool_t ok = CorrectEffic(hit, iq); // efficiency correction
1096 if (!ok) {
1097 // Remove hit
1098 fCbmHitsArray->RemoveLast();
1099 continue;
1100 } */
1101 // AZ-140824 - shift coordinates
1102 hit->SetY(hit->GetY() + 0.30);
1103 }
1104 FairRootManager::Instance()->SetUseFairLinks(kTRUE);
1105 hit->ResetLinks();
1106 hit->SetLinks(bmnHit->GetLinks());
1107 FairRootManager::Instance()->SetUseFairLinks(kFALSE);
1108 hit->SetRefIndex(bmnHit->GetRefIndex());
1109 hit->fDigiF =
1110 bmnHit->GetLowerClusterIndex() + 4000000;
1111 hit->fDigiB =
1112 bmnHit->GetUpperClusterIndex() + 5000000;
1113 // AZ-221022
1114 StripCluster* upp = (StripCluster*)fBmnVspUpperClusters->UncheckedAt(hit->fDigiB % 1000000);
1115 StripCluster* low = (StripCluster*)fBmnVspLowerClusters->UncheckedAt(hit->fDigiF % 1000000);
1116 Double_t uppq = (upp) ? upp->TotalSignal : 0.0;
1117 Double_t lowq = (low) ? low->TotalSignal : 0.0;
1118 Double_t corrq = (uppq - lowq) / (uppq + lowq);
1119 hit->SetSignalDiv(corrq);
1120 uppq = (upp) ? upp->Width : 0.0; // N of strips
1121 lowq = (low) ? low->Width : 0.0;
1122 corrq = (uppq - lowq) / (uppq + lowq);
1123 hit->SetTimeStamp(corrq); // just using available space
1124 hit->SetStrips((low) ? low->MeanPosition : 0.0,
1125 (upp) ? upp->MeanPosition : 0.0); // AZ-190823
1126 // hit->SetWidths (low->Width, upp->Width); //AZ-190823
1127 } // for (Int_t iHit = 0; iHit < nhits
1128}
i
int i
Definition P4_F32vec4.h:22
f
float f
Definition P4_F32vec4.h:21
kSILICON
@ kSILICON
Definition BmnDetectorList.h:25
kGEM
@ kGEM
Definition BmnDetectorList.h:15
kVSP
@ kVSP
Definition BmnDetectorList.h:21
BmnGemStripLayer::IsPointInsideStripLayer
Bool_t IsPointInsideStripLayer(Double_t x, Double_t y)
Definition BmnGemStripLayer.cxx:254
BmnGemStripModule::GetStripLayer
BmnGemStripLayer & GetStripLayer(Int_t num)
Definition BmnGemStripModule.h:71
BmnGemStripStationSet::GetStation
BmnGemStripStation * GetStation(Int_t station_num)
Definition BmnGemStripStationSet.h:45
BmnGemStripStationSet::GetZStationPosition
Double_t GetZStationPosition(Int_t station_num)
Definition BmnGemStripStationSet.cxx:132
BmnGemStripStation::GetModule
BmnGemStripModule * GetModule(Int_t module_num)
Definition BmnGemStripStation.cxx:183
BmnHit::GetModule
Int_t GetModule()
Definition BmnHit.h:77
BmnHit::IsPseudo
Bool_t IsPseudo()
Definition BmnHit.h:145
BmnHit::GetUpperClusterIndex
Int_t GetUpperClusterIndex()
Definition BmnHit.h:125
BmnHit::GetLowerClusterIndex
Int_t GetLowerClusterIndex()
Definition BmnHit.h:129
BmnHit::GetStation
Short_t GetStation() const
Definition BmnHit.h:45
BmnSiliconStationSet::GetZStationPosition
Double_t GetZStationPosition(Int_t station_num)
Definition BmnSiliconStationSet.cxx:80
BmnTask
BmnTask.
Definition BmnTask.h:13
BmnToCbmHitConverter::Init
virtual InitStatus Init()
Definition BmnToCbmHitConverter.cxx:320
BmnToCbmHitConverter::SetFixedErrors
void SetFixedErrors(Float_t dXgem=0.015, Float_t dYgem=0.058, Float_t dXsil=0.003, Float_t dYsil=0.021)
Definition BmnToCbmHitConverter.cxx:378
BmnToCbmHitConverter::ApplyAlignment
void ApplyAlignment(CbmStsHit *hit)
Definition BmnToCbmHitConverter.cxx:394
BmnToCbmHitConverter::CorrectEffic
Bool_t CorrectEffic(CbmStsHit *hit, int iq)
Definition BmnToCbmHitConverter.cxx:853
BmnToCbmHitConverter::OnlineWrite
virtual void OnlineWrite(const std::unique_ptr< TTree > &resultTree)
Write task resul to tree.
Definition BmnToCbmHitConverter.cxx:369
BmnToCbmHitConverter::CorrectHitMC
void CorrectHitMC(CbmStsHit *hit, Double_t pmom, int iq)
Definition BmnToCbmHitConverter.cxx:920
BmnToCbmHitConverter::Exec
virtual void Exec(Option_t *opt)
Definition BmnToCbmHitConverter.cxx:87
BmnToCbmHitConverter::ReadCorrections
void ReadCorrections(int irun, int igem, std::map< int, std::vector< Double_t > > *locCorX, std::map< int, std::vector< Double_t > > *locCorY, int &nx, int &ny)
Definition BmnToCbmHitConverter.cxx:731
BmnVSPStationSet::GetZStationPosition
Double_t GetZStationPosition(Int_t station_num)
Definition BmnVSPStationSet.cxx:79
CbmMCTrack::GetP
Double_t GetP() const
Definition CbmMCTrack.h:68
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition CbmMCTrack.h:56
CbmStsHit::SetStrips
void SetStrips(Double_t s1, Double_t s2)
Definition CbmStsHit.h:87
CbmStsHit::GetStationNr
virtual Int_t GetStationNr() const
Definition CbmStsHit.h:66
CbmStsHit::SetSignalDiv
void SetSignalDiv(Double_t sigDiv)
Definition CbmStsHit.h:86
CbmStsHit::fDigiB
Int_t fDigiB
Definition CbmStsHit.h:93
CbmStsHit::GetSectorNr
Int_t GetSectorNr() const
Definition CbmStsHit.h:68
CbmStsHit::GetSystemId
Int_t GetSystemId() const
Definition CbmStsHit.h:64
CbmStsHit::fDigiF
Int_t fDigiF
Definition CbmStsHit.h:92
StripCluster::TotalSignal
Double_t TotalSignal
Definition BmnStripData.h:181
StripCluster::MeanPosition
Double_t MeanPosition
Definition BmnStripData.h:180
run
-clang-format