BmnRoot
Loading...
Searching...
No Matches
BmnLambdaEmbedding.cxx
Go to the documentation of this file.
1#include "BmnLambdaEmbedding.h"
2#include "CbmStsPoint.h"
3#include "UniRun.h"
4
5// Defines GCC_DIAGNOSTIC_AWARE if GCC 4.7 or above
6#define GCC_DIAGNOSTIC_AWARE 1
7#if GCC_DIAGNOSTIC_AWARE
8# pragma GCC diagnostic ignored "-Wunknown-pragmas"
9#endif
10
11BmnLambdaEmbedding::BmnLambdaEmbedding()
12: fInfo(nullptr),
13 fSim(nullptr),
14 fReco(nullptr),
15 fLambdaSim(nullptr),
16 fMCTracks(nullptr),
17 fVertices(nullptr),
18 fGemPoints(nullptr),
19 fGemDigits(nullptr),
20 fGemMatch(nullptr),
21 fCscPoints(nullptr),
22 fCscDigits(nullptr),
23 fCscMatch(nullptr),
24 fSiliconPoints(nullptr),
25 fSiliconDigits(nullptr),
26 fSiliconMatch(nullptr),
27 fADC32(nullptr),
28 fADC128(nullptr),
29 fSync(nullptr),
30 fLambdaStore(nullptr),
31 fHeader(nullptr),
32 fMon(nullptr)
33{
34 // Useful tools
35 fInfo = new BmnLambdaMisc(); // Initialize useful tools to work with mapping ...
36}
37
38BmnLambdaEmbedding::BmnLambdaEmbedding(TString raw, TString sim, TString reco, TString out, Int_t nEvs, TString lambdaStore)
39: fInfo(nullptr),
40 fSim(nullptr),
41 fReco(nullptr),
42 fLambdaSim(nullptr),
43 fMCTracks(nullptr),
44 fVertices(nullptr),
45 fGemPoints(nullptr),
46 fGemDigits(nullptr),
47 fGemMatch(nullptr),
48 fCscPoints(nullptr),
49 fCscDigits(nullptr),
50 fCscMatch(nullptr),
51 fSiliconPoints(nullptr),
52 fSiliconDigits(nullptr),
53 fSiliconMatch(nullptr),
54 fADC32(nullptr),
55 fADC128(nullptr),
56 fSync(nullptr),
57 fLambdaStore(nullptr),
58 fWrittenStores(nullptr),
59 fHeader(nullptr),
60 fMon(nullptr),
61 normUtils(nullptr),
62 fPdgDecParticle(3122),
63 fPdgPostive(2212),
64 fPdgNegative(-211)
65{
66 fDataFileName = raw;
67 fDigiFileName = out;
68 fStorePath = "outputDataStore";
69
70 // Initialize steering flags by default ...
71 doLambdaStore = kTRUE;
72 doListOfEventsWithReconstructedVertex = kTRUE;
73 doSimulateLambdaThroughSetup = kTRUE;
74 doRawRootConvertion = kTRUE;
75 doEmbedding = kTRUE;
76 doDecode = kTRUE;
77 doPrintStoreInfo = kFALSE;
78 doEmbeddingMonitor = kTRUE;
79
80 isGemEmbedded = kTRUE;
81 isSilEmbedded = kTRUE;
82 isCscEmbedded = kFALSE;
83
84 fEvents = nEvs;
85
86 // Open simu file
87 fSim = new TChain("bmndata");
88 fSim->Add(sim.Data());
89 fSim->SetBranchAddress("MCTrack", &fMCTracks);
90
91 // Open dst file
92 fReco = new TChain("bmndata");
93 fReco->Add(reco.Data());
94 fReco->SetBranchAddress("BmnVertex", &fVertices);
95 fReco->SetBranchAddress("DstEventHeader.", &fHeader);
96
97 // Useful tools
98 fInfo = new BmnLambdaMisc(); // Initialize useful tools to work with mapping ...
99
100 fNstores = 10;
101 fNeventsToBeSimulated = 50;
102 fStoreToProcess = -1;
103
104 if (!lambdaStore.IsNull()) {
105 fWrittenStores = new TClonesArray*[fNstores];
106 for (Int_t iStore = 0; iStore < fNstores; iStore++)
107 fWrittenStores[iStore] = nullptr;
108
109 for (Int_t iStore = 0; iStore < fNstores; iStore++)
110 fWrittenStores[iStore] = new TClonesArray("BmnParticleStore");
111
112 TChain* ch = new TChain("bmndata");
113 ch->Add(lambdaStore.Data());
114
115 TClonesArray* tmp = nullptr;
116 ch->SetBranchAddress("BmnParticleStore", &tmp);
117
118 for (Int_t iStore = 0; iStore < fNstores; iStore++) {
119 ch->GetEntry(iStore);
120
121 fWrittenStores[iStore] = (TClonesArray*) (tmp->UncheckedAt(0))->Clone();
122 }
123
124 delete ch;
125 } else
126 fLambdaStore = new TClonesArray("BmnParticleStore");
127
128 fRunId = -1;
129 fFieldScale = -1.;
130
131 // Default cut values ...
132 fZmin = -3.;
133 fZmax = +3.;
134 fEtaMin = 0.;
135 fEtaMax = 100.;
136 fPhiMin = 0.;
137 fPhiMax = 360.;
138 fMomMin = 0.;
139 fNHitsProton = 4;
140 fNHitsPion = 4;
141 isUseRealSignal = kTRUE;
142
143 fSignal.push_back(1000);
144 fSignal.push_back(1500);
145 fSignal.push_back(1000);
146
147 // Resetting default values of flags depending on an user's scenario if recognized
148 // Scenario 1 (just selecting Lambdas ...)
149 if ((!out.Contains("digi") && !out.IsNull()) && raw.IsNull() && !sim.IsNull() && reco.IsNull() && lambdaStore.IsNull()) {
150 cout << "Recognized scenario 1" << endl;
151 DoLambdaStore(kTRUE);
152 DoPrintStoreInfo(kTRUE);
153 SetUseRealSignals(kFALSE);
154 DoListOfEventsWithReconstructedVertex(kFALSE);
155 DoSimulateLambdaThroughSetup(kFALSE);
156 DoRawRootConvertion(kFALSE);
157 DoEmbeddingMonitor(kFALSE);
158 DoEmbedding(kFALSE);
159 DoDecode(kFALSE);
160 } // Scenario 2 (passing all stores (or selected one) in file to simulation ...)
161 else if (raw.IsNull() && sim.IsNull() && !reco.IsNull() && out.IsNull() && !lambdaStore.IsNull()) {
162 cout << "Recognized scenario 2 (or 3)" << endl;
163 DoLambdaStore(kFALSE);
164 DoPrintStoreInfo(kFALSE);
165 SetUseRealSignals(kFALSE);
166 DoSimulateLambdaThroughSetup(kTRUE);
167 DoRawRootConvertion(kFALSE);
168 DoEmbeddingMonitor(kFALSE);
169 DoEmbedding(kFALSE);
170 DoDecode(kFALSE);
171 }
172 // Scenario 4 (namely, embedding )
173 else if (!raw.IsNull() && (out.Contains("digi") && !out.IsNull()) && sim.IsNull() && !reco.IsNull() && lambdaStore.IsNull()) {
174 cout << "Recognized scenario 4" << endl;
175 DoLambdaStore(kFALSE);
176 DoPrintStoreInfo(kFALSE);
177 DoRawRootConvertion(kTRUE);
178 DoSimulateLambdaThroughSetup(kFALSE);
179 SetUseRealSignals(kTRUE);
180 DoEmbeddingMonitor(kTRUE);
181 DoEmbedding(kTRUE);
182 DoDecode(kTRUE);
183 }
184}
185
186void BmnLambdaEmbedding::Embedding() {
187 fFileNamePart = (fPdgDecParticle == 3122) ? "lambda" : (fPdgDecParticle == 310) ? "kaon" : "unknown";
188
189 // 0. Get a function from current dst to be used for rescaling ...
190 if (isUseRealSignal) {
191 TString digiFile = "";
192
193 TString eos_host = "root://ncm.jinr.ru/";
194 TString digiDir = "/eos/nica/bmn/exp/digi/";
195 TString run = "run7";
196 TString release = "20.02.0";
197 TString dirList[2] = {"4720-5186_BMN_Krypton", "3590-4707_BMN_Argon"};
198
199 fReco->GetEntry();
200
201 Int_t idx = (UniRun::GetRun(7, fHeader->GetRunId())->GetBeamParticle() == "Ar") ? 1 :
202 (UniRun::GetRun(7, fHeader->GetRunId())->GetBeamParticle() == "Kr") ? 0 : -1;
203
204 TString command = TString::Format("xrdcp %s/%s/%s/%s/%s/bmn_run%d_digi.root %s",
205 eos_host.Data(), digiDir.Data(), run.Data(), release.Data(), dirList[idx].Data(), fHeader->GetRunId(), gSystem->GetFromPipe("pwd").Data());
206
207 gSystem->Exec(command.Data());
208
209 TFile* file = new TFile(Form("bmn_run%d_digi.root", fHeader->GetRunId()));
210
211 if (file->IsOpen())
212 digiFile = file->GetName();
213
214 delete file;
215
216 if (!digiFile.IsNull())
217 normUtils = new SignalNormalizationUtils(digiFile, TString(fReco->GetFile()->GetName()));
218 }
219
220 // 1. Create a store with lambdas ...
221 if (doLambdaStore) {
222 if (!fDigiFileName.Contains("digi") && !fDigiFileName.IsNull()) {
223 TClonesArray* store = CreateLambdaStore();
224
225 TFile* f = new TFile(fDigiFileName.Data(), "recreate");
226 TTree* tree = new TTree("bmndata", "bmndata");
227
228 TClonesArray* tmp = new TClonesArray("BmnParticleStore");
229 tree->Branch("BmnParticleStore", &tmp);
230
231 for (Int_t iEntry = 0; iEntry < store->GetEntriesFast(); iEntry++) {
232 tmp->Delete();
233 TClonesArray* arr = (TClonesArray*) store->UncheckedAt(iEntry);
234
235 new ((*tmp)[tmp->GetEntriesFast()]) BmnParticleStore(*((BmnParticleStore*) arr));
236 tree->Fill();
237 }
238
239 tree->Write();
240
241 delete tmp;
242 delete tree;
243 delete f;
244
245 } else
246 CreateLambdaStore();
247 }
248
249 if (doPrintStoreInfo)
250 PrintStoreInfo();
251
252 // 2. Create list of eventId for reconstructed events where the primary vertex is assumed to be defined ...
253 map <UInt_t, TVector3> EventIdsVpMap;
254 if (doListOfEventsWithReconstructedVertex)
255 EventIdsVpMap = ListOfEventsWithReconstructedVp();
256
257 // 3. Get information on mag. field in the run ...
258 if (doSimulateLambdaThroughSetup) {
259 if (fRunId != -1) {
260 Double_t map_current = 55.87;
261 UniRun* runInfo = UniRun::GetRun(7, fRunId);
262 if (runInfo)
263 fFieldScale = *runInfo->GetFieldVoltage() / map_current;
264 }
265 }
266
267 // 4. Loop over store with lambdas ...
268 //Int_t nThreads = 1;
269#if defined(_OPENMP)
270 nThreads = fNstores;
271 omp_set_num_threads(nThreads);
272#endif
273
274 if (doSimulateLambdaThroughSetup) {
275 Int_t low = (fStoreToProcess != -1) ? fStoreToProcess : 0;
276 Int_t up = (fStoreToProcess != -1) ? (fStoreToProcess + 1) : fNstores;
277
278#pragma omp parallel for
279 for (Int_t iLambda = low; iLambda < up; iLambda++) {
280 BmnParticleStore* lambda = (BmnParticleStore*) fWrittenStores[iLambda];
281 Double_t eta = lambda->GetEta();
282 Double_t phi = lambda->GetPhi();
283 Double_t p = lambda->GetP();
284 Double_t pdg = lambda->GetPdg();
285
286 Int_t iVertex = 0;
287
288 for (auto it : EventIdsVpMap) {
289 // Get reconstructed primary vertex ...
290 TVector3 Vp = it.second;
291
292 SimulateLambdaPassing(pdg, p, TVector2(eta, phi), Vp, iLambda, iVertex);
293
294 iVertex++;
295 }
296 }
297 }
298
299 // 5. Find at least one lambda to be reconstructed for a given Vp
300 map <UInt_t, pair <TVector3, TVector3>> GoodLambdaForEachVertex; // <evId ---> (<iLambda, iVertex, iEvent>, <Vp>)
301
302 map <UInt_t, BmnParticlePair> embMonitorMap; // evId --> particlePair (see physics/particles/BmnParticlePair.h(*.cxx))
303
304 Int_t iVertex = 0;
305 for (auto it : EventIdsVpMap) {
306 UInt_t id = it.first;
307 TVector3 lambda(-1, -1, -1); // <iLambda, iVertex, iEvent>, no appropriate any lambda found yet for current eventId ...
308 TVector3 Vp = it.second;
309
310 GoodLambdaForEachVertex[id] = make_pair(lambda, Vp); // Put info for a given eventId
311
312 BmnParticlePair aPairToBeWritten;
313 embMonitorMap[id] = aPairToBeWritten;
314
315 for (Int_t iLambda = 0; iLambda < fNstores; iLambda++) {
316 BmnParticlePair currentPair;
317
318 Int_t eve = FindReconstructableLambdaFromStore(iLambda, iVertex, currentPair);
319
320 // Put current pair into the embMonitorMap ...
321 auto itEmbInfoMap = embMonitorMap.find(id);
322 if (itEmbInfoMap != embMonitorMap.end())
323 (*itEmbInfoMap).second = currentPair;
324
325 if (eve != -1) {
326 TVector3 tmp = TVector3(iLambda, iVertex, eve);
327
328 // Trying to find id to be replaced ...
329 auto itr = GoodLambdaForEachVertex.find(id);
330 if (itr != GoodLambdaForEachVertex.end())
331 (*itr).second.first = tmp;
332
333 break;
334 }
335 }
336 iVertex++;
337 }
338
339 // 6. Monitor events with embedded lambdas ...
340 if (doEmbeddingMonitor) {
341 fMon = new BmnLambdaEmbeddingMonitor();
342 TFile* file = new TFile(Form("embedding_%d.root", fRunId), "recreate");
343 TTree* tree = new TTree("bmndata", "bmndata");
344
345 tree->Branch("EmbeddingMonitor.", &fMon);
346 for (auto it : GoodLambdaForEachVertex) {
347 UInt_t id = it.first;
348
349 TVector3 par = it.second.first;
350
351 fMon->SetId(id);
352 fMon->SetStoreVertexEvent(par);
353
354 auto itEmbInfoMap = embMonitorMap.find(id);
355 // if (itEmbInfoMap != embMonitorMap.end()) {
356 BmnParticlePair pair = (*itEmbInfoMap).second;
357
358 fMon->SetProtonP(pair.GetMomPart1());
359 fMon->SetPionP(pair.GetMomPart2());
360
361 fMon->SetTxProton(pair.GetTxPart1());
362 fMon->SetTxPion(pair.GetTxPart2());
363
364 fMon->SetTyProton(pair.GetTyPart1());
365 fMon->SetTyPion(pair.GetTyPart2());
366
367 fMon->SetNHitsProton(pair.GetNHitsPart1());
368 fMon->SetNHitsPion(pair.GetNHitsPart2());
369 // }
370
371 if ((Int_t) par[0] != -1 && (Int_t) par[1] != -1 && (Int_t) par[2] != -1)
372 fMon->IsEmbedded(kTRUE);
373
374 else
375 fMon->IsEmbedded(kFALSE);
376
377 tree->Fill();
378 }
379
380 tree->Write();
381 file->Close();
382 }
383
384 // 7. Create digi-arrays corresponding to a certain eventId ... <evId --> lambdaDigs>
385 map <UInt_t, vector < BmnStripDigit>> digsFromLambdasGem; // GEM
386 map <UInt_t, vector < BmnStripDigit>> digsFromLambdasSilicon; // SILICON
387 map <UInt_t, vector < BmnStripDigit>> digsFromLambdasCsc; // CSC
388
389 if (doEmbedding)
390 for (auto it : GoodLambdaForEachVertex) {
391 TVector3 par = it.second.first;
392 Int_t lambda = par[0];
393 Int_t vertex = par[1];
394 Int_t event = par[2];
395
396 if (lambda == -1 || vertex == -1 || event == -1)
397 continue;
398
399 if (isGemEmbedded)
400 digsFromLambdasGem[it.first] = GetDigitsFromLambda(TString::Format("%s/%s%d_vertex%d.root", fStorePath.Data(), fFileNamePart.Data(), lambda, vertex), event, "GEM");
401 if (isSilEmbedded)
402 digsFromLambdasSilicon[it.first] = GetDigitsFromLambda(TString::Format("%s/%s%d_vertex%d.root", fStorePath.Data(), fFileNamePart.Data(), lambda, vertex), event, "SILICON");
403 if (isCscEmbedded)
404 digsFromLambdasCsc[it.first] = GetDigitsFromLambda(TString::Format("%s/%s%d_vertex%d.root", fStorePath.Data(), fFileNamePart.Data(), lambda, vertex), event, "CSC");
405 }
406
407 cout << "GEM# " << digsFromLambdasGem.size() << " SILICON# " << digsFromLambdasSilicon.size() << " CSC# " << digsFromLambdasCsc.size() << endl;
408
409 // 7a. Collecting all MC-digits together ...
410 if (isUseRealSignal) {
411 TTree* tree = new TTree("bmndata", "bmndata");
412
413 TClonesArray* digitsGem = new TClonesArray("BmnGemStripDigit");
414 tree->Branch("BmnGemStripDigit", &digitsGem);
415
416 TClonesArray* digitsSil = new TClonesArray("BmnSiliconDigit");
417 tree->Branch("BmnSiliconDigit", &digitsSil);
418
419 for (auto digs : digsFromLambdasGem) {
420 digitsGem->Delete();
421
422 vector <BmnStripDigit> digitsInEvent = digs.second;
423
424 for (BmnStripDigit digit : digitsInEvent) {
425 BmnStripDigit* dig = &digit;
426
427 new ((*digitsGem) [digitsGem->GetEntriesFast()]) BmnStripDigit(dig);
428 }
429 tree->Fill();
430 }
431
432 for (auto digs : digsFromLambdasSilicon) {
433 digitsSil->Delete();
434
435 vector <BmnStripDigit> digitsInEvent = digs.second;
436
437 for (BmnStripDigit digit : digitsInEvent) {
438 BmnStripDigit* dig = &digit;
439
440 new ((*digitsSil) [digitsSil->GetEntriesFast()]) BmnStripDigit(dig);
441 }
442 tree->Fill();
443 }
444
445 normUtils->SetMcDataSet((TChain*) tree);
446 delete tree;
447 delete digitsGem;
448 delete digitsSil;
449
450 // Recalculating all signals according to the scal function got ...
451 TF1* scaleFuncGem = normUtils->GetRescaleFunction("GEM");
452 TF1* scaleFuncSil = normUtils->GetRescaleFunction("SILICON");
453
454 if (!scaleFuncGem || !scaleFuncSil) {
455 cout << "Something is wrong, exiting ..." << endl;
456 throw;
457 }
458
459 for (auto& digs : digsFromLambdasGem) {
460 vector <BmnStripDigit> &digitsInEvent = digs.second;
461
462 for (BmnStripDigit& digit : digitsInEvent)
463 digit.SetStripSignal(scaleFuncGem->Eval(digit.GetStripSignal()));
464 }
465
466 for (auto& digs : digsFromLambdasSilicon) {
467 vector <BmnStripDigit> &digitsInEvent = digs.second;
468
469 for (BmnStripDigit& digit : digitsInEvent)
470 digit.SetStripSignal(scaleFuncSil->Eval(digit.GetStripSignal()));
471 }
472 }
473
474 // 8. Make correspondence between evId and lambda digits with info on channel and serial ...
475 map <UInt_t, map < pair <Int_t, Int_t>, Long_t>> evIdGemChannelSerial; // GEM --- <evId ---> <digi index + ch, serial>>
476 map <UInt_t, map < pair <Int_t, Int_t>, Long_t>> evIdCscChannelSerial; // CSC --- <evId ---> <digi index + ch, serial>>
477 map <UInt_t, map < vector <Int_t>, Long_t>> evIdSiliconChannelSerial; // SILICON --- <evId ---> <digi index + ch + sample, serial>>
478
479 if (doEmbedding) {
480 // GEM
481 if (isGemEmbedded)
482 for (auto it : digsFromLambdasGem) {
483 vector <BmnStripDigit> digits = it.second;
484 map <pair <Int_t, Int_t>, Long_t> tmpMap = GetGemChannelSerialFromDigi(digits);
485 evIdGemChannelSerial[it.first] = tmpMap;
486 }
487
488 // SILICON
489 if (isSilEmbedded)
490 for (auto it : digsFromLambdasSilicon) {
491 vector <BmnStripDigit> digits = it.second;
492 map <vector <Int_t>, Long_t> tmpMap = GetSiliconChannelSerialFromDigi(digits);
493 evIdSiliconChannelSerial[it.first] = tmpMap;
494 }
495
496 // CSC
497 if (isCscEmbedded)
498 for (auto it : digsFromLambdasCsc) {
499 vector <BmnStripDigit> digits = it.second;
500 map <pair <Int_t, Int_t>, Long_t> tmpMap = GetCscChannelSerialFromDigi(digits);
501 evIdCscChannelSerial[it.first] = tmpMap;
502 }
503 }
504
505 // 9. Do *.data --> *raw.root convertion for requested events ...
506 DoRawRootFromBinaryData(fDataFileName);
507
508 // 10. Loop over *raw.root to embed lambda digits ...
509 TString rawRoot = "";
510 if (doEmbedding)
511 rawRoot = AddInfoToRawFile(digsFromLambdasGem, evIdGemChannelSerial,
512 digsFromLambdasSilicon, evIdSiliconChannelSerial,
513 digsFromLambdasCsc, evIdCscChannelSerial);
514
515 // 11. Start decoding ...
516 if (!rawRoot.IsNull() && doDecode)
517 StartDecodingWithEmbeddedLambdas(rawRoot);
518}
519
520Int_t BmnLambdaEmbedding::FindReconstructableLambdaFromStore(Int_t iLambda, Int_t iVertex, BmnParticlePair& pairFromLambda) {
521 TString fileName = TString::Format("%s/%s%d_vertex%d.root", fStorePath.Data(), fFileNamePart.Data(), iLambda, iVertex);
522 TChain ch("bmndata");
523 if (ch.Add(fileName.Data()) == 0)
524 return -1;
525
526 // Open lambda file
527 TClonesArray* tracks = nullptr;
528 TClonesArray* gemPoints = nullptr;
529 TClonesArray* siliconPoints = nullptr;
530 TClonesArray* cscPoints = nullptr;
531
532 ch.SetBranchAddress("MCTrack", &tracks);
533 ch.SetBranchAddress("StsPoint", &gemPoints);
534 ch.SetBranchAddress("SiliconPoint", &siliconPoints);
535 ch.SetBranchAddress("CSCPoint", &cscPoints);
536
537 for (Int_t iEv = 0; iEv < ch.GetEntries(); iEv++) {
538 ch.GetEntry(iEv);
539
540 // Loop over protons ...
541 // A proton we are looking for is single !
542 Int_t nHitsPerProton = 0;
543 Double_t Pprot = 0., TxProt = -1., TyProt = -1.;
544
545 for (Int_t iTrack = 0; iTrack < tracks->GetEntriesFast(); iTrack++) {
546 CbmMCTrack* track = (CbmMCTrack*) tracks->UncheckedAt(iTrack);
547
548 // Skip primary tracks
549 if (track->GetMotherId() == -1)
550 continue;
551
552 // Skipping particles not been protons
553 if (track->GetPdgCode() != fPdgPostive)
554 continue;
555
556 Int_t id = track->GetMotherId();
557
558 // Get secondary protons and pions from lambda decay
559 if (((CbmMCTrack*) tracks->UncheckedAt(id))->GetPdgCode() != fPdgDecParticle)
560 continue;
561
562 if (track->GetP() < 0.5)
563 continue;
564
565 const Int_t nStats = 10; // 6 GEMs + 3 SILICONs + 1 CSC [<0 ... 5>, <6, 7, 8>, <9>]
566 Int_t pointsOnGemsSiliconsAndCsc[nStats];
567 for (Int_t iStat = 0; iStat < nStats; iStat++)
568 pointsOnGemsSiliconsAndCsc[iStat] = 0;
569
570 // Loop over gem points ...
571 if (isGemEmbedded)
572 for (Int_t iGem = 0; iGem < gemPoints->GetEntriesFast(); iGem++) {
573 CbmStsPoint* point = (CbmStsPoint*) gemPoints->UncheckedAt(iGem);
574 if (point->GetTrackID() != iTrack)
575 continue;
576 pointsOnGemsSiliconsAndCsc[point->GetStation()]++;
577 }
578
579 // Loop over silicon points ...
580 if (isSilEmbedded)
581 for (Int_t iSilicon = 0; iSilicon < siliconPoints->GetEntriesFast(); iSilicon++) {
582 FairMCPoint* point = (FairMCPoint*) siliconPoints->UncheckedAt(iSilicon);
583 if (point->GetTrackID() != iTrack)
584 continue;
585 pointsOnGemsSiliconsAndCsc[DefineSiliconStatByZpoint(point->GetZ())]++;
586 }
587
588 // Loop over csc points ...
589 if (isCscEmbedded)
590 for (Int_t iCsc = 0; iCsc < cscPoints->GetEntriesFast(); iCsc++) {
591 BmnCSCPoint* point = (BmnCSCPoint*) cscPoints->UncheckedAt(iCsc);
592 if (point->GetTrackID() != iTrack)
593 continue;
594 pointsOnGemsSiliconsAndCsc[9]++;
595 }
596
597 // Here we decide on a track whether being accepted or not according a list of detectors to be used when embedding MC-data ...
598 Bool_t isPrelimAccepted[3] = {kTRUE, kTRUE, kTRUE}; // GEM, SILICON, CSC ...
599 Int_t pointCounter;
600
601 // Checking GEMs ...
602 if (isGemEmbedded) {
603 pointCounter = 0;
604 for (Int_t iStat = 0; iStat < 6; iStat++) {
605 if (pointsOnGemsSiliconsAndCsc[iStat] != 0)
606 pointCounter++;
607 }
608 if (pointCounter == 0)
609 isPrelimAccepted[0] = kFALSE;
610 }
611
612 // Checking SILICONs ...
613 if (isSilEmbedded) {
614 pointCounter = 0;
615 for (Int_t iStat = 6; iStat < 9; iStat++) {
616 if (pointsOnGemsSiliconsAndCsc[iStat] != 0)
617 pointCounter++;
618 }
619 if (pointCounter == 0)
620 isPrelimAccepted[1] = kFALSE;
621 }
622
623 // Checking CSC ...
624 if (isCscEmbedded) {
625 pointCounter = 0;
626 for (Int_t iStat = 9; iStat < nStats; iStat++) {
627 if (pointsOnGemsSiliconsAndCsc[iStat] != 0)
628 pointCounter++;
629 }
630 if (pointCounter == 0)
631 isPrelimAccepted[2] = kFALSE;
632 }
633
634 Bool_t doWeAcceptEvent = (isPrelimAccepted[0] && isPrelimAccepted[1] && isPrelimAccepted[2]) ? kTRUE : kFALSE;
635
636 if (!doWeAcceptEvent)
637 continue;
638
639 // Loop over statArray to avoid multiple Z ...
640 Bool_t isMultiplePointsOnZ = kFALSE;
641 for (Int_t iStat = 0; iStat < nStats; iStat++) {
642 if (pointsOnGemsSiliconsAndCsc[iStat] == 2) {
643 isMultiplePointsOnZ = kTRUE;
644 break;
645 } else
646 nHitsPerProton += pointsOnGemsSiliconsAndCsc[iStat];
647 }
648 if (isMultiplePointsOnZ || nHitsPerProton < 4 || nHitsPerProton > nStats)
649 continue;
650
651 Pprot = track->GetP();
652 TxProt = track->GetPx() / track->GetPz();
653 TyProt = track->GetPy() / track->GetPz();
654 }
655
656 // Skipping the next loop if not satisfied condition ...
657 if (nHitsPerProton < fNHitsProton || Abs(Pprot) < FLT_EPSILON)
658 continue;
659
660 // Loop over pions ...
661 Int_t nHitsPerPion = 0;
662 Double_t Ppion = 0., TxPion = -1., TyPion = -1.;
663
664 for (Int_t jTrack = 0; jTrack < tracks->GetEntriesFast(); jTrack++) {
665 CbmMCTrack* track = (CbmMCTrack*) tracks->UncheckedAt(jTrack);
666
667 // Skip primary tracks
668 if (track->GetMotherId() == -1)
669 continue;
670
671 // Skipping particles not been pions ...
672 if (track->GetPdgCode() != fPdgNegative)
673 continue;
674
675 Int_t id = track->GetMotherId();
676
677 // Get secondary protons and pions from lambda decay
678 if (((CbmMCTrack*) tracks->UncheckedAt(id))->GetPdgCode() != fPdgDecParticle)
679 continue;
680
681 if (track->GetP() < 0.5)
682 continue;
683
684 const Int_t nStats = 9; // 6 GEMs + 3 SILICONs [<0 ... 5>, <6, 7, 8>]
685 Int_t pointsOnGemsAndSilicon[nStats];
686 for (Int_t iStat = 0; iStat < nStats; iStat++)
687 pointsOnGemsAndSilicon[iStat] = 0;
688
689 // Loop over gem points ...
690 if (isGemEmbedded)
691 for (Int_t iGem = 0; iGem < gemPoints->GetEntriesFast(); iGem++) {
692 CbmStsPoint* point = (CbmStsPoint*) gemPoints->UncheckedAt(iGem);
693 if (point->GetTrackID() != jTrack)
694 continue;
695 pointsOnGemsAndSilicon[point->GetStation()]++;
696 }
697
698 // Loop over silicon points ...
699 if (isSilEmbedded)
700 for (Int_t iSilicon = 0; iSilicon < siliconPoints->GetEntriesFast(); iSilicon++) {
701 FairMCPoint* point = (FairMCPoint*) siliconPoints->UncheckedAt(iSilicon);
702 if (point->GetTrackID() != jTrack)
703 continue;
704 pointsOnGemsAndSilicon[DefineSiliconStatByZpoint(point->GetZ())]++;
705 }
706
707 // Here we decide on a track whether being accepted or not according a list of detectors to be used when embedding MC-data ...
708 Bool_t isPrelimAccepted[2] = {kTRUE, kTRUE}; // GEM, SILICON
709 Int_t pointCounter;
710
711 // Checking GEMs ...
712 if (isGemEmbedded) {
713 pointCounter = 0;
714 for (Int_t iStat = 0; iStat < 6; iStat++) {
715 if (pointsOnGemsAndSilicon[iStat] != 0)
716 pointCounter++;
717 }
718 if (pointCounter == 0)
719 isPrelimAccepted[0] = kFALSE;
720 }
721
722 // Checking SILICONs ...
723 if (isSilEmbedded) {
724 pointCounter = 0;
725 for (Int_t iStat = 6; iStat < nStats; iStat++) {
726 if (pointsOnGemsAndSilicon[iStat] != 0)
727 pointCounter++;
728 }
729 if (pointCounter == 0)
730 isPrelimAccepted[1] = kFALSE;
731 }
732
733 Bool_t doWeAcceptEvent = (isPrelimAccepted[0] && isPrelimAccepted[1]) ? kTRUE : kFALSE;
734
735 if (!doWeAcceptEvent)
736 continue;
737
738 // Loop over statArray to avoid multiple Z ...
739 Bool_t isMultiplePointsOnZ = kFALSE;
740 for (Int_t iStat = 0; iStat < nStats; iStat++) {
741 if (pointsOnGemsAndSilicon[iStat] == 2) {
742 isMultiplePointsOnZ = kTRUE;
743 break;
744 } else
745 nHitsPerPion += pointsOnGemsAndSilicon[iStat];
746 }
747 if (isMultiplePointsOnZ || nHitsPerPion < 4 || nHitsPerPion > nStats)
748 continue;
749
750 Ppion = track->GetP();
751 TxPion = track->GetPx() / track->GetPz();
752 TyPion = track->GetPy() / track->GetPz();
753 }
754
755 if (nHitsPerPion < fNHitsPion || Abs(Ppion) < FLT_EPSILON)
756 continue;
757
758 pairFromLambda.SetMomPair(Pprot, Ppion);
759 pairFromLambda.SetTxPair(TxProt, TxPion);
760 pairFromLambda.SetTyPair(TyProt, TyPion);
761 pairFromLambda.SetNHitsPair(nHitsPerProton, nHitsPerPion);
762
763 return iEv;
764 }
765 return -1;
766}
767
768void BmnLambdaEmbedding::DoRawRootFromBinaryData(TString raw) {
769 BmnRawDataDecoder* decoder = new BmnRawDataDecoder(raw, "", fEvents);
770 decoder->SetPeriodId(7);
771 decoder->SetRunId(fRunId);
772
773 if (doRawRootConvertion)
774 decoder->ConvertRawToRoot();
775 fRawRootFileName = decoder->GetRootFileName();
776
777 delete decoder;
778}
779
780void BmnLambdaEmbedding::StartDecodingWithEmbeddedLambdas(TString raw) {
781 BmnRawDataDecoder* rdd = new BmnRawDataDecoder("", "", fEvents);
782 BmnDecoder * decoder = rdd->GetDecoder();
783 decoder->SetPeriodId(7);
784 decoder->SetRunId(fRunId);
785 decoder->SetRawRootFile(raw);
786 decoder->SetDigiRootFile(fDigiFileName);
787 decoder->SetBmnSetup(kBMNSETUP);
788
789 std::map<DetectorId, bool> setup; // flags to determine BM@N setup
790 setup.insert(std::make_pair(kBC, 0)); // TRIGGERS
791 setup.insert(std::make_pair(kSILICON, 1)); // SILICON
792 setup.insert(std::make_pair(kGEM, 1)); // GEM
793
794 decoder->SetDetectorSetup(setup);
795
796 decoder->SetMSCMapping("MSC_map_Run7.txt");
797 decoder->SetTrigPlaceMapping("Trig_PlaceMap_Run7.txt");
798 decoder->SetTrigChannelMapping("Trig_map_Run7.txt");
799 decoder->SetSiliconMapping("SILICON_map_run7.txt");
800 decoder->SetGemMapping("GEM_map_run7.txt");
801 decoder->SetCSCMapping("CSC_map_period7.txt");
802
803 decoder->DecodeDataToDigi();
804
805 delete decoder;
806}
807
808map <UInt_t, TVector3> BmnLambdaEmbedding::ListOfEventsWithReconstructedVp() {
809 map <UInt_t, TVector3> EventIdsVpMap;
810 for (Int_t iEntry = 0; iEntry < fEvents; iEntry++) {
811 fReco->GetEntry(iEntry);
812
813 if (iEntry == 0)
814 fRunId = fHeader->GetRunId();
815
816 CbmVertex* vertex = (CbmVertex*) fVertices->UncheckedAt(0);
817
818 // Used cuts on primary vertex
819 if (vertex->GetZ() < fZmin || vertex->GetZ() > fZmax)
820 continue;
821
822 EventIdsVpMap[fHeader->GetEventId()] = TVector3(vertex->GetX(), vertex->GetY(), vertex->GetZ());
823 }
824 return EventIdsVpMap;
825}
826
827vector <BmnStripDigit> BmnLambdaEmbedding::GetDigitsFromLambda(TString lambdaEve, Int_t evNum, TString det) {
828 Bool_t isSilicon = (det.Contains("SILICON")) ? kTRUE : kFALSE;
829 Bool_t isGem = (det.Contains("GEM")) ? kTRUE : kFALSE;
830 Bool_t isCsc = (det.Contains("CSC")) ? kTRUE : kFALSE;
831
832 vector <BmnStripDigit> digits;
833
834 // Open file with simulated lambdas
835 fLambdaSim = new TChain("bmndata");
836 fLambdaSim->Add(lambdaEve.Data());
837 fLambdaSim->SetBranchAddress("MCTrack", &fMCTracks);
838
839 fLambdaSim->SetBranchAddress("SiliconPoint", &fSiliconPoints);
840 fLambdaSim->SetBranchAddress("BmnSiliconDigit", &fSiliconDigits);
841 fLambdaSim->SetBranchAddress("BmnSiliconDigitMatch", &fSiliconMatch);
842
843 fLambdaSim->SetBranchAddress("StsPoint", &fGemPoints);
844 fLambdaSim->SetBranchAddress("BmnGemStripDigit", &fGemDigits);
845 fLambdaSim->SetBranchAddress("BmnGemStripDigitMatch", &fGemMatch);
846
847 fLambdaSim->SetBranchAddress("CSCPoint", &fCscPoints);
848 fLambdaSim->SetBranchAddress("BmnCSCDigit", &fCscDigits);
849 fLambdaSim->SetBranchAddress("BmnCSCDigitMatch", &fCscMatch);
850
851 fLambdaSim->GetEntry(evNum);
852
853 Int_t idxPi = -1, idxProt = -1;
854
855 // Looking for protons and pions from lambda decay
856 for (Int_t iTrack = 0; iTrack < fMCTracks->GetEntriesFast(); iTrack++) {
857 CbmMCTrack* track = (CbmMCTrack*) fMCTracks->UncheckedAt(iTrack);
858
859 // Skip primary tracks
860 if (track->GetMotherId() == -1)
861 continue;
862
863 // Skip particles not been protons or pions
864 if (track->GetPdgCode() != fPdgPostive && track->GetPdgCode() != fPdgNegative)
865 continue;
866
867 Int_t pdg = track->GetPdgCode();
868 Int_t id = track->GetMotherId();
869
870 // Get secondary protons and pions from lambda decay
871 if (((CbmMCTrack*) fMCTracks->UncheckedAt(id))->GetPdgCode() != fPdgDecParticle)
872 continue;
873
874 if (pdg > 0)
875 idxProt = iTrack;
876 else
877 idxPi = iTrack;
878 }
879
880 // Loop over GEM-digits
881 if (isGem)
882 for (Int_t iDig = 0; iDig < fGemDigits->GetEntriesFast(); iDig++) {
883 BmnMatch* digiMatch = (BmnMatch*) fGemMatch->UncheckedAt(iDig);
884
885 Int_t idxPoint = digiMatch->GetMatchedLink().GetIndex();
886 FairMCPoint* point = (FairMCPoint*) fGemPoints->UncheckedAt(idxPoint);
887 if (point->GetTrackID() != idxProt && point->GetTrackID() != idxPi)
888 continue;
889
890 BmnStripDigit* digi = (BmnStripDigit*) fGemDigits->UncheckedAt(iDig);
891 // digi->SetName((point->GetTrackID() == idxProt) ? "PROTON" : "PION");
892 digits.push_back(*digi);
893 }
894
895 // Loop over CSC-digits
896 if (isCsc)
897 for (Int_t iDig = 0; iDig < fCscDigits->GetEntriesFast(); iDig++) {
898 BmnMatch* digiMatch = (BmnMatch*) fCscMatch->UncheckedAt(iDig);
899
900 Int_t idxPoint = digiMatch->GetMatchedLink().GetIndex();
901 FairMCPoint* point = (FairMCPoint*) fCscPoints->UncheckedAt(idxPoint);
902 if (point->GetTrackID() != idxProt && point->GetTrackID() != idxPi)
903 continue;
904
905 BmnStripDigit* digi = (BmnStripDigit*) fCscDigits->UncheckedAt(iDig);
906 digits.push_back(*digi);
907 }
908
909 // Loop over SILICON-digits
910 if (isSilicon)
911 for (Int_t iDig = 0; iDig < fSiliconDigits->GetEntriesFast(); iDig++) {
912 BmnMatch* digiMatch = (BmnMatch*) fSiliconMatch->UncheckedAt(iDig);
913
914 Int_t idxPoint = digiMatch->GetMatchedLink().GetIndex();
915 FairMCPoint* point = (FairMCPoint*) fSiliconPoints->UncheckedAt(idxPoint);
916 if (point->GetTrackID() != idxProt && point->GetTrackID() != idxPi)
917 continue;
918
919 BmnStripDigit* digi = (BmnStripDigit*) fSiliconDigits->UncheckedAt(iDig);
920 digits.push_back(*digi);
921 }
922
923 delete fLambdaSim;
924
925 return digits;
926}
927
928map <pair <Int_t, Int_t>, Long_t> BmnLambdaEmbedding::GetGemChannelSerialFromDigi(vector <BmnStripDigit> digits) {
929 map <Int_t, Int_t> digiToChannel; // (digi index in vector ---> channel)
930 map <pair <Int_t, Int_t>, Long_t> digiChannelToSerial; // (digi index in vector + channel ---> serial)
931
932 const Int_t nDigs = digits.size();
933
934 // Fill digiToChannel map ...
935 Long_t serial[nDigs];
936
937 for (Int_t iDigi = 0; iDigi < nDigs; iDigi++) {
938 serial[iDigi] = 0x0;
939 BmnStripDigit* dig = &digits[iDigi];
940 digiToChannel[iDigi] = fInfo->GemDigiToChannel(dig, serial[iDigi]);
941 }
942
943 // Fill digiChannelToSerial map ...
944 for (auto it : digiToChannel) {
945
946 // Skip digits with no channel found ...
947 if (it.second == -1)
948 continue;
949
950 Long_t value = (serial[it.first] == 0) ? fInfo->GemDigiChannelToSerial(make_pair(digits[it.first], it.second)) : serial[it.first];
951 digiChannelToSerial[make_pair(it.first, it.second)] = value;
952 }
953
954 return digiChannelToSerial;
955}
956
957map <pair <Int_t, Int_t>, Long_t> BmnLambdaEmbedding::GetCscChannelSerialFromDigi(vector <BmnStripDigit> digits) {
958 map <pair <Int_t, Int_t>, Long_t> digiChannelToSerial; // (digi index in vector + channel ---> serial)
959
960 for (size_t iDigi = 0; iDigi < digits.size(); iDigi++) {
961 BmnStripDigit* dig = &digits[iDigi];
962
963 Long_t serial = 0x0;
964
965 digiChannelToSerial[make_pair(iDigi, fInfo->CscDigiToChannel(dig, serial))] = serial;
966 }
967
968 return digiChannelToSerial;
969}
970
971map <vector <Int_t>, Long_t> BmnLambdaEmbedding::GetSiliconChannelSerialFromDigi(vector <BmnStripDigit> digits) {
972 map <vector <Int_t>, Long_t> digiChannelToSerial; // (digi index in vector + channel + sample ---> serial)
973
974 for (size_t iDigi = 0; iDigi < digits.size(); iDigi++) {
975 BmnStripDigit* dig = &digits[iDigi];
976
977 Int_t chan = -1, sample = -1;
978 Long_t serial = 0x0;
979
980 // Get channel, sample and serial information ...
981 fInfo->SiliconDigiToChannelSampleSerial(dig, chan, sample, serial);
982 vector <Int_t> tmp{(int)iDigi, chan, sample};
983
984 digiChannelToSerial[tmp] = serial;
985 }
986
987 return digiChannelToSerial;
988}
989
990TString BmnLambdaEmbedding::AddInfoToRawFile(map <UInt_t, vector <BmnStripDigit>> digsGem, map <UInt_t, map <pair <Int_t, Int_t>, Long_t>> evIdChannelSerialGem,
991 map <UInt_t, vector < BmnStripDigit>> digsSilicon, map <UInt_t, map < vector <Int_t>, Long_t>> evIdChannelSerialSilicon,
992 map <UInt_t, vector <BmnStripDigit>> digsCsc, map <UInt_t, map <pair <Int_t, Int_t>, Long_t>> evIdChannelSerialCsc) {
993
994 // Create output file ...
995 TString tmp = fRawRootFileName;
996 tmp = tmp.ReplaceAll(".root", "");
997 tmp += "_withLambdaEmbedded.root";
998
999 TFile* f = new TFile(tmp.Data(), "recreate");
1000 TTree* t = new TTree("BMN_RAW", "BMN_RAW");
1001
1002 TClonesArray* gemOrCsc = new TClonesArray("BmnADCDigit");
1003 TClonesArray* silicon = new TClonesArray("BmnADCDigit");
1004 TClonesArray* sync = new TClonesArray("BmnSyncDigit");
1005
1006 BmnEventHeader* header = new BmnEventHeader();
1007
1008 t->Branch("BmnEventHeader.", &header);
1009 t->Branch("SYNC", &sync);
1010 t->Branch("ADC32", &gemOrCsc);
1011 t->Branch("ADC128", &silicon);
1012
1013 TChain* ch = new TChain("BMN_RAW");
1014 ch->Add(fRawRootFileName.Data());
1015
1016 BmnEventHeader* fEventHeader = nullptr;
1017
1018 ch->SetBranchAddress("ADC32", &fADC32);
1019 ch->SetBranchAddress("ADC128", &fADC128);
1020 ch->SetBranchAddress("SYNC", &fSync);
1021 ch->SetBranchAddress("BmnEventHeader.", &fEventHeader);
1022
1023 // Loop over existing *raw.root file ...
1024 for (Int_t iEntry = 0; iEntry < fEvents; iEntry++) {
1025 if (iEntry % 1000 == 0)
1026 cout << "Embedding, event# " << iEntry << endl;
1027 // Clear TClonesArrays from the previous event ...
1028 gemOrCsc->Delete();
1029 silicon->Delete();
1030 sync->Delete();
1031
1032 ch->GetEntry(iEntry);
1033
1034 UInt_t eventId = fEventHeader->GetEventId();
1035
1036 header->SetRunId(fEventHeader->GetRunId());
1037 header->SetEventId(eventId);
1038 header->SetEventType(fEventHeader->GetEventType());
1039
1040 // Get SYNC digits ....
1041 for (Int_t iSync = 0; iSync < fSync->GetEntriesFast(); iSync++) {
1042 BmnSyncDigit* dig = (BmnSyncDigit*) fSync->UncheckedAt(iSync);
1043
1044 new ((*sync)[sync->GetEntriesFast()]) BmnSyncDigit(dig->GetSerial(), dig->GetEvent(), dig->GetTime_sec(), dig->GetTime_ns());
1045 }
1046
1047 // Looking for current eventId in the digi-map ...
1048 // EventId should be the same for three maps (digsGem, digsSilicon, digsCsc)
1049 Bool_t isEventIdFoundInTheMap = kTRUE;
1050
1051 auto itDigGem = digsGem.find(eventId);
1052 auto itDigSilicon = digsSilicon.find(eventId);
1053 auto itDigCsc = digsCsc.find(eventId);
1054
1055 if (itDigGem == digsGem.end())
1056 isEventIdFoundInTheMap = kFALSE;
1057
1058 // Array with digits for certain eventId
1059 // GEM
1060 vector <BmnStripDigit> digitsGem;
1061 map <pair <Int_t, Int_t>, Long_t> idxChanSerGem;
1062
1063 // SILICON
1064 vector <BmnStripDigit> digitsSilicon;
1065 map <vector <Int_t>, Long_t> idxChanSampleSerSilicon;
1066
1067 // CSC
1068 vector <BmnStripDigit> digitsCsc;
1069 map <pair <Int_t, Int_t>, Long_t> idxChanSerCsc;
1070
1071 if (isEventIdFoundInTheMap) {
1072 // GEM
1073 if (isGemEmbedded) {
1074 digitsGem = itDigGem->second;
1075 idxChanSerGem = evIdChannelSerialGem.find(eventId)->second;
1076 }
1077
1078 // SILICON
1079 if (isSilEmbedded) {
1080 digitsSilicon = itDigSilicon->second;
1081 idxChanSampleSerSilicon = evIdChannelSerialSilicon.find(eventId)->second;
1082 }
1083
1084 // CSC
1085 if (isCscEmbedded) {
1086 digitsCsc = itDigCsc->second;
1087 idxChanSerCsc = evIdChannelSerialCsc.find(eventId)->second;
1088 }
1089 }
1090
1091 // GEM and CSC
1092 if (isGemEmbedded || isCscEmbedded)
1093 for (Int_t iAdc32 = 0; iAdc32 < fADC32->GetEntriesFast(); iAdc32++) {
1094 BmnADCDigit* adc32 = (BmnADCDigit*) fADC32->UncheckedAt(iAdc32);
1095
1096 UInt_t channel = adc32->GetChannel(); // ---> (0 ... 63)
1097 UInt_t serial = adc32->GetSerial();
1098
1099 // GEM:
1100 // Looking for a pair <channel, serial> in the correspondence map ...
1101 if (isGemEmbedded)
1102 for (auto itCorr : idxChanSerGem) {
1103 Int_t channelTmp = itCorr.first.second; // ---> (0 ... 2047)
1104
1105 Int_t channelBuff = Int_t(channelTmp / 32);
1106 Long_t serialBuff = itCorr.second;
1107
1108 if ((Int_t)channel != channelBuff || serial != serialBuff)
1109 continue;
1110
1111 Int_t sample = Int_t(channelTmp % 32);
1112 Int_t signal = Int_t(digitsGem[itCorr.first.first].GetStripSignal());
1113
1114 if (isUseRealSignal)
1115 adc32->GetShortValue()[sample] += signal * 16;
1116 else
1117 adc32->GetShortValue()[sample] += fSignal[0] * 16;
1118
1119 adc32->SetAsEmbedded(kTRUE);
1120 }
1121
1122 // CSC:
1123 // Looking for a pair <channel, serial> in the correspondence map ...
1124 if (isCscEmbedded)
1125 for (auto itCorr : idxChanSerCsc) {
1126 Int_t channelTmp = itCorr.first.second; // ---> (0 ... 2047)
1127
1128 Int_t channelBuff = Int_t(channelTmp / 32);
1129 Long_t serialBuff = itCorr.second;
1130
1131 if ((Int_t)channel != channelBuff || serial != serialBuff)
1132 continue;
1133
1134 Int_t sample = Int_t(channelTmp % 32);
1135 Int_t signal = Int_t(digitsCsc[itCorr.first.first].GetStripSignal() * 16.);
1136
1137 if (isUseRealSignal)
1138 adc32->GetShortValue()[sample] += signal;
1139 else
1140 adc32->GetShortValue()[sample] += fSignal[2] * 16;
1141
1142 adc32->SetAsEmbedded(kTRUE);
1143 }
1144
1145 // Write ADC-digit to the output file ...
1146 new ((*gemOrCsc)[gemOrCsc->GetEntriesFast()]) BmnADCDigit(adc32->GetSerial(), adc32->GetChannel(), adc32->GetNSamples(), adc32->GetShortValue(), adc32->IsEmbedded());
1147 }
1148
1149 // SILICON
1150 if (isSilEmbedded)
1151 for (Int_t iAdc128 = 0; iAdc128 < fADC128->GetEntriesFast(); iAdc128++) {
1152 BmnADCDigit* adc128 = (BmnADCDigit*) fADC128->UncheckedAt(iAdc128);
1153
1154 UInt_t channel = adc128->GetChannel();
1155 UInt_t serial = adc128->GetSerial();
1156
1157 // Looking for an entry with a given <channel, sample, serial> in the correspondence map ...
1158 for (auto itCorr : idxChanSampleSerSilicon) {
1159 Int_t chan = itCorr.first[1];
1160 Long_t ser = itCorr.second;
1161
1162 if ((Int_t)channel != chan || serial != ser)
1163 continue;
1164
1165 Int_t idx = itCorr.first[0];
1166 Int_t signal = Int_t(digitsSilicon[idx].GetStripSignal() * 16.);
1167
1168 Int_t sample = itCorr.first[2];
1169
1170 // x'-layer has a negative polarity!
1171 Int_t sign = (digitsSilicon[idx].GetStripLayer() == 0) ? +1 : -1;
1172 if (isUseRealSignal)
1173 adc128->GetShortValue()[sample] += sign * signal;
1174 else
1175 adc128->GetShortValue()[sample] += sign * fSignal[1] * 16;
1176
1177 adc128->SetAsEmbedded(kTRUE);
1178 }
1179
1180 // Write ADC-digit to the output file ...
1181 new ((*silicon)[silicon->GetEntriesFast()]) BmnADCDigit(adc128->GetSerial(), adc128->GetChannel(), adc128->GetNSamples(), adc128->GetShortValue(), adc128->IsEmbedded());
1182 }
1183
1184 t->Fill();
1185 }
1186
1187 delete ch;
1188
1189 t->Write();
1190 f->Close();
1191
1192 return tmp;
1193}
1194
1195// Create store of possible primary lambdas to be used for embedding ...
1196
1197TClonesArray* BmnLambdaEmbedding::CreateLambdaStore() {
1198 cout << fNstores << endl;
1199 for (Int_t iEntry = 0; iEntry < fSim->GetEntries(); iEntry++) {
1200 fSim->GetEntry(iEntry);
1201 if (iEntry % 1000 == 0) {
1202 cout << "Event# " << iEntry << " processed" << endl;
1203 cout << "Size of store# " << fLambdaStore->GetEntriesFast() << endl;
1204 }
1205
1206 if (fLambdaStore->GetEntriesFast() == fNstores) // Maximum size of current store
1207 break;
1208
1209 for (Int_t iTrack = 0; iTrack < fMCTracks->GetEntriesFast(); iTrack++) {
1210 CbmMCTrack* track = (CbmMCTrack*) fMCTracks->UncheckedAt(iTrack);
1211
1212 // Get primary lambdas only ...
1213 if (track->GetPdgCode() != fPdgDecParticle || track->GetMotherId() != -1)
1214 continue;
1215
1216 // Get lambdas in forward direction ...
1217 if (track->GetPz() < 0.)
1218 continue;
1219
1220 Double_t p = track->GetP();
1221 Double_t Tx = track->GetPx() / track->GetPz();
1222 Double_t Ty = track->GetPy() / track->GetPz();
1223
1224 Double_t Pz = p / TMath::Sqrt(1 + Tx * Tx + Ty * Ty);
1225 Double_t eta = 0.5 * TMath::Log((p + Pz) / (p - Pz));
1226
1227 TVector3 tmp(Tx * Pz, Ty * Pz, Pz);
1228 Double_t phi = tmp.Phi() * TMath::RadToDeg() + 180.;
1229
1230 if (eta < fEtaMin || eta > fEtaMax || p < fMomMin || phi < fPhiMin || phi > fPhiMax)
1231 continue;
1232
1233 new ((*fLambdaStore)[fLambdaStore->GetEntriesFast()]) BmnParticleStore(fPdgDecParticle, p, Tx, Ty);
1234 }
1235 }
1236 return fLambdaStore;
1237}
1238
1239void BmnLambdaEmbedding::PrintStoreInfo() {
1240 cout << "N_particles_store# " << fLambdaStore->GetEntriesFast() << endl;
1241
1242 for (Int_t iLambda = 0; iLambda < fLambdaStore->GetEntriesFast(); iLambda++) {
1243 BmnParticleStore* lambda = (BmnParticleStore*) fLambdaStore->UncheckedAt(iLambda);
1244 cout << "P = " << lambda->GetP() << " Eta = " << lambda->GetEta() << " Phi = " << lambda->GetPhi() << endl;
1245 }
1246}
1247
1248BmnLambdaEmbedding::~BmnLambdaEmbedding() {
1249 delete fSim;
1250 delete fReco;
1251
1252 delete fInfo;
1253
1254 if (fMon)
1255 delete fMon;
1256
1257 if (normUtils)
1258 delete normUtils;
1259}
1260
1261void BmnLambdaEmbedding::SimulateLambdaPassing(Int_t pdg, Double_t P, TVector2 pos, TVector3 Vp, Int_t iLambda, Int_t iVertex) {
1262 // Here there is a list of arguments to be passed to the macro below ...
1263
1264 TString namePart = (pdg == 3122) ? "lambda" : (pdg == 310) ? "kaon" : "unknown";
1265 TString outFile = TString::Format("%s/%s%d_vertex%d.root", fStorePath.Data(), namePart.Data(), iLambda, iVertex);
1266
1267 TString macroName = "SimulateLambdaPassing.C";
1268
1269 TString gPathConfig = gSystem->Getenv("VMCWORKDIR");
1270 TString gPathFull = gPathConfig + "/macro/embedding/" + macroName;
1271
1272 TString argList = TString::Format("(%d, %d, %G, %G, %G, %G, %G, %G, %G, \"%s\")", pdg, fNeventsToBeSimulated, fFieldScale, P, pos.X(), pos.Y(), Vp.X(), Vp.Y(), Vp.Z(), outFile.Data());
1273 TString exeCommand = TString::Format("root -b -q '%s%s'", gPathFull.Data(), argList.Data());
1274
1275 gSystem->Exec(exeCommand.Data());
1276}
1277
1278Int_t BmnLambdaEmbedding::DefineSiliconStatByZpoint(Double_t z) {
1279 // Boardings are valid for RUN7 only!!!
1280 const Double_t z1 = 14.;
1281 const Double_t z2 = 20.;
1282
1283 Int_t shift = 6; // nGems in RUN7
1284 Int_t tmp;
1285
1286 if (z < z1)
1287 tmp = 0;
1288 else if (z > z1 && z < z2)
1289 tmp = 1;
1290 else if (z > z2)
1291 tmp = 2;
1292 else
1293 return -1;
1294
1295 return shift + tmp;
1296}
f
float f
Definition P4_F32vec4.h:21
kBC
@ kBC
Definition BmnDetectorList.h:31
kSILICON
@ kSILICON
Definition BmnDetectorList.h:25
kGEM
@ kGEM
Definition BmnDetectorList.h:15
kBMNSETUP
@ kBMNSETUP
Definition BmnEnums.h:90
BmnADCDigit::GetNSamples
UInt_t GetNSamples() const
Definition BmnADCDigit.h:35
BmnADCDigit::GetChannel
UShort_t GetChannel() const
Definition BmnADCDigit.h:33
BmnADCDigit::GetSerial
UInt_t GetSerial() const
Definition BmnADCDigit.h:31
BmnADCDigit::SetAsEmbedded
void SetAsEmbedded(Bool_t flag)
Definition BmnADCDigit.h:53
BmnADCDigit::IsEmbedded
Bool_t IsEmbedded()
Definition BmnADCDigit.h:55
BmnADCDigit::GetShortValue
Short_t * GetShortValue() const
Definition BmnADCDigit.h:39
BmnDecoder::SetBmnSetup
void SetBmnSetup(BmnSetup v)
Definition BmnDecoder.h:235
BmnDecoder::SetRunId
void SetRunId(UInt_t v)
Definition BmnDecoder.h:124
BmnDecoder::SetRawRootFile
void SetRawRootFile(TString filename)
Definition BmnDecoder.h:270
BmnDecoder::SetDetectorSetup
void SetDetectorSetup(std::map< DetectorId, bool > setup)
Definition BmnDecoder.h:229
BmnDecoder::SetDigiRootFile
void SetDigiRootFile(TString filename)
Definition BmnDecoder.h:272
BmnDecoder::SetMSCMapping
void SetMSCMapping(TString map)
Definition BmnDecoder.h:225
BmnDecoder::SetPeriodId
void SetPeriodId(UInt_t v)
Definition BmnDecoder.h:126
BmnDecoder::SetTrigPlaceMapping
void SetTrigPlaceMapping(TString map)
Definition BmnDecoder.h:160
BmnDecoder::SetGemMapping
void SetGemMapping(TString map)
Definition BmnDecoder.h:174
BmnDecoder::DecodeDataToDigi
BmnStatus DecodeDataToDigi()
Definition BmnDecoder.cxx:341
BmnDecoder::SetTrigChannelMapping
void SetTrigChannelMapping(TString file)
Definition BmnDecoder.h:168
BmnDecoder::SetCSCMapping
void SetCSCMapping(TString map)
Definition BmnDecoder.h:166
BmnDecoder::SetSiliconMapping
void SetSiliconMapping(TString map)
Definition BmnDecoder.h:162
BmnEventHeader::GetEventId
UInt_t GetEventId()
Definition BmnEventHeader.h:91
BmnEventHeader::GetEventType
BmnEventType GetEventType()
Definition BmnEventHeader.h:104
BmnLambdaEmbedding::SetUseRealSignals
void SetUseRealSignals(Bool_t flag)
Definition BmnLambdaEmbedding.h:224
BmnLambdaEmbedding::DoRawRootConvertion
void DoRawRootConvertion(Bool_t flag)
Definition BmnLambdaEmbedding.h:194
BmnLambdaEmbedding::DoSimulateLambdaThroughSetup
void DoSimulateLambdaThroughSetup(Bool_t flag)
Definition BmnLambdaEmbedding.h:190
BmnLambdaEmbedding::DoDecode
void DoDecode(Bool_t flag)
Definition BmnLambdaEmbedding.h:202
BmnLambdaEmbedding::DoLambdaStore
void DoLambdaStore(Bool_t flag)
Definition BmnLambdaEmbedding.h:182
BmnLambdaEmbedding::DoEmbeddingMonitor
void DoEmbeddingMonitor(Bool_t flag)
Definition BmnLambdaEmbedding.h:210
BmnLambdaEmbedding::CreateLambdaStore
TClonesArray * CreateLambdaStore()
Definition BmnLambdaEmbedding.cxx:1197
BmnLambdaEmbedding::DoListOfEventsWithReconstructedVertex
void DoListOfEventsWithReconstructedVertex(Bool_t flag)
Definition BmnLambdaEmbedding.h:186
BmnLambdaEmbedding::DoPrintStoreInfo
void DoPrintStoreInfo(Bool_t flag)
Definition BmnLambdaEmbedding.h:206
BmnLambdaEmbedding::DoEmbedding
void DoEmbedding(Bool_t flag)
Definition BmnLambdaEmbedding.h:198
BmnLambdaMisc::GemDigiToChannel
Int_t GemDigiToChannel(BmnStripDigit *, Long_t &)
Definition BmnLambdaMisc.cxx:408
BmnLambdaMisc::SiliconDigiToChannelSampleSerial
void SiliconDigiToChannelSampleSerial(BmnStripDigit *, Int_t &, Int_t &, Long_t &)
Definition BmnLambdaMisc.cxx:317
BmnLambdaMisc::CscDigiToChannel
Int_t CscDigiToChannel(BmnStripDigit *, Long_t &)
Definition BmnLambdaMisc.cxx:273
BmnLambdaMisc::GemDigiChannelToSerial
Long_t GemDigiChannelToSerial(pair< BmnStripDigit, Int_t >)
Definition BmnLambdaMisc.cxx:540
BmnMatch::GetMatchedLink
const BmnLink & GetMatchedLink() const
Definition BmnMatch.h:39
BmnParticlePair::GetNHitsPart1
Int_t GetNHitsPart1(TString det="")
Definition BmnParticlePair.h:134
BmnParticlePair::SetTxPair
void SetTxPair(Double_t val1, Double_t val2)
Definition BmnParticlePair.h:80
BmnParticlePair::GetNHitsPart2
Int_t GetNHitsPart2(TString det="")
Definition BmnParticlePair.h:148
BmnParticlePair::GetMomPart1
Double_t GetMomPart1()
Definition BmnParticlePair.h:190
BmnParticlePair::GetTxPart1
Double_t GetTxPart1()
Definition BmnParticlePair.h:218
BmnParticlePair::GetTyPart1
Double_t GetTyPart1()
Definition BmnParticlePair.h:226
BmnParticlePair::SetNHitsPair
void SetNHitsPair(vector< Int_t > part1, vector< Int_t > part2)
Definition BmnParticlePair.h:118
BmnParticlePair::SetTyPair
void SetTyPair(Double_t val1, Double_t val2)
Definition BmnParticlePair.h:85
BmnParticlePair::SetMomPair
void SetMomPair(Double_t val1, Double_t val2)
Definition BmnParticlePair.h:70
BmnParticlePair::GetTxPart2
Double_t GetTxPart2()
Definition BmnParticlePair.h:222
BmnParticlePair::GetMomPart2
Double_t GetMomPart2()
Definition BmnParticlePair.h:194
BmnParticlePair::GetTyPart2
Double_t GetTyPart2()
Definition BmnParticlePair.h:230
BmnRawDataDecoder::SetPeriodId
void SetPeriodId(UInt_t v)
Definition BmnRawDataDecoder.h:93
BmnRawDataDecoder::GetDecoder
BmnDecoder * GetDecoder() const
Definition BmnRawDataDecoder.h:193
BmnRawDataDecoder::SetRunId
void SetRunId(UInt_t v)
Definition BmnRawDataDecoder.h:84
BmnSyncDigit::GetSerial
UInt_t GetSerial() const
Definition BmnSyncDigit.h:17
BmnSyncDigit::GetTime_sec
Long64_t GetTime_sec() const
Definition BmnSyncDigit.h:25
BmnSyncDigit::GetTime_ns
Long64_t GetTime_ns() const
Definition BmnSyncDigit.h:29
BmnSyncDigit::GetEvent
Long64_t GetEvent() const
Definition BmnSyncDigit.h:21
CbmMCTrack::GetPy
Double_t GetPy() const
Definition CbmMCTrack.h:59
CbmMCTrack::GetPz
Double_t GetPz() const
Definition CbmMCTrack.h:60
CbmMCTrack::GetPx
Double_t GetPx() const
Definition CbmMCTrack.h:58
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition CbmMCTrack.h:57
CbmMCTrack::GetP
Double_t GetP() const
Definition CbmMCTrack.h:68
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition CbmMCTrack.h:56
CbmStsPoint::GetStation
Int_t GetStation() const
Definition CbmStsPoint.h:81
CbmVertex::GetZ
Double_t GetZ() const
Definition CbmVertex.h:60
CbmVertex::GetX
Double_t GetX() const
Definition CbmVertex.h:58
CbmVertex::GetY
Double_t GetY() const
Definition CbmVertex.h:59
DstEventHeader::GetEventId
UInt_t GetEventId()
Definition DstEventHeader.h:53
SignalNormalizationUtils::SetMcDataSet
void SetMcDataSet(TChain *chainMC)
Definition BmnLambdaEmbedding.h:103
SignalNormalizationUtils::GetRescaleFunction
TF1 * GetRescaleFunction(TString det)
Definition BmnLambdaEmbedding.h:119
UniRun::GetRun
static UniRun * GetRun(int period_number, int run_number)
get run from the database
Definition UniRun.cxx:177
UniRun::GetBeamParticle
TString GetBeamParticle()
get beam particle of the current run
Definition UniRun.h:113
UniRun::GetFieldVoltage
double * GetFieldVoltage()
get field voltage of the current run
Definition UniRun.h:125
centrality_clusterization.int
int
Definition centrality_clusterization.py:210
publish_test_results.header
dict header
Definition publish_test_results.py:7
run
-clang-format
setup
Definition setup.py:1