CbmRoot
Loading...
Searching...
No Matches
PairAnalysis.cxx
Go to the documentation of this file.
1
2// PairAnalysis Main Analysis class //
3// //
4// Authors:
5// * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
6// Julian Book <Julian.Book@cern.ch>
7/*
8
9Framework to perform event, single track and pair selections.
10
11Convention for the signs of the track in fTracks:
12The names are available via the function TrackClassName(Int_t i)
13
14 0: SE+ (same event +)
15 1: SE- (same event -)
16
17
18Convention for the signs of the pair in fPairCandidates:
19The names are available via the function PairClassName(Int_t i)
20
21 0: SE+ SE+ (same event like sign +)
22 1: SE+ SE- (same event unlike sign)
23 2: SE- SE- (same event like sign -)
24
25 3: SE+ ME+ (mixed event like sign +)
26 4: SE- ME+ (mixed event unlike sign -+)
27 5: SE+ ME- (mixed event unlike sign +-)
28 6: SE- ME- (mixed event like sign -)
29
30 7: SE+ SE- (same event track rotation)
31
32
33Some options to configure your analysis instance:
34 SetLegPdg(Int_t pdgLeg1, Int_t pdgLeg2) // define pair production of interest
35 SetRefitWithMassAssump(kTRUE) // whether tracks should be refitted accordingly
36 SetNoPairing() // look only at track level
37 SetProcessLS(kFALSE) // switch off like-sign pair calculations
38 SetUseKF(kTRUE) // use Kalman-Filter mathematics for vertexing
39
40
41Some options to configure the output:
42 SetHistogramManager(PairAnalysisHistos * const histos) // default histogram manager
43 SetHistogramArray(PairAnalysisHF * const histoarray) // histogram matrix
44 SetCutQA() // switch on basic qa histograms
45
46
47Some options to add background estimators:
48 SetMixingHandler(PairAnalysisMixingHandler *mix) // event mixing
49 SetTrackRotator(PairAnalysisTrackRotator * const rot) // track rotation
50
51
52*/
53// //
55
56#include "CbmHit.h"
57#include "CbmMCTrack.h"
58#include "CbmMatch.h"
59#include "CbmRichRing.h"
60#include "CbmTrack.h"
61
62#include "FairMCPoint.h"
63
64#include <TDatabasePDG.h>
65#include <TGrid.h>
66#include <TList.h>
67#include <TLorentzVector.h>
68#include <TMath.h>
69#include <TObject.h>
70#include <TRandom3.h>
71#include <TString.h>
72
73#include "PairAnalysisEvent.h"
74#include "PairAnalysisHistos.h"
75#include "PairAnalysisMC.h"
76#include "PairAnalysisPair.h"
77#include "PairAnalysisPairKF.h"
78#include "PairAnalysisPairLV.h"
79#include "PairAnalysisTrack.h"
82//#include "PairAnalysisDebugTree.h"
86//#include "PairAnalysisV0Cuts.h"
87#include "PairAnalysis.h"
88#include "PairAnalysisHistos.h"
89
91
92 const char* PairAnalysis::fgkTrackClassNames[2] = {"+", "-"};
93
94const char* PairAnalysis::fgkPairClassNames[8] = {"SE++", "SE+-", "SE--", "ME++", "ME-+", "ME+-", "ME--", "TR+-"};
95
96//________________________________________________________________
97PairAnalysis::PairAnalysis() : PairAnalysis("PairAnalysis", "PairAnalysis")
98{
99 //
100 // Default constructor
101 //
102}
103
104//________________________________________________________________
105PairAnalysis::PairAnalysis(const char* name, const char* title)
106 : TNamed(name, title)
107 , fEventFilter("EventFilter")
108 , fTrackFilter("TrackFilter")
109 , fPairPreFilterLegs("PairPreFilterLegs")
110 , fPairPreFilter("PairPreFilter")
111 , fFinalTrackFilter("FinalTrackFilter")
112 , fPairFilter("PairFilter")
113 , fTrackFilterMC("TrackFilterMC")
114 , fPairFilterMC("PairFilterMC")
115 , fUsedVars(new TBits(PairAnalysisVarManager::kNMaxValuesMC))
116 , fPairCandidates(new TObjArray(8))
117{
118 //
119 // Named constructor
120 //
121}
122
123//________________________________________________________________
125{
126 //
127 // Default destructor
128 //
129 if (fQAmonitor) delete fQAmonitor;
130 if (fHistos) delete fHistos;
131 if (fUsedVars) delete fUsedVars;
133 if (fMixing) delete fMixing;
134 if (fSignalsMC) delete fSignalsMC;
135 if (fHistoArray) delete fHistoArray;
136 if (fCutStepHistos) delete fCutStepHistos;
137}
138
139//________________________________________________________________
141{
142 //
143 // Initialise objects
144 //
145
146 // if(GetHasMC()) PairAnalysisMC::Instance()->SetHasMC(GetHasMC());
147
149
150 // compress the MC signal array
151 if (fSignalsMC) fSignalsMC->Compress();
152
153 if (fMixing) fMixing->Init(this);
154 if (fHistoArray) {
155 // fHistoArray->SetSignalsMC(fSignalsMC);
156 fHistoArray->Init();
157 }
158
159 // for internal train wagons
160 if (!fEventProcess) {
161 PairAnalysisPairLegCuts* trk2leg = new PairAnalysisPairLegCuts("trk2leg", "trk2leg");
162 // move all track cuts (if any) into pair leg cuts
163 TIter listIterator(fTrackFilter.GetCuts());
164 while (AnalysisCuts* thisCut = (AnalysisCuts*) listIterator()) {
165 trk2leg->GetLeg1Filter().AddCuts((AnalysisCuts*) thisCut->Clone());
166 trk2leg->GetLeg2Filter().AddCuts((AnalysisCuts*) thisCut->Clone());
167 }
168 TIter listIterator2(fFinalTrackFilter.GetCuts());
169 while (AnalysisCuts* thisCut = (AnalysisCuts*) listIterator2()) {
170 trk2leg->GetLeg1Filter().AddCuts((AnalysisCuts*) thisCut->Clone());
171 trk2leg->GetLeg2Filter().AddCuts((AnalysisCuts*) thisCut->Clone());
172 }
173 // add pair leg cuts to pair filter
174 fPairFilter.AddCuts(trk2leg);
175 }
176
177 // initialize simple cut qa
178 if (fCutQA) {
179 fQAmonitor = new PairAnalysisCutQa(Form("QAcuts_%s", GetName()), "QAcuts");
186 fQAmonitor->Init();
187 }
188
189 if (fHistos) { (*fUsedVars) |= (*fHistos->GetUsedVars()); }
190
191 // initialize cut step histograms
195 else
198 else
200 fCutStepHistos->SetName(Form("CutSteps_%s", GetName()));
201 }
202}
203
204//________________________________________________________________
205
206void PairAnalysis::Process(TObjArray* arr)
207{
208 //
209 // Process the pair array
210 //
211
212 // set pair arrays
213 fPairCandidates = arr;
214
215 // fill pair and pair leg histograms
216 if (fHistos) FillHistograms(0x0, kTRUE);
217
218 // apply cuts and fill output TODO: OBSOLETE!
219 // if (fHistos) FillHistogramsFromPairArray();
220
222}
223
224//________________________________________________________________
226{
227 //
228 // Process the events
229 //
230
232 if (!ev1->GetPrimaryVertex()) {
233 Fatal("Process", "No vertex found!");
234 //Error("Process","No vertex found!");
235 return kFALSE;
236 }
237
240 TDatabasePDG::Instance()->GetParticle(fPdgMother)->Mass());
241
245
247 if (fMixing) {
250 }
251
252
253 //in case we have MC load the MC event and process the MC particles
254 // why do not apply the event cuts first ????
255 // if (PairAnalysisMC::Instance()->ConnectMCEvent()){
256 if (PairAnalysisMC::Instance()->HasMC()) { ProcessMC(); }
257
259 if (!fPairCandidates->UncheckedAt(0)) { InitPairCandidateArrays(); }
260 else {
261 ClearArrays();
262 }
263
265 UInt_t selectedMask = (1 << fEventFilter.GetCuts()->GetEntries()) - 1;
266
268 UInt_t cutmask = fEventFilter.IsSelected(ev1);
269 if (fCutQA) fQAmonitor->FillAll(ev1);
270 if (fCutQA) fQAmonitor->Fill(cutmask, ev1);
271 if (ev1 && cutmask != selectedMask) return kFALSE;
272
274 FillTrackArrays(ev1);
275
276 // prefilter track
277 if ((fPreFilterAllSigns || fPreFilterUnlikeOnly) && (fPairPreFilter.GetCuts()->GetEntries() > 0))
278 PairPreFilter(0, 1, fTracks[0], fTracks[1]);
279
280 // remove tracks from arrays that DO NOT pass 2nd cut iteration
282
283 // create SE pairs and fill pair candidate arrays
284 if (!fNoPairing) {
285 for (Int_t itrackArr1 = 0; itrackArr1 < 2; ++itrackArr1) {
286 for (Int_t itrackArr2 = itrackArr1; itrackArr2 < 2; ++itrackArr2) {
287 if (!fProcessLS && GetPairIndex(itrackArr1, itrackArr2) != static_cast<Int_t>(EPairType::kSEPM)) continue;
288 FillPairArrays(itrackArr1, itrackArr2);
289 }
290 }
291 // add track rotated pairs
293 }
294
295 //process event mixing
296 if (fMixing) {
297 fMixing->Fill(ev1, this);
298 // FillHistograms(0x0,kTRUE);
299 }
300
301 // fill candidate variables
302 Double_t ntracks = fTracks[0].GetEntriesFast() + fTracks[1].GetEntriesFast();
303 Double_t npairs = PairArray(static_cast<Int_t>(EPairType::kSEPM))->GetEntriesFast();
306
307 //in case there is a histogram manager, fill the QA histograms
308 if (fHistos) FillHistograms(ev1);
309 // fill histo array with event information only
310 if (fHistoArray) fHistoArray->Fill(0, const_cast<Double_t*>(PairAnalysisVarManager::GetData()), 0x0, 0x0);
311
312 // clear arrays
314
315 return kTRUE;
316}
317
318//________________________________________________________________
320{
321 //
322 // Process the MC data
323 //
325
326 // fill mc true and rec event
327 // if (fHistos) FillHistogramsMC(papaMC->GetMCEvent(), ev1);
328
329 // there are mc tracks
330 if (!papaMC->GetNMCTracks()) return;
331
332 // signals to be studied
333 if (!fSignalsMC) return;
334 Int_t nSignals = fSignalsMC->GetEntriesFast();
335 if (!nSignals) return;
336
337 //loop over all MC data and Fill the HF, CF containers and histograms if they exist
338
339 Bool_t bFillHF = kFALSE;
340 Bool_t bFillHist = kFALSE;
341 for (Int_t isig = 0; isig < nSignals; isig++) {
342 PairAnalysisSignalMC* sigMC = (PairAnalysisSignalMC*) fSignalsMC->UncheckedAt(isig);
343 if (!sigMC->GetFillPureMCStep()) continue;
344 TString sigName = sigMC->GetName();
345 if (fHistos && !bFillHist) {
346 bFillHist |= fHistos->HasHistClass(Form("Pair_%s_MCtruth", sigName.Data()));
347 bFillHist |= fHistos->HasHistClass(Form("Track.Leg_%s_MCtruth", sigName.Data()));
348 bFillHist |= fHistos->HasHistClass(Form("Track.%s_%s_MCtruth", fgkPairClassNames[1], sigName.Data()));
349 }
350 if (fHistoArray && !bFillHF) {
351 bFillHF |= fHistoArray->HasHistClass(Form("Pair_%s_MCtruth", sigName.Data()));
352 bFillHF |= fHistoArray->HasHistClass(Form("Track.Leg_%s_MCtruth", sigName.Data()));
353 bFillHF |= fHistoArray->HasHistClass(Form("Track.%s_%s_MCtruth", fgkPairClassNames[1], sigName.Data()));
354 }
355 }
356
357 // check if there is anything to fill
358 if (!bFillHF && !bFillHist) return;
359
360 // initialize 2D arrays of labels for particles from each MC signal
361 Int_t** labels1; // labels for particles satisfying branch 1
362 Int_t** labels2; // labels for particles satisfying branch 2
363 Int_t** labels12; // labels for particles satisfying both branches
364 labels1 = new Int_t*[nSignals];
365 labels2 = new Int_t*[nSignals];
366 labels12 = new Int_t*[nSignals];
367 Int_t* indexes1 = new Int_t[nSignals];
368 Int_t* indexes2 = new Int_t[nSignals];
369 Int_t* indexes12 = new Int_t[nSignals];
370 for (Int_t isig = 0; isig < nSignals; ++isig) {
371 *(labels1 + isig) = new Int_t[papaMC->GetNMCTracks()];
372 *(labels2 + isig) = new Int_t[papaMC->GetNMCTracks()];
373 *(labels12 + isig) = new Int_t[papaMC->GetNMCTracks()];
374 for (Int_t ip = 0; ip < papaMC->GetNMCTracks(); ++ip) {
375 labels1[isig][ip] = -1;
376 labels2[isig][ip] = -1;
377 labels12[isig][ip] = -1;
378 }
379 indexes1[isig] = 0;
380 indexes2[isig] = 0;
381 indexes12[isig] = 0;
382 }
383
384 Bool_t truth1 = kFALSE;
385 Bool_t truth2 = kFALSE;
386 // selection of particles
387 UInt_t selectedMask = (1 << fTrackFilterMC.GetCuts()->GetEntries()) - 1;
388 // get event data (this contains MC event informations as well)
389 Double_t* values = PairAnalysisVarManager::GetData();
390 // loop over the MC tracks
391 for (Int_t ipart = 0; ipart < papaMC->GetNMCTracks(); ++ipart) {
392
393 //get MC particle
394 CbmMCTrack* mctrk = papaMC->GetMCTrackFromMCEvent(ipart);
395
396 // fill variables
397 PairAnalysisVarManager::Fill(mctrk, values);
398
399 //apply track cuts
400 UInt_t cutmask = fTrackFilterMC.IsSelected(values);
401 //fill cut QA
402 if (fCutQA) fQAmonitor->FillAll(mctrk);
403 if (fCutQA) fQAmonitor->Fill(cutmask, mctrk);
404
406 if (fTrackFilterMC.GetHistogramList()->GetSize()) FillCutStepHistogramsMC(&fTrackFilterMC, cutmask, ipart, values);
407
408 // rejection
409 if (cutmask != selectedMask) continue;
410
411 // loop over signals
412 for (Int_t isig = 0; isig < nSignals; ++isig) {
413 PairAnalysisSignalMC* sigMC = (PairAnalysisSignalMC*) fSignalsMC->UncheckedAt(isig);
414 // NOTE: Some signals can be satisfied by many particles and this leads to high
415 // computation times (e.g. secondary electrons from the GEANT transport). Be aware of this!!
416
417 // Proceed only if this signal is required in the pure MC step
418 if (!sigMC->GetFillPureMCStep()) continue;
419
420 truth1 = papaMC->IsMCTruth(ipart, (PairAnalysisSignalMC*) fSignalsMC->UncheckedAt(isig), 1);
422 truth2 = papaMC->IsMCTruth(ipart, (PairAnalysisSignalMC*) fSignalsMC->UncheckedAt(isig), 2);
423
424 // particles satisfying both branches are treated separately to avoid double counting during pairing
425 if (truth1 && truth2) {
426 labels12[isig][indexes12[isig]] = ipart;
427 ++indexes12[isig];
428 }
429 else {
430 if (truth1) {
431 labels1[isig][indexes1[isig]] = ipart;
432 ++indexes1[isig];
433 }
434 if (truth2) {
435 labels2[isig][indexes2[isig]] = ipart;
436 ++indexes2[isig];
437 }
438 }
439 }
440 } // end loop over MC particles
441
446
447 selectedMask = (1 << fPairFilterMC.GetCuts()->GetEntries()) - 1;
448 // loop over signals
449 for (Int_t isig = 0; isig < nSignals; ++isig) {
450
451 // mix the particles which satisfy only one of the signal branches
452 for (Int_t i1 = 0; i1 < indexes1[isig]; ++i1) {
453 CbmMCTrack* part1 = papaMC->GetMCTrackFromMCEvent(labels1[isig][i1]);
454 Int_t mLabel1 = part1->GetMotherId();
455
456 // (e.g. single electrons only, no pairs)
457 if (!indexes2[isig]) FillMCHistograms(labels1[isig][i1], -1, isig);
458 // if(sigMC->IsSingleParticle()) continue;
459
464 for (Int_t i2 = i1; i2 < indexes2[isig]; ++i2) {
465 // for(Int_t i2=0;i2<indexes2[isig];++i2) { /// old slower way
466
467 CbmMCTrack* part2 = papaMC->GetMCTrackFromMCEvent(labels2[isig][i2]);
468 Int_t mLabel2 = part2->GetMotherId();
469 CbmMCTrack* mother = 0x0;
470
472 if (mLabel1 == mLabel2) mother = papaMC->GetMCTrackFromMCEvent(mLabel1);
473
474 // selection of MCtruth pairs
475 UInt_t cutmask = 0;
476 if (mother) cutmask = fPairFilterMC.IsSelected(mother);
477 else {
482 pair->SetMCTracks(part1, part2);
483 cutmask = fPairFilterMC.IsSelected(pair);
484 }
485 //apply MC truth pair cuts
486 if (cutmask != selectedMask) continue;
487
488 if (bFillHF) fHistoArray->Fill(labels1[isig][i1], labels2[isig][i2], isig);
489 if (bFillHist) {
490 if (FillMCHistograms(labels1[isig][i1], labels2[isig][i2], isig)) break;
491 }
492 }
493 }
494
495 // mix the particles which satisfy both branches
496 for (Int_t i1 = 0; i1 < indexes12[isig]; ++i1) {
497 for (Int_t i2 = 0; i2 < i1; ++i2) {
498
499 CbmMCTrack* part1 = papaMC->GetMCTrackFromMCEvent(labels12[isig][i1]);
500 CbmMCTrack* part2 = papaMC->GetMCTrackFromMCEvent(labels12[isig][i2]);
501 Int_t mLabel1 = part1->GetMotherId();
502 Int_t mLabel2 = part2->GetMotherId();
503
505 CbmMCTrack* mother = 0x0;
506 if (mLabel1 == mLabel2) mother = papaMC->GetMCTrackFromMCEvent(mLabel1);
507
508 // selection of MCtruth pairs
509 UInt_t cutmask = 0;
510 if (mother) cutmask = fPairFilterMC.IsSelected(mother);
511 else {
513 pair->SetMCTracks(part1, part2);
514 cutmask = fPairFilterMC.IsSelected(pair);
515 }
516 //apply MC truth pair cuts
517 if (cutmask != selectedMask) continue;
518
519 if (bFillHF) fHistoArray->Fill(labels12[isig][i1], labels12[isig][i2], isig);
520 if (bFillHist) {
521 if (FillMCHistograms(labels12[isig][i1], labels12[isig][i2], isig)) break;
522 }
523 }
524 }
525 } // end loop over signals
526
527 // release the memory
528 delete pair;
529 for (Int_t isig = 0; isig < nSignals; ++isig) {
530 delete[] * (labels1 + isig);
531 delete[] * (labels2 + isig);
532 delete[] * (labels12 + isig);
533 }
534 delete[] labels1;
535 delete[] labels2;
536 delete[] labels12;
537 delete[] indexes1;
538 delete[] indexes2;
539 delete[] indexes12;
540}
541
542//________________________________________________________________
544{
545 //
546 // Fill Histogram information for tracks after prefilter
547 // ignore mixed events - for prefilter, only single tracks +/- are relevant
548 //
549
550 TString className, className2;
551 Double_t* values = PairAnalysisVarManager::GetData();
553
554 //Fill track information, separately for the track array candidates
555 for (Int_t i = 0; i < 2; ++i) {
556 className.Form("Pre_%s", fgkTrackClassNames[i]);
557 if (!fHistos->GetHistogramList()->FindObject(className.Data())) continue;
558 Int_t ntracks = tracks[i]->GetEntriesFast();
559 for (Int_t itrack = 0; itrack < ntracks; ++itrack) {
560 PairAnalysisVarManager::Fill(tracks[i]->UncheckedAt(itrack), values);
561 fHistos->FillClass(className, values);
562 }
563 }
564}
565
566
567//________________________________________________________________
569{
570 //
571 // Fill Histogram information for MCEvents
572 //
573 return;
574 /*
575 Double_t *values=PairAnalysisVarManager::GetData();
576 PairAnalysisVarManager::SetFillMap(fUsedVars);
577
578 // Fill event information
579 PairAnalysisVarManager::Fill(ev1, values); // ESD/AOD information
580 PairAnalysisVarManager::Fill(ev, values); // MC truth info
581 if (fHistos->GetHistogramList()->FindObject("MCEvent"))
582 fHistos->FillClass("MCEvent", values);
583 */
584}
585
586
587//________________________________________________________________
588void PairAnalysis::FillHistograms(const PairAnalysisEvent* ev, Bool_t pairInfoOnly)
589{
590 //
591 // Fill Histogram information for event, pairs, tracks, hits
592 //
593
594 TString className, className2, className3;
595 TString sigName;
596 Double_t* values = PairAnalysisVarManager::GetData(); //NEW CHANGED
598
599 //Fill event information
600 if (ev) {
601 if (fHistos && fHistos->HasHistClass("Event")) fHistos->FillClass("Event", values);
602 if (fHistoArray && fHistoArray->HasHistClass("Event")) fHistoArray->FillClass("Event", values);
603 }
604
605 //Fill track information, separately for the track array candidates
606 Int_t nsig = (fSignalsMC ? fSignalsMC->GetEntriesFast() : 0);
607 TBits trkClassMC(nsig);
608 TBits trkClassMChf(nsig);
609 TBits fillMC(nsig);
610 PairAnalysisMC* mc = (nsig ? PairAnalysisMC::Instance() : 0x0);
611 PairAnalysisSignalMC* sigMC = 0x0;
612 if (!pairInfoOnly) {
613 className2.Form("Track.%s", fgkPairClassNames[1]); // unlike sign, SE only
614 Bool_t mergedtrkClass = fHistos->HasHistClass(className2);
615 Bool_t mergedtrkClass2 = (fHistoArray && fHistoArray->HasHistClass(className2));
616 // check mc signal filling
617 for (Int_t isig = 0; isig < nsig; isig++) {
618 sigName = className2 + "_" + fSignalsMC->At(isig)->GetName();
619 trkClassMC.SetBitNumber(isig, fHistos->HasHistClass(sigName));
620 trkClassMChf.SetBitNumber(isig, fHistoArray && fHistoArray->HasHistClass(sigName));
621 }
622 // loop over both track arrays
623 for (Int_t i = 0; i < fLegTypes; ++i) {
624 className.Form("Track.%s", fgkTrackClassNames[i]);
625 Bool_t trkClass = fHistos->HasHistClass(className);
626 Bool_t trkClass2 = (fHistoArray && fHistoArray->HasHistClass(className));
627 Bool_t fill = (mergedtrkClass || mergedtrkClass2 || trkClass || trkClass2);
628 // skip stuff not to be filled
629 if (!fill && !trkClassMC.CountBits() && !trkClassMChf.CountBits()) continue;
630 Int_t ntracks = fTracks[i].GetEntriesFast();
631 // loop over all tracks
632 for (Int_t itrack = 0; itrack < ntracks; ++itrack) {
633 PairAnalysisTrack* track = static_cast<PairAnalysisTrack*>(fTracks[i].UncheckedAt(itrack));
634 PairAnalysisVarManager::Fill(track, values);
635 if (trkClass) fHistos->FillClass(className, values);
636 if (trkClass2) fHistoArray->FillClass(className, values);
637 if (mergedtrkClass) fHistos->FillClass(className2, values);
638 if (mergedtrkClass2) fHistoArray->FillClass(className2, values); //TODO: check only SE?
639 // check and do mc signal filling
640 for (Int_t isig = 0; isig < nsig; isig++) {
641 if (!trkClassMC.TestBitNumber(isig) && !trkClassMChf.TestBitNumber(isig)) continue;
642 // test if track does correspond to the signal
643 fillMC.SetBitNumber(isig, mc->IsMCTruth(track, (PairAnalysisSignalMC*) fSignalsMC->At(isig), 1)
644 || mc->IsMCTruth(track, (PairAnalysisSignalMC*) fSignalsMC->At(isig), 2));
645 sigName = className2 + "_" + fSignalsMC->At(isig)->GetName();
646 // fill histos
647 if (fillMC.TestBitNumber(isig)) {
648 if (trkClassMC.TestBitNumber(isig)) fHistos->FillClass(sigName, values);
649 if (trkClassMChf.TestBitNumber(isig)) fHistoArray->FillClass(sigName, values);
650 }
651 }
652
653 //Fill tracks hit information
654 FillHistogramsHits(ev, &fillMC, track, kFALSE, values);
655
656 } // track loop
657 } // loop leg type
658 } // not pair info only
659
660 //Fill Pair information, separately for all pair candidate arrays and the legs
661 TBits legClassMC(nsig);
662 TBits legClassMChf(nsig);
663 TBits pairClassMC(nsig);
664 TBits pairClassMChf(nsig);
665 TObjArray arrLegs(100);
666 for (Int_t i = 0; i < (fNTypes - 1); ++i) {
667 className.Form("Pair.%s", fgkPairClassNames[i]);
668 className2.Form("Track.Legs.%s", fgkPairClassNames[i]);
669 Bool_t pairClass = fHistos->HasHistClass(className);
670 Bool_t pairClass2 = (fHistoArray && fHistoArray->HasHistClass(className));
671 Bool_t legClass = fHistos->HasHistClass(className2);
672 Bool_t legClass2 = (fHistoArray && fHistoArray->HasHistClass(className2));
673
674 Bool_t legClassHits = kFALSE;
675
676 // check leg hits and mc signal filling
677 if (i == static_cast<Int_t>(EPairType::kSEPM)) {
678
679 // loop over all detectors and check for hit histos
680 for (ECbmModuleId idet = ECbmModuleId::kRef; idet < ECbmModuleId::kNofSystems; ++idet) {
681 className3 + Form("Hit.Legs.") + PairAnalysisHelper::GetDetName(idet);
682 legClassHits = fHistos->HasHistClass(className3);
683 if (legClassHits) break;
684 }
685
686
687 for (Int_t isig = 0; isig < nsig; isig++) {
688 sigName = Form("Pair_%s", fSignalsMC->At(isig)->GetName());
689 pairClassMC.SetBitNumber(isig, fHistos->HasHistClass(sigName));
690 pairClassMChf.SetBitNumber(isig, fHistoArray && fHistoArray->HasHistClass(sigName));
691 // Printf("fill %s: %d histos %d",sigName.Data(),pairClassMC.TestBitNumber(isig),fHistos->HasHistClass(sigName));
692 sigName = Form("Track.Legs_%s", fSignalsMC->At(isig)->GetName());
693 legClassMC.SetBitNumber(isig, fHistos->HasHistClass(sigName));
694 legClassMChf.SetBitNumber(isig, fHistoArray && fHistoArray->HasHistClass(sigName));
695
696 // check mc signal leg hits filling
697 if (!legClassHits) {
698 // loop over all detectors
699 for (ECbmModuleId idet = ECbmModuleId::kRef; idet < ECbmModuleId::kNofSystems; ++idet) {
700 className3 = Form("Hit.Legs.") + PairAnalysisHelper::GetDetName(idet);
701 sigName = className3 + "_" + fSignalsMC->At(isig)->GetName();
702 legClassHits = fHistos->HasHistClass(className3);
703 if (legClassHits) break; // abort when at least something should be filled
704 }
705 }
706 }
707 }
708
709 Bool_t fill = (pairClass || pairClass2 || legClass || legClass2 || legClassHits);
710 if (!fill && !legClassMC.CountBits() && !legClassMChf.CountBits() && !pairClassMC.CountBits()
711 && !pairClassMChf.CountBits())
712 continue;
713
714 // loop over all pairs
715 UInt_t selectedMask = (1 << fPairFilter.GetCuts()->GetEntries()) - 1; // for internal train wagons only
716 Int_t npairs = PairArray(i)->GetEntriesFast();
717 for (Int_t ipair = 0; ipair < npairs; ++ipair) {
718 PairAnalysisPair* pair = static_cast<PairAnalysisPair*>(PairArray(i)->UncheckedAt(ipair));
719
720 //fill pair information
721 if (pairClass || pairClass2 || pairClassMC.CountBits() || pairClassMChf.CountBits()) {
722
723 // if(i==3)
724 // printf("train: %d \t pair type: %s \t p:%p d1:%p d2:%p \n",!fEventProcess,fgkPairClassNames[i],pair,
725 // pair->GetFirstDaughter()->GetGlobalTrack(),pair->GetSecondDaughter()->GetGlobalTrack());
726 // in case of internal train wagon do the cut selections
727 if (!fEventProcess) {
728 UInt_t cutMask = fPairFilter.IsSelected(pair);
729 // apply cuts
730 if (cutMask != selectedMask) continue;
731 }
732
733 // fill histograms
734 PairAnalysisVarManager::Fill(pair, values);
735 if (pairClass) fHistos->FillClass(className, values);
736 if (pairClass2) fHistoArray->FillClass(className, values);
737
738 // check mc filling
739 fillMC.ResetAllBits();
740 if (i == static_cast<Int_t>(EPairType::kSEPM)) {
741 for (Int_t isig = 0; isig < nsig; isig++) {
742 // next if nothing needs to be filled
743 if (!pairClassMC.TestBitNumber(isig) && !pairClassMChf.TestBitNumber(isig)) continue;
744 sigMC = (PairAnalysisSignalMC*) fSignalsMC->At(isig);
745 // next if single particle signal
746 if (sigMC->IsSingleParticle()) continue;
747 Bool_t isMCtruth = mc->IsMCTruth(pair, sigMC);
748 fillMC.SetBitNumber(isig, isMCtruth);
749 if (!isMCtruth) continue;
750 sigName = Form("Pair_%s", sigMC->GetName());
751 if (pairClassMC.TestBitNumber(isig)) fHistos->FillClass(sigName, values);
752 if (pairClassMChf.TestBitNumber(isig)) fHistoArray->FillClass(sigName, values);
753 // Double_t wght = sigMC->GetWeight(values);
754 // if(wght!= values[PairAnalysisVarManager::kWeight])
755 // Warning("FillHistograms","pair weight from tracks (%.3e) not matching MC weight (%.3e) for signal %s",
756 // values[PairAnalysisVarManager::kWeight],wght,sigMC->GetName());
757 }
758 }
759 }
760
761 //fill leg information, don't fill the information twice
762 if (legClass || legClass2 || legClassMC.CountBits() || legClassMChf.CountBits() || legClassHits) {
763 PairAnalysisTrack* d1 = pair->GetFirstDaughter();
764 PairAnalysisTrack* d2 = pair->GetSecondDaughter();
765 if (!arrLegs.FindObject(d1)) {
767 if (legClass) fHistos->FillClass(className2, values);
768 if (legClass2) fHistoArray->FillClass(className2, values);
769 // mc signal filling
770 if (i == static_cast<Int_t>(EPairType::kSEPM)) {
771 for (Int_t isig = 0; isig < nsig; isig++) {
772 if (!fillMC.TestBitNumber(isig)) continue;
773 sigMC = (PairAnalysisSignalMC*) fSignalsMC->At(isig);
774 sigName = Form("Track.Legs_%s", sigMC->GetName());
775 if (legClassMC.TestBitNumber(isig)) fHistos->FillClass(sigName, values);
776 if (legClassMChf.TestBitNumber(isig)) fHistoArray->FillClass(sigName, values);
777 }
778 }
779
780 //Fill leg hits information
781 if (legClassHits) FillHistogramsHits(ev, &fillMC, d1, kTRUE, values);
782
783 arrLegs.Add(d1);
784 }
785 if (!arrLegs.FindObject(d2)) {
787 if (legClass) fHistos->FillClass(className2, values);
788 if (legClass2) fHistoArray->FillClass(className2, values);
789 // mc signal filling
790 if (i == static_cast<Int_t>(EPairType::kSEPM)) {
791 for (Int_t isig = 0; isig < nsig; isig++) {
792 if (!fillMC.TestBitNumber(isig)) continue;
793 sigMC = (PairAnalysisSignalMC*) fSignalsMC->At(isig);
794 sigName = Form("Track.Legs_%s", sigMC->GetName());
795 if (legClassMC.TestBitNumber(isig)) fHistos->FillClass(sigName, values);
796 if (legClassMChf.TestBitNumber(isig)) fHistoArray->FillClass(sigName, values);
797 }
798 }
799
800 //Fill leg hits information
801 if (legClassHits) FillHistogramsHits(ev, &fillMC, d2, kTRUE, values);
802
803 arrLegs.Add(d2);
804 }
805 }
806 }
807 if (legClass || legClass2) arrLegs.Clear();
808 }
809}
810//________________________________________________________________
811void PairAnalysis::FillHistogramsPair(PairAnalysisPair* pair, Bool_t fromPreFilter /*=kFALSE*/)
812{
813 //
814 // Fill Histogram information for pairs and the track in the pair
815 // NOTE: in this funtion the leg information may be filled multiple
816 // times. This funtion is used in the track rotation pairing
817 // and those legs are not saved!
818 //
819 TString className, className2;
820 Double_t* values = PairAnalysisVarManager::GetData();
822
823 //Fill Pair information, separately for all pair candidate arrays and the legs
824 TObjArray arrLegs(100);
825 const Int_t type = pair->GetType();
826 if (fromPreFilter) {
827 className.Form("RejPair.%s", fgkPairClassNames[type]);
828 className2.Form("RejTrack.%s", fgkPairClassNames[type]);
829 }
830 else {
831 className.Form("Pair.%s", fgkPairClassNames[type]);
832 className2.Form("Track.Legs.%s", fgkPairClassNames[type]);
833 }
834
835 Bool_t pairClass = fHistos->HasHistClass(className);
836 Bool_t pairClass2 = (fHistoArray && fHistoArray->HasHistClass(className));
837 Bool_t legClass = fHistos->HasHistClass(className2);
838 Bool_t legClass2 = (fHistoArray && fHistoArray->HasHistClass(className2));
839
840 //fill pair information
841 if (pairClass || pairClass2) {
842 PairAnalysisVarManager::Fill(pair, values);
843 if (pairClass) fHistos->FillClass(className, values);
844 if (pairClass2) fHistoArray->FillClass(className, values);
845 }
846
847 if (legClass || legClass2) {
848 PairAnalysisTrack* d1 = pair->GetFirstDaughter();
850 if (legClass) fHistos->FillClass(className2, values);
851 if (legClass2) fHistoArray->FillClass(className2, values);
852
853 PairAnalysisTrack* d2 = pair->GetSecondDaughter();
855 if (legClass) fHistos->FillClass(className2, values);
856 if (legClass2) fHistoArray->FillClass(className2, values);
857 }
858}
859
860//________________________________________________________________
862{
863 //
864 // select tracks and fill track candidate arrays
865 // use track arrays 0 and 1
866 //
867
869 PairAnalysisMC* papaMC = 0x0;
870 if (fHasMC && fSignalsMC) papaMC = PairAnalysisMC::Instance();
871
872 // get event data
873 Double_t* values = PairAnalysisVarManager::GetData();
874
875 Int_t ntracks = ev->GetNumberOfTracks();
876 UInt_t selectedMask = (1 << fTrackFilter.GetCuts()->GetEntries()) - 1;
877 for (Int_t itrack = 0; itrack < ntracks; ++itrack) {
878
879 //get particle
880 PairAnalysisTrack* particle = ev->GetTrack(itrack);
881
882 // adapt mass hypothesis accordingly (they were initialized with PDG11)
884
885 // fill variables
886 PairAnalysisVarManager::Fill(particle, values);
887
888 //apply track cuts
889 UInt_t cutmask = fTrackFilter.IsSelected(values);
890 //UInt_t cutmask=fTrackFilter.IsSelected(particle);
891 //fill cut QA
892 if (fCutQA) fQAmonitor->FillAll(particle);
893 if (fCutQA) fQAmonitor->Fill(cutmask, particle);
894
895 // if raw spectra before any cuts are requested then fill
896 if (fHistos && fHistos->HasHistClass("Track.noCuts")) {
898 // PairAnalysisVarManager::Fill(particle, values);
899 fHistos->FillClass("Track.noCuts", values);
900 }
901
903 if (fTrackFilter.GetHistogramList()->GetSize()) FillCutStepHistograms(&fTrackFilter, cutmask, particle, values);
904
905 // rejection
906 if (cutmask != selectedMask) continue;
907
908 // store signal weights in the tracks - ATTENTION later signals should be more specific
909 if (fHasMC && fSignalsMC) {
910 // printf("particle %p: pdg: %.0f \t mother: %.0f grand: %.0f \n", particle, values[PairAnalysisVarManager::kPdgCode],
911 // values[PairAnalysisVarManager::kPdgCodeMother], values[PairAnalysisVarManager::kPdgCodeGrandMother]);
912 for (Int_t isig = 0; isig < fSignalsMC->GetEntriesFast(); isig++) {
914 Double_t wght = sigMC->GetWeight(values);
915 Bool_t useMCweight = (TMath::Abs(wght - 1.) > 1.0e-15);
916 if (!useMCweight) continue;
917 // printf("\t temp. particle weight: %f \t signal weight for %s is %f \n",particle->GetWeight(),sigMC->GetName(),sigMC->GetWeight());
918 if (papaMC->IsMCTruth(particle, sigMC, 1) || papaMC->IsMCTruth(particle, sigMC, 2)) {
919 // printf("particle weight: %f \t signal weight for %s is %f \n",particle->GetWeight(),sigMC->GetName(),sigMC->GetWeight());
920 // if(sig->GetWeight(values) != 1.0) particle->SetWeight( sig->GetWeight(values) );
921 particle->SetWeight(wght);
922 // printf("no weight set for %s, use default (1.)",sigMC->GetName());
923 }
924 }
925 // printf(" -----> final particle weight: %f \n",particle->GetWeight());
926
927 Double_t wght = particle->GetWeight();
928 Bool_t hasMCweight = (TMath::Abs(wght - 1.) > 1.0e-15);
929
932 Int_t label = particle->GetLabel();
933 while (!hasMCweight) {
934
935 Int_t motherLabel = papaMC->GetMothersLabel(label);
936 if (motherLabel == -1) break; // stop if primary particle is reached
937
939 for (Int_t isig = 0; isig < fSignalsMC->GetEntriesFast(); isig++) {
941 Double_t wghtMC = sigMC->GetWeight(values);
942 Bool_t useMCweight = (TMath::Abs(wghtMC - 1.) > 1.0e-15);
943 if (!useMCweight) continue; // signal has no weights applied
944
945 if (papaMC->IsMCTruth(motherLabel, sigMC, 1) || papaMC->IsMCTruth(motherLabel, sigMC, 2)) {
946 particle->SetWeight(wghtMC);
947 hasMCweight = kTRUE;
948 }
949 }
950
952 label = motherLabel;
953 }
954
955 } //end: store mc weights
956
957 //fill selected particle into the corresponding track arrays
958 Short_t charge = particle->Charge();
959 if (charge > 0) fTracks[0].Add(particle); // positive tracks
960 else if (charge < 0)
961 fTracks[1].Add(particle); // negative tracks
962 }
963}
964
965//________________________________________________________________
966void PairAnalysis::PairPreFilter(Int_t arr1, Int_t arr2, TObjArray& arrTracks1, TObjArray& arrTracks2)
967{
968 //
969 // Prefilter tracks from pairs
970 // Needed for datlitz rejections
971 // remove all tracks from the Single track arrays that pass the cuts in this filter
972 //
973
980 Double_t mass = TDatabasePDG::Instance()->GetParticle(fPdgLeg1)->Mass();
981 Double_t pt = gRandom->Exp(3.); // return t from exp(-t/3.)
982 if (pt < 0.075) pt = 0.075;
983 Double_t eta = -TMath::Log(TMath::Tan((gRandom->Uniform(2.5, 25.) / 180. * TMath::Pi()) / 2));
984 Double_t phi = gRandom->Uniform(TMath::TwoPi());
986 t1->SetPdgCode(fPdgLeg1);
987 t1->GetMomentum()->SetPtEtaPhiM(pt, eta, phi, mass);
988 fTracks[0].Add(t1); // positive tracks
990 mass = TDatabasePDG::Instance()->GetParticle(fPdgLeg2)->Mass();
991 pt = gRandom->Exp(3.); // return t from exp(-t/3.)
992 if (pt < 0.075) pt = 0.075;
993 eta = -TMath::Log(TMath::Tan((gRandom->Uniform(2.5, 25.) / 180. * TMath::Pi()) / 2));
994 phi = gRandom->Uniform(TMath::TwoPi());
996 t2->SetPdgCode(fPdgLeg2);
997 t2->GetMomentum()->SetPtEtaPhiM(pt, eta, phi, mass);
998 fTracks[1].Add(t2); // negative tracks
1000
1004 UInt_t selectedMaskLeg = (1 << fPairPreFilterLegs.GetCuts()->GetEntries()) - 1;
1005 Bool_t isLeg1selected = kTRUE;
1006 Bool_t isLeg2selected = kTRUE;
1007
1008
1009 Int_t ntrack1 = arrTracks1.GetEntriesFast();
1010 Int_t ntrack2 = arrTracks2.GetEntriesFast();
1011
1013 Bool_t* bTracks1 = new Bool_t[ntrack1];
1014 for (Int_t itrack1 = 0; itrack1 < ntrack1; ++itrack1)
1015 bTracks1[itrack1] = kFALSE;
1016 Bool_t* bTracks2 = new Bool_t[ntrack2];
1017 for (Int_t itrack2 = 0; itrack2 < ntrack2; ++itrack2)
1018 bTracks2[itrack2] = kFALSE;
1019
1021 bTracks1[ntrack1 - 1] = kTRUE;
1022 bTracks2[ntrack2 - 1] = kTRUE;
1023
1024 // candiate
1025 PairAnalysisPair* candidate;
1026 if (fUseKF) candidate = new PairAnalysisPairKF();
1027 else
1028 candidate = new PairAnalysisPairLV();
1029 candidate->SetKFUsage(fUseKF);
1030
1031 UInt_t selectedMask = (1 << fPairPreFilter.GetCuts()->GetEntries()) - 1;
1032
1033 Int_t nRejPasses = 1; //for fPreFilterUnlikeOnly and no set flag
1034 if (fPreFilterAllSigns) nRejPasses = 3;
1035
1036 // loop over rejection passes: OS (+ 2xLS)
1037 for (Int_t iRP = 0; iRP < nRejPasses; ++iRP) {
1038
1039 // default rejection pass OS
1040 Int_t arr1RP = arr1, arr2RP = arr2;
1041 TObjArray* arrTracks1RP = &arrTracks1;
1042 TObjArray* arrTracks2RP = &arrTracks2;
1043 Bool_t* bTracks1RP = bTracks1;
1044 Bool_t* bTracks2RP = bTracks2;
1045
1046 // change for LS rejection passes
1047 switch (iRP) {
1048 case 1:
1049 arr1RP = arr1;
1050 arr2RP = arr1;
1051 arrTracks1RP = &arrTracks1;
1052 arrTracks2RP = &arrTracks1;
1053 bTracks1RP = bTracks1;
1054 bTracks2RP = bTracks1;
1055 break;
1056 case 2:
1057 arr1RP = arr2;
1058 arr2RP = arr2;
1059 arrTracks1RP = &arrTracks2;
1060 arrTracks2RP = &arrTracks2;
1061 bTracks1RP = bTracks2;
1062 bTracks2RP = bTracks2;
1063 break;
1064 default:; //nothing to do
1065 }
1066
1067 Int_t ntrack1RP = (*arrTracks1RP).GetEntriesFast();
1068 Int_t ntrack2RP = (*arrTracks2RP).GetEntriesFast();
1069
1070 Int_t pairIndex = GetPairIndex(arr1RP, arr2RP);
1071 // loop over all tracks in both arrays
1072 for (Int_t itrack1 = 0; itrack1 < ntrack1RP; ++itrack1) {
1073 Int_t end = ntrack2RP;
1074 if (arr1RP == arr2RP) end = itrack1;
1075
1076 TObject* track1 = (*arrTracks1RP).UncheckedAt(itrack1);
1077 if (!track1) continue;
1078
1080 // pre-pairing filter for the first track candidate in the pairing step
1081 // kBothLegs demands that both track candidates fullfill the cuts
1082 // kAnyLeg demands that at least one candidate fullfills the cuts
1083 // kOneLeg demands exactly one candidate to fullfill the cuts
1084 if (selectedMaskLeg) {
1085 isLeg1selected = (fPairPreFilterLegs.IsSelected(static_cast<PairAnalysisTrack*>(track1))
1086 == selectedMaskLeg); // check cuts for the first track candidate
1087 if (fCutType == ECutType::kBothLegs && !isLeg1selected) continue;
1088 }
1089
1090
1091 for (Int_t itrack2 = 0; itrack2 < end; ++itrack2) {
1092 TObject* track2 = (*arrTracks2RP).UncheckedAt(itrack2);
1093 if (!track2) continue;
1094
1096 if (selectedMaskLeg) {
1097 if (fCutType != ECutType::kAnyLeg || (fCutType == ECutType::kAnyLeg && !isLeg1selected)) {
1098 isLeg2selected = (fPairPreFilterLegs.IsSelected(static_cast<PairAnalysisTrack*>(track2))
1099 == selectedMaskLeg); // check cuts for the second track candidate
1100 if (fCutType == ECutType::kBothLegs && !isLeg2selected) continue; // the second track did not fullfill cuts
1101 if (fCutType == ECutType::kAnyLeg && !isLeg2selected)
1102 continue; // the first and second track did not fullfill cuts
1103 if (fCutType == ECutType::kOneLeg && isLeg1selected == isLeg2selected)
1104 continue; // both tracks have the same value
1105 }
1106 }
1107
1109 candidate->SetTracks(static_cast<PairAnalysisTrack*>(track1), fPdgLeg1, static_cast<PairAnalysisTrack*>(track2),
1110 fPdgLeg2);
1111
1112 candidate->SetType(pairIndex);
1113 candidate->SetLabel(PairAnalysisMC::Instance()->GetLabelMotherWithPdg(candidate, fPdgMother));
1114
1116 Bool_t testParticle = (track1 == t1 || track1 == t2 || track2 == t1 || track2 == t2);
1117
1119 UInt_t cutmask = fPairPreFilter.IsSelected(candidate);
1120
1122 if (!testParticle) {
1123 if (fCutQA) fQAmonitor->FillAll(candidate, 1);
1124 if (fCutQA) fQAmonitor->Fill(cutmask, candidate, 1);
1125 }
1126
1128 if (cutmask != selectedMask) { continue; }
1129
1131 if (testParticle) {
1132 // set variable to randomrejection probability to 1
1134 continue;
1135 }
1136
1138 if (fHistos) FillHistogramsPair(candidate, kTRUE); // kTRUE: fromPrefilter
1139
1141 bTracks1RP[itrack1] = kTRUE;
1142 bTracks2RP[itrack2] = kTRUE;
1143 }
1144 }
1145
1146 } // end rejection passes
1147
1148 //clear surplus candiate
1149 delete candidate;
1150
1151 //remove the tracks from the Track arrays
1152 for (Int_t itrack1 = 0; itrack1 < ntrack1; ++itrack1) {
1153 if (bTracks1[itrack1]) arrTracks1.AddAt(0x0, itrack1);
1154 }
1155 for (Int_t itrack2 = 0; itrack2 < ntrack2; ++itrack2) {
1156 if (bTracks2[itrack2]) arrTracks2.AddAt(0x0, itrack2);
1157 }
1158
1159 // clean up
1160 delete[] bTracks1;
1161 delete[] bTracks2;
1162
1163 //compress the track arrays
1164 arrTracks1.Compress();
1165 arrTracks2.Compress();
1166
1167 /*
1168 //apply leg cuts after the pre filter
1169 if ( fFinalTrackFilter.GetCuts()->GetEntries()>0 ) {
1170 selectedMask=(1<<fFinalTrackFilter.GetCuts()->GetEntries())-1;
1171 //loop over tracks from array 1
1172 for (Int_t itrack=0; itrack<arrTracks1.GetEntriesFast();++itrack){
1173 //test cuts
1174 UInt_t cutMask=fFinalTrackFilter.IsSelected(arrTracks1.UncheckedAt(itrack));
1175
1176 //apply cut
1177 if (cutMask!=selectedMask) arrTracks1.AddAt(0x0,itrack);
1178 }
1179 arrTracks1.Compress();
1180
1181 //in case of like sign don't loop over second array
1182 if (arr1==arr2) {
1183 arrTracks2=arrTracks1;
1184 } else {
1185
1186 //loop over tracks from array 2
1187 for (Int_t itrack=0; itrack<arrTracks2.GetEntriesFast();++itrack){
1188 //test cuts
1189 UInt_t cutMask=fFinalTrackFilter.IsSelected(arrTracks2.UncheckedAt(itrack));
1190 //apply cut
1191 if (cutMask!=selectedMask) arrTracks2.AddAt(0x0,itrack);
1192 }
1193 arrTracks2.Compress();
1194
1195 }
1196 }
1197 //For unlike-sign monitor track-cuts:
1198 if (arr1!=arr2&&fHistos) {
1199 TObjArray *unlikesignArray[2] = {&arrTracks1,&arrTracks2};
1200 FillHistogramsTracks(unlikesignArray);
1201 }
1202 */
1203}
1204
1205//________________________________________________________________
1206void PairAnalysis::FilterTrackArrays(TObjArray& arrTracks1, TObjArray& arrTracks2)
1207{
1208 //
1209 // select tracks and adapt track candidate arrays
1210 // second and final track selection
1211 //
1212
1213 // get event data
1214 Double_t* values = PairAnalysisVarManager::GetData();
1215
1216 //apply leg cuts after the pre filter
1217 if (fFinalTrackFilter.GetCuts()->GetEntries() < 1) return;
1218
1219 UInt_t selectedMask = (1 << fFinalTrackFilter.GetCuts()->GetEntries()) - 1;
1220 //loop over tracks from array 1
1221 for (Int_t itrack = 0; itrack < arrTracks1.GetEntriesFast(); ++itrack) {
1222
1223 //get particle
1224 PairAnalysisTrack* particle = static_cast<PairAnalysisTrack*>(arrTracks1.UncheckedAt(itrack));
1225
1226 // fill variables
1227 PairAnalysisVarManager::Fill(particle, values);
1228
1229 //apply cuts
1230 UInt_t cutmask = fFinalTrackFilter.IsSelected(values);
1231 // UInt_t cutmask=fFinalTrackFilter.IsSelected(particle);
1232 //fill cut QA
1233 if (fCutQA) fQAmonitor->FillAll(particle, 1);
1234 if (fCutQA) fQAmonitor->Fill(cutmask, particle, 1);
1235
1237 if (fFinalTrackFilter.GetHistogramList()->GetSize())
1238 FillCutStepHistograms(&fFinalTrackFilter, cutmask, particle, values);
1239
1240 // rejection
1241 if (cutmask != selectedMask) arrTracks1.AddAt(0x0, itrack);
1242 }
1243 arrTracks1.Compress();
1244
1245 //loop over tracks from array 2
1246 for (Int_t itrack = 0; itrack < arrTracks2.GetEntriesFast(); ++itrack) {
1247
1248 //get particle
1249 PairAnalysisTrack* particle = static_cast<PairAnalysisTrack*>(arrTracks2.UncheckedAt(itrack));
1250
1251 // fill variables
1252 PairAnalysisVarManager::Fill(particle, values);
1253
1254 //apply cuts
1255 UInt_t cutmask = fFinalTrackFilter.IsSelected(values);
1256 // UInt_t cutmask=fFinalTrackFilter.IsSelected(particle);
1257
1258 //fill cut QA
1259 if (fCutQA) fQAmonitor->FillAll(particle, 1);
1260 if (fCutQA) fQAmonitor->Fill(cutmask, particle, 1);
1261
1263 if (fFinalTrackFilter.GetHistogramList()->GetSize())
1264 FillCutStepHistograms(&fFinalTrackFilter, cutmask, particle, values);
1265
1266 // rejection
1267 if (cutmask != selectedMask) arrTracks2.AddAt(0x0, itrack);
1268 }
1269 arrTracks2.Compress();
1270}
1271
1272//________________________________________________________________
1273void PairAnalysis::FillPairArrays(Int_t arr1, Int_t arr2)
1274{
1275 //
1276 // select pairs and fill pair candidate arrays
1277 //
1278
1279 TObjArray arrTracks1 = fTracks[arr1];
1280 TObjArray arrTracks2 = fTracks[arr2];
1281
1282 //process pre filter if set
1283 // if ((!fPreFilterAllSigns) && (!fPreFilterUnlikeOnly) && ( fPairPreFilter.GetCuts()->GetEntries()>0 )) PairPreFilter(arr1, arr2, arrTracks1, arrTracks2);
1284
1285 Int_t pairIndex = GetPairIndex(arr1, arr2);
1286
1287 Int_t ntrack1 = arrTracks1.GetEntriesFast();
1288 Int_t ntrack2 = arrTracks2.GetEntriesFast();
1289
1294 PairAnalysisPair* candidate;
1295 if (fUseKF && pairIndex <= static_cast<Int_t>(EPairType::kSEMM)) candidate = new PairAnalysisPairKF();
1296 else
1297 candidate = new PairAnalysisPairLV();
1298 candidate->SetKFUsage(fUseKF);
1299 candidate->SetType(pairIndex);
1300
1301 UInt_t selectedMask = (1 << fPairFilter.GetCuts()->GetEntries()) - 1;
1302
1303 for (Int_t itrack1 = 0; itrack1 < ntrack1; ++itrack1) {
1304 Int_t end = ntrack2;
1305 if (arr1 == arr2) end = itrack1;
1306 for (Int_t itrack2 = 0; itrack2 < end; ++itrack2) {
1307 //create the pair (direct pointer to the memory by this daughter reference are kept also for ME)
1308 candidate->SetTracks(&(*static_cast<PairAnalysisTrack*>(arrTracks1.UncheckedAt(itrack1))), fPdgLeg1,
1309 &(*static_cast<PairAnalysisTrack*>(arrTracks2.UncheckedAt(itrack2))), fPdgLeg2);
1310 // TODO: maybe set here the mother pdg code and remove fPdgMother
1311 Int_t label = PairAnalysisMC::Instance()->GetLabelMotherWithPdg(candidate, fPdgMother);
1312 candidate->SetLabel(label);
1313 if (label > -1) candidate->SetPdgCode(fPdgMother);
1314 else
1315 candidate->SetPdgCode(0);
1316
1317 //pair cuts
1318 UInt_t cutMask = fPairFilter.IsSelected(candidate);
1319
1320 // cut qa
1321 if (pairIndex == static_cast<Int_t>(EPairType::kSEPM) && fCutQA) {
1322 fQAmonitor->FillAll(candidate);
1323 fQAmonitor->Fill(cutMask, candidate);
1324 }
1325
1326 //apply cut
1327 if (cutMask != selectedMask) continue;
1328
1329 // if(pairIndex==3)
1330 // printf("fill pair array \t train: %d \t pair type: %s \t p:%p d1:%p d2:%p \n",!fEventProcess,fgkPairClassNames[pairIndex],candidate,
1331 // candidate->GetFirstDaughter()->GetGlobalTrack(),candidate->GetSecondDaughter()->GetGlobalTrack());
1332 //histogram array for the pair
1334
1335 //add the candidate to the candidate array
1336 PairArray(pairIndex)->Add(candidate);
1337 //get a new candidate
1338 if (fUseKF && pairIndex <= static_cast<Int_t>(EPairType::kSEMM)) candidate = new PairAnalysisPairKF();
1339 else
1340 candidate = new PairAnalysisPairLV();
1341 candidate->SetKFUsage(fUseKF);
1342 candidate->SetType(pairIndex);
1343 }
1344 }
1345 //delete the surplus candidate
1346 delete candidate;
1347}
1348
1349//________________________________________________________________
1351{
1352 //
1353 // rotate pairs and fill into pair candidate arrays
1354 //
1355
1356 Int_t ntrack1 = fTracks[0].GetEntriesFast();
1357 Int_t ntrack2 = fTracks[1].GetEntriesFast();
1358
1359 PairAnalysisPair* candidate;
1360 if (fUseKF) candidate = new PairAnalysisPairKF();
1361 else
1362 candidate = new PairAnalysisPairLV();
1363 candidate->SetKFUsage(fUseKF);
1364 candidate->SetType(static_cast<Int_t>(EPairType::kSEPMRot));
1365
1366 UInt_t selectedMask = (1 << fPairFilter.GetCuts()->GetEntries()) - 1;
1367 // loop over track arrays
1368 for (Int_t itrack1 = 0; itrack1 < ntrack1; ++itrack1) {
1369 for (Int_t itrack2 = 0; itrack2 < ntrack2; ++itrack2) {
1370
1371 // loop over iterations
1372 for (Int_t irot = 0; irot < fTrackRotator->GetIterations(); ++irot) {
1373 // build candidate
1374 candidate->SetTracks(&(*static_cast<PairAnalysisTrack*>(fTracks[0].UncheckedAt(itrack1))), fPdgLeg1,
1375 &(*static_cast<PairAnalysisTrack*>(fTracks[1].UncheckedAt(itrack2))), fPdgLeg2);
1376 // rotate the candidates daughter track
1377 candidate->RotateTrack(fTrackRotator);
1378
1379 //pair cuts
1380 UInt_t cutMask = fPairFilter.IsSelected(candidate);
1381
1382 //apply cut
1383 if (cutMask != selectedMask) continue;
1384
1385 //histogram array for the pair
1386 if (fHistoArray) fHistoArray->Fill(static_cast<Int_t>(EPairType::kSEPMRot), candidate);
1387
1388 if (fHistos) FillHistogramsPair(candidate);
1389 if (fStoreRotatedPairs) {
1390 if (fUseKF)
1391 PairArray(static_cast<Int_t>(EPairType::kSEPMRot))
1392 ->Add(static_cast<PairAnalysisPairKF*>(candidate->Clone()));
1393 else
1394 PairArray(static_cast<Int_t>(EPairType::kSEPMRot))
1395 ->Add(static_cast<PairAnalysisPairLV*>(candidate->Clone()));
1396 // if(fUseKF) PairArray(kSEPMRot)->Add(new PairAnalysisPairKF(*candidate);
1397 // else PairArray(kSEPMRot)->Add(new PairAnalysisPairLV(*candidate);
1398 }
1399
1400 } // end of iterations
1401 } //arr0
1402 } //arr1
1403
1404 //delete the surplus candidate
1405 delete candidate;
1406}
1407
1408//________________________________________________________________
1410{
1411 //
1412 // Add an MC signal to the signals list
1413 //
1414 if (!fSignalsMC) {
1415 fSignalsMC = new TObjArray();
1416 fSignalsMC->SetOwner();
1417 }
1418 // sort mc signal (first single particle, then pair signals)
1419 if (signal->IsSingleParticle()) fSignalsMC->AddAtFree(signal);
1420 else
1421 fSignalsMC->AddAtAndExpand(signal, fSignalsMC->GetLast() < 10 ? 10 : fSignalsMC->GetLast() + 1);
1422 //fSignalsMC->Add(signal);
1423}
1424
1425//________________________________________________________________
1426Bool_t PairAnalysis::FillMCHistograms(Int_t label1, Int_t label2, Int_t nSignal)
1427{
1428 //
1429 // fill QA MC TRUTH histograms for pairs and legs of all added mc signals
1430 //
1431 PairAnalysisSignalMC* sigMC = (PairAnalysisSignalMC*) fSignalsMC->At(nSignal);
1432
1433 TString className = Form("Pair_%s_MCtruth", sigMC->GetName());
1434 TString className2 = Form("Track.Legs_%s_MCtruth", sigMC->GetName());
1435 TString className3 = Form("Track.%s_%s_MCtruth", fgkPairClassNames[1], sigMC->GetName());
1436 Bool_t pairClass = fHistos->HasHistClass(className.Data());
1437 Bool_t legClass = fHistos->HasHistClass(className2.Data());
1438 Bool_t trkClass = fHistos->HasHistClass(className3.Data());
1439 // printf("fill signal %d: pair %d legs %d trk %d \n",nSignal,pairClass,legClass,trkClass);
1440 if (!pairClass && !legClass && !trkClass) return kFALSE;
1441
1443 CbmMCTrack* part1 = papaMC->GetMCTrackFromMCEvent(label1);
1444 CbmMCTrack* part2 = papaMC->GetMCTrackFromMCEvent(label2);
1445 if (!part1 && !part2) return kFALSE;
1446 if (part1 && part2) {
1447 // fill only unlike sign (and only SE)
1448 if (part1->GetCharge() * part2->GetCharge() > 0) return kFALSE;
1449 }
1450
1451 Int_t mLabel1 = papaMC->GetMothersLabel(label1);
1452 Int_t mLabel2 = papaMC->GetMothersLabel(label2);
1453 // printf("leg/mother labels: %d/%d %d/%d \t part %p,%p \n",label1,mLabel1,label2,mLabel2,part1,part2);
1454
1455 // check the same mother option
1456 if (sigMC->GetMothersRelation() == PairAnalysisSignalMC::EBranchRelation::kSame && mLabel1 != mLabel2) return kFALSE;
1457 if (sigMC->GetMothersRelation() == PairAnalysisSignalMC::EBranchRelation::kDifferent && mLabel1 == mLabel2)
1458 return kFALSE;
1459
1460 // fill event values
1461 Double_t* values = PairAnalysisVarManager::GetData(); //NEW CHANGED
1464 // fill the leg variables
1465 if (legClass || trkClass) {
1466 if (part1) PairAnalysisVarManager::Fill(part1, values);
1468 if (part1 && trkClass) fHistos->FillClass(className3, values);
1469 if (part1 && part2 && legClass) fHistos->FillClass(className2, values);
1470 if (part2) PairAnalysisVarManager::Fill(part2, values);
1472 if (part2 && trkClass) fHistos->FillClass(className3, values);
1473 if (part1 && part2 && legClass) fHistos->FillClass(className2, values);
1474 }
1475
1477 FairMCPoint* pnt = NULL;
1478 TString className4;
1479 // loop over daughters
1480 for (Int_t ipart = 0; ipart < 2; ipart++) {
1481
1482 CbmMCTrack* part = (!ipart ? part1 : part2);
1483 Int_t label = (!ipart ? label1 : label2);
1484 if (!part) continue;
1485 if (ipart && sigMC->IsSingleParticle()) continue;
1486
1487 // leg and no leg loop
1488 for (Int_t ileg = 0; ileg < 2; ileg++) {
1489
1490 // loop over all detectors
1491 for (ECbmModuleId idet = ECbmModuleId::kRef; idet < ECbmModuleId::kNofSystems; ++idet) {
1492 className4 = Form("Hit.%s", (ileg ? "Legs." : ""));
1493 className4 += PairAnalysisHelper::GetDetName(idet) + "_" + sigMC->GetName() + "_MCtruth";
1494 if (!fHistos->HasHistClass(className4)) continue;
1495
1496 Int_t npnts = part->GetNPoints(idet);
1497 //printf("track %p(%d) \t has %d %s mc points \n",part,label1,npnts,PairAnalysisHelper::GetDetName(idet).Data());
1498 if (!npnts) continue;
1499
1500 TClonesArray* points = PairAnalysisVarManager::GetCurrentEvent()->GetPoints(idet); // get point array
1501 Int_t psize = points->GetSize();
1502 if (!points || psize < 1) continue;
1503
1504 Int_t nfnd = 0;
1505 for (Int_t idx = 0; idx < psize; idx++) {
1506 if (nfnd == npnts) break; // all points found
1507
1508 pnt = static_cast<FairMCPoint*>(points->At(idx));
1509 if (pnt->GetTrackID() == label) {
1510 // printf("det %s \t point index: %d/%d found! \n",PairAnalysisHelper::GetDetName(idet).Data(),idx,psize);
1511 nfnd++; // found point
1512 PairAnalysisVarManager::Fill(pnt, values);
1513 fHistos->FillClass(className4, values);
1514 }
1515 }
1516
1517 } //idet loop
1518 } // leg or not loop
1519 } // daughters loop
1520
1521 //fill pair information
1522 if (pairClass && part1 && part2) {
1523 PairAnalysisVarManager::FillVarMCParticle(part1, part2, values);
1525 fHistos->FillClass(className, values);
1526 }
1527
1528 return kTRUE;
1529}
1530
1531//________________________________________________________________
1532void PairAnalysis::FillHistogramsFromPairArray(Bool_t pairInfoOnly /*=kFALSE*/)
1533{
1534 //
1535 // Fill Histogram information for tracks and pairs
1536 //
1537
1538 TString className, className2;
1539 Double_t* values = PairAnalysisVarManager::GetData();
1541
1542 //Fill event information
1543 if (!pairInfoOnly) {
1544 if (fHistos->GetHistogramList()->FindObject("Event")) {
1545 fHistos->FillClass("Event", PairAnalysisVarManager::GetData());
1546 }
1547 }
1548
1549 UInt_t selectedMask = (1 << fPairFilter.GetCuts()->GetEntries()) - 1;
1550
1551 //Fill Pair information, separately for all pair candidate arrays and the legs
1552 TObjArray arrLegs(100);
1553 for (Int_t i = 0; i < (fNTypes - 1); ++i) {
1554 Int_t npairs = PairArray(i)->GetEntriesFast();
1555 if (npairs < 1) continue;
1556
1557 className.Form("Pair.%s", fgkPairClassNames[i]);
1558 className2.Form("Track.Legs_%s", fgkPairClassNames[i]);
1559 Bool_t pairClass = fHistos->HasHistClass(className);
1560 Bool_t legClass = fHistos->HasHistClass(className2);
1561
1562 // if (!pairClass&&!legClass) continue;
1563 for (Int_t ipair = 0; ipair < npairs; ++ipair) {
1564 PairAnalysisPair* pair = static_cast<PairAnalysisPair*>(PairArray(i)->UncheckedAt(ipair));
1565
1566 // apply cuts
1567 UInt_t cutMask = fPairFilter.IsSelected(pair);
1568
1569 // cut qa
1570 if (i == static_cast<Int_t>(EPairType::kSEPM) && fCutQA) {
1571 fQAmonitor->FillAll(pair);
1572 fQAmonitor->Fill(cutMask, pair);
1573 }
1574
1575 //apply cut
1576 if (cutMask != selectedMask) continue;
1577
1578 //histogram array for the pair
1579 if (fHistoArray) fHistoArray->Fill(i, pair);
1580
1581 // fill map
1583
1584 //fill pair information
1585 if (pairClass) {
1586 PairAnalysisVarManager::Fill(pair, values);
1587 fHistos->FillClass(className, values);
1588 }
1589
1590 //fill leg information, don't fill the information twice
1591 if (legClass) {
1592 PairAnalysisTrack* d1 = pair->GetFirstDaughter();
1593 PairAnalysisTrack* d2 = pair->GetSecondDaughter();
1594 if (!arrLegs.FindObject(d1)) {
1595 PairAnalysisVarManager::Fill(d1, values);
1596 fHistos->FillClass(className2, values);
1597 arrLegs.Add(d1);
1598 }
1599 if (!arrLegs.FindObject(d2)) {
1600 PairAnalysisVarManager::Fill(d2, values);
1601 fHistos->FillClass(className2, values);
1602 arrLegs.Add(d2);
1603 }
1604 }
1605 }
1606 if (legClass) arrLegs.Clear();
1607 }
1608}
1609
1610//________________________________________________________________
1612 const Double_t* values)
1613{
1614
1616 PairAnalysisMC* papaMC = 0x0;
1617 if (fHasMC && fSignalsMC) papaMC = PairAnalysisMC::Instance();
1618
1620 TString className;
1621 TString classNameMC;
1622 TString classNamePM = Form("Track.%s", fgkPairClassNames[1]);
1623
1624 PairAnalysisHistos histo;
1625 AnalysisCuts* cuts;
1626 TIter next(filter->GetCuts());
1627 Int_t iCut = 0;
1628
1630 Int_t nsig = (fSignalsMC ? fSignalsMC->GetEntriesFast() : 1);
1631 PairAnalysisSignalMC* sigMC;
1632 for (Int_t isig = 0; isig < nsig; isig++) {
1633 if (fSignalsMC) {
1634 sigMC = (PairAnalysisSignalMC*) fSignalsMC->At(isig);
1635 classNameMC = classNamePM + "_" + sigMC->GetName();
1636 // printf("fill cut details for %s \n",classNameMC.Data());
1637 }
1638 // check if machtes mc signal
1639 Bool_t isMCtruth = fSignalsMC && (papaMC->IsMCTruth(trk, sigMC, 1) || papaMC->IsMCTruth(trk, sigMC, 2));
1640 if (isig && !isMCtruth) continue;
1641
1643 // if(sigMC->GetWeight(values) != 1.0) trk->SetWeight( sigMC->GetWeight(values) );
1644
1646 next.Reset();
1647 iCut = 0;
1648
1650 while ((cuts = (AnalysisCuts*) next())) {
1652 UInt_t cutRef = (1 << (iCut + 1)) - 1; // increasing cut match
1653 // printf(" fill cut %s for track %p in hist details \n",cuts->GetName(),trk);
1654
1656 if ((cutmask & cutRef) == cutRef) {
1657
1658 // printf(" track %p passed cut \n",trk);
1660 histo.SetHistogramList(*(THashList*) filter->GetHistogramList()->FindObject(cuts->GetName()), kFALSE);
1661
1663 if (!isig) {
1664 histo.FillClass(classNamePM, values);
1665 for (Int_t i = 0; i < fLegTypes; ++i) {
1666 className.Form("Track.%s", fgkTrackClassNames[i]);
1667 histo.FillClass(className, values);
1668 }
1669 }
1671 if (isMCtruth) histo.FillClass(classNameMC, values);
1672 }
1673 iCut++;
1674 }
1675 }
1676}
1677
1678//________________________________________________________________
1679void PairAnalysis::FillCutStepHistogramsMC(AnalysisFilter* filter, UInt_t cutmask, Int_t label, const Double_t* values)
1680{
1681
1683 PairAnalysisMC* papaMC = 0x0;
1684 if (fHasMC && fSignalsMC) papaMC = PairAnalysisMC::Instance();
1685
1687 TString className;
1688 TString classNameMC;
1689 TString classNamePM = Form("Track.%s", fgkPairClassNames[1]);
1690
1691 PairAnalysisHistos histo;
1692 AnalysisCuts* cuts;
1693 TIter next(filter->GetCuts());
1694 Int_t iCut = 0;
1695
1697 Int_t nsig = (fSignalsMC ? fSignalsMC->GetEntriesFast() : 0);
1698 PairAnalysisSignalMC* sigMC;
1699 for (Int_t isig = 0; isig < nsig; isig++) {
1700
1701 sigMC = (PairAnalysisSignalMC*) fSignalsMC->At(isig);
1702 classNameMC = classNamePM + "_" + sigMC->GetName() + "_MCtruth";
1703 // printf("fill cut details for %s \n",classNameMC.Data());
1704
1705 // Proceed only if this signal is required in the pure MC step
1706 if (!sigMC->GetFillPureMCStep()) continue;
1707
1708 // check if matches mc signal
1709 Bool_t isMCtruth = fSignalsMC && (papaMC->IsMCTruth(label, sigMC, 1) || papaMC->IsMCTruth(label, sigMC, 2));
1710 if (isig && !isMCtruth) continue;
1711
1713 // if(sigMC->GetWeight(values) != 1.0) trk->SetWeight( sigMC->GetWeight(values) );
1714
1716 next.Reset();
1717 iCut = 0;
1718
1720 while ((cuts = (AnalysisCuts*) next())) {
1722 UInt_t cutRef = (1 << (iCut + 1)) - 1; // increasing cut match
1723
1725 if ((cutmask & cutRef) == cutRef) {
1726
1728 histo.SetHistogramList(*(THashList*) filter->GetHistogramList()->FindObject(cuts->GetName()), kFALSE);
1729
1731 if (isMCtruth) histo.FillClass(classNameMC, values);
1732 }
1733 iCut++;
1734 }
1735 }
1736}
1737
1738//________________________________________________________________
1740 Bool_t trackIsLeg, Double_t* values)
1741{
1742 //
1743 // Fill Histogram information for hits and hits of legs
1744 //
1745 TString className;
1746
1747 Int_t nsig = (fSignalsMC ? fSignalsMC->GetEntriesFast() : 0);
1748 TBits hitClassMC(nsig);
1749 TBits hitClassMChf(nsig);
1750
1751 TString sigName;
1752 PairAnalysisMC* mc = (nsig ? PairAnalysisMC::Instance() : 0x0);
1753
1754
1755 // loop over all detectors
1756 for (ECbmModuleId idet = ECbmModuleId::kRef; idet < ECbmModuleId::kNofSystems; ++idet) {
1757
1758 // detectors implemented
1759 switch (idet) {
1760 case ECbmModuleId::kMvd:
1761 case ECbmModuleId::kSts:
1763 case ECbmModuleId::kTrd:
1765 case ECbmModuleId::kTof: /* */ break;
1766 default: continue;
1767 }
1768
1770 className.Form("Hit.%s", (trackIsLeg ? "Legs." : ""));
1771 className += PairAnalysisHelper::GetDetName(idet); // detector hit
1772
1773 Bool_t hitClass = fHistos->HasHistClass(className);
1774 Bool_t hitClass2 = (fHistoArray && fHistoArray->HasHistClass(className));
1775
1776 // check mc signal filling
1777 for (Int_t isig = 0; isig < nsig; isig++) {
1778 sigName = className + "_" + fSignalsMC->At(isig)->GetName();
1779 hitClassMC.SetBitNumber(isig, fHistos->HasHistClass(sigName));
1780 hitClassMChf.SetBitNumber(isig, fHistoArray && fHistoArray->HasHistClass(sigName));
1781 }
1782 if (!hitClass && !hitClass2 && !hitClassMC.CountBits() && !hitClassMChf.CountBits()) continue;
1783
1784 // get hit array
1785 TClonesArray* hits = ev->GetHits(idet);
1786 if (!hits || hits->GetSize() < 1) continue;
1787
1788 // get matched track and mc track index/id
1789 CbmTrackMatchNew* tmtch = track->GetTrackMatch(idet);
1790 Int_t mctrk = (tmtch ? tmtch->GetMatchedLink().GetIndex() : -1);
1791 // Printf("mc track id via track match (%p) link: %d",tmtch,mctrk);
1792
1793 // get detector tracks
1794 CbmTrack* trkl = 0x0;
1795 CbmRichRing* ring = 0x0;
1796 switch (idet) {
1797 case ECbmModuleId::kMvd:
1798 case ECbmModuleId::kSts:
1800 case ECbmModuleId::kTrd: trkl = track->GetTrack(idet); break;
1801 case ECbmModuleId::kRich: ring = track->GetRichRing(); break;
1802 case ECbmModuleId::kTof: /* */ break;
1803 default: continue;
1804 }
1805
1806 // get number of hits
1807 Int_t nhits = 0;
1808 switch (idet) {
1809 case ECbmModuleId::kMvd:
1810 if (trkl) nhits = static_cast<CbmStsTrack*>(trkl)->GetNofMvdHits();
1811 break;
1812 case ECbmModuleId::kSts:
1813 if (trkl) nhits = static_cast<CbmStsTrack*>(trkl)->GetNofStsHits();
1814 break;
1816 case ECbmModuleId::kTrd:
1817 if (trkl) nhits = trkl->GetNofHits();
1818 break;
1819 case ECbmModuleId::kTof:
1820 nhits = 1; /* one is maximum */
1821 break;
1823 if (ring) nhits = ring->GetNofHits();
1824 break;
1825 default: continue;
1826 }
1827
1828 // loop over all reconstructed hits
1829 CbmHit* hit = 0x0;
1830 CbmMatch* mtch = 0x0;
1831 FairMCPoint* pnt = 0x0;
1832 for (Int_t ihit = 0; ihit < nhits; ihit++) {
1833 Int_t idx = -1;
1834 switch (idet) {
1835 case ECbmModuleId::kMvd: idx = static_cast<CbmStsTrack*>(trkl)->GetMvdHitIndex(ihit); break;
1836 case ECbmModuleId::kSts: idx = static_cast<CbmStsTrack*>(trkl)->GetStsHitIndex(ihit); break;
1838 case ECbmModuleId::kTrd: idx = trkl->GetHitIndex(ihit); break;
1839 case ECbmModuleId::kTof: hit = track->GetTofHit(); break;
1840 case ECbmModuleId::kRich: idx = ring->GetHit(ihit); break;
1841 default: continue;
1842 }
1843 // get hit
1844 if (idet != ECbmModuleId::kTof && idx > -1) { hit = dynamic_cast<CbmHit*>(hits->At(idx)); }
1845 if (!hit) continue;
1846
1847 // fill rec hit variables
1848 PairAnalysisVarManager::Fill(hit, values);
1849 Bool_t trueHit = kTRUE;
1850 Bool_t fakeHit = kTRUE;
1851 // access to mc points
1852 if ((mtch = hit->GetMatch()) && ev->GetPoints(idet)) {
1853 Int_t nlinks = mtch->GetNofLinks();
1854
1855 // check if mc point corresponds to the matched track (true or fake pnt)
1856 // DEFINITION: always a fake hit if you link to >1 mc points?
1857 // trueHit = (pnt->GetTrackID() == mctrk && mctrk>-1 && nlinks==1);
1858
1859 // fill MC hit variables
1860 // NOTE: the sum of all linked mc points is stored, you have to normlize to the mean
1861 // loop over all linked mc points
1862 for (Int_t iLink = 0; iLink < nlinks; iLink++) {
1863 pnt = static_cast<FairMCPoint*>(ev->GetPoints(idet)->At(mtch->GetLink(iLink).GetIndex()));
1864
1865 // Fill the MC hit variables
1866 if (!iLink) PairAnalysisVarManager::Fill(pnt, values);
1867 else
1869
1870 // hit type defintion
1871 if (!pnt) trueHit = kFALSE;
1872 else if (mc) {
1873 Int_t lbl = pnt->GetTrackID();
1874 Int_t lblM = mc->GetMothersLabel(lbl);
1875 Int_t lblG = mc->GetMothersLabel(lblM);
1876 if (lbl != mctrk && lblM != mctrk && lblG != mctrk) trueHit = kFALSE;
1877 else
1878 fakeHit = kFALSE;
1879 }
1880
1881 } //end links
1882
1883 } //end match found
1884
1885 // fill rec hit histos
1886 if (hitClass) fHistos->FillClass(className, values);
1887 if (hitClass2) fHistoArray->FillClass(className, values);
1888 // fill true, distorted or fake hit histos
1889 if (trueHit) {
1890 if (hitClass) fHistos->FillClass(className + "_true", values);
1891 if (hitClass2) fHistoArray->FillClass(className + "_true", values);
1892 }
1893 else if (fakeHit) {
1894 if (hitClass) fHistos->FillClass(className + "_fake", values);
1895 if (hitClass2) fHistoArray->FillClass(className + "_fake", values);
1896 }
1897 else {
1898 if (hitClass) fHistos->FillClass(className + "_dist", values);
1899 if (hitClass2) fHistoArray->FillClass(className + "_dist", values);
1900 }
1901 // check and fill mc signal histos
1902 for (Int_t isig = 0; isig < nsig; isig++) {
1903 sigName = className + "_" + fSignalsMC->At(isig)->GetName();
1904 if (fillMC->TestBitNumber(isig)) { // track is MC signal truth
1905 if (hitClassMC.TestBitNumber(isig)) fHistos->FillClass(sigName, values);
1906 if (hitClassMChf.TestBitNumber(isig)) fHistoArray->FillClass(sigName, values);
1907 }
1908 }
1909
1910 // reset pointer (needed for tof)
1911 hit = 0x0;
1912 } // rec hit loop
1913 } // det loop
1914}
TClonesArray * tracks
TClonesArray * points
ECbmModuleId
Definition CbmDefs.h:39
@ kMvd
Micro-Vertex Detector.
@ kTrd
Transition Radiation Detector.
@ kTof
Time-of-flight Detector.
@ kSts
Silicon Tracking System.
@ kRef
Reference plane.
@ kMuch
Muon detection system.
@ kNofSystems
For loops over active systems.
@ kRich
Ring-Imaging Cherenkov Detector.
static vector< vector< QAHit > > hits
ClassImp(PairAnalysis) const char *PairAnalysis
TList * GetCuts() const
THashList * GetHistogramList()
virtual void AddCuts(AnalysisCuts *cuts)
virtual UInt_t IsSelected(Double_t *const values)
CbmMatch * GetMatch() const
Definition CbmHit.h:75
double GetCharge() const
Charge of the associated particle.
int32_t GetMotherId() const
Definition CbmMCTrack.h:69
int32_t GetNPoints(ECbmModuleId detId) const
const CbmLink & GetLink(int32_t i) const
Definition CbmMatch.h:39
int32_t GetNofLinks() const
Definition CbmMatch.h:42
const CbmLink & GetMatchedLink() const
Definition CbmMatch.h:41
uint32_t GetHit(int32_t i) const
Definition CbmRichRing.h:39
int32_t GetNofHits() const
Definition CbmRichRing.h:37
virtual int32_t GetNofHits() const
Definition CbmTrack.h:58
int32_t GetHitIndex(int32_t iHit) const
Definition CbmTrack.h:59
void Fill(UInt_t mask, TObject *obj, UInt_t addIdx=0)
void AddTrackFilter2(AnalysisFilter *trkFilter2)
void AddTrackFilter(AnalysisFilter *trkFilter)
void AddPrePairFilter(AnalysisFilter *pairFilter)
void AddPairFilter(AnalysisFilter *pairFilter)
void AddEventFilter(AnalysisFilter *eventFilter)
void FillAll(TObject *obj, UInt_t addIdx=0)
void AddTrackFilterMC(AnalysisFilter *trkFilterMC)
TClonesArray * GetHits(ECbmModuleId det) const
CbmVertex * GetPrimaryVertex() const
Int_t GetNumberOfTracks() const
PairAnalysisTrack * GetTrack(UInt_t pos)
TClonesArray * GetPoints(ECbmModuleId det) const
void Fill(Int_t pairIndex, const PairAnalysisPair *particle)
void FillClass(const char *histClass, const Double_t *values)
Int_t GetLabelMotherWithPdg(const PairAnalysisPair *pair, Int_t pdgMother)
static PairAnalysisMC * Instance()
Int_t GetMothersLabel(Int_t daughterLabel) const
CbmMCTrack * GetMCTrackFromMCEvent(Int_t label) const
Bool_t IsMCTruth(const PairAnalysisPair *pair, const PairAnalysisSignalMC *signalMC) const
void Init(const PairAnalysis *papa=0x0)
Int_t FindBin(const Double_t values[], TString *dim=0x0)
void Fill(const PairAnalysisEvent *ev, PairAnalysis *papa)
void SetPdgCode(Int_t pdgCode)
void SetType(Char_t type)
virtual void SetTracks(PairAnalysisTrack *const particle1, Int_t pid1, PairAnalysisTrack *const particle2, Int_t pid2)=0
virtual void RotateTrack(PairAnalysisTrackRotator *rot)=0
void SetLabel(Int_t label)
void SetKFUsage(Bool_t KFUsage)
EBranchRelation GetMothersRelation() const
TLorentzVector * GetMomentum()
CbmTrackMatchNew * GetTrackMatch(ECbmModuleId det) const
CbmTrack * GetTrack(ECbmModuleId det) const
void SetMassHypo(Int_t pdg1, Int_t pdg2, Bool_t refitMassAssump)
void SetWeight(Double_t wght)
CbmTofHit * GetTofHit() const
void SetPdgCode(Int_t pdg)
Double_t GetWeight() const
Short_t Charge() const
CbmRichRing * GetRichRing() const
static void SetFillMap(TBits *map)
static void SetEvent(PairAnalysisEvent *const ev)
static void FillVarMCParticle(const CbmMCTrack *p1, const CbmMCTrack *p2, Double_t *const values)
static void SetValue(ValueTypes var, Double_t val)
static PairAnalysisEvent * GetCurrentEvent()
static void Fill(const TObject *particle, Double_t *const values)
static void FillSum(const TObject *particle, Double_t *const values)
AnalysisFilter fPairPreFilterLegs
void FillCutStepHistograms(AnalysisFilter *filter, UInt_t cutmask, PairAnalysisTrack *trk, const Double_t *values)
PairAnalysisTrackRotator * fTrackRotator
Pair candidate arrays.
static const char * fgkPairClassNames[8]
AnalysisFilter fFinalTrackFilter
TObjArray * fSignalsMC
TObjArray * fPairCandidates
Selected track candidates.
PairAnalysisHF * fHistoArray
Bool_t fPreFilterUnlikeOnly
void FillHistograms(const PairAnalysisEvent *ev, Bool_t pairInfoOnly=kFALSE)
void FillHistogramsHits(const PairAnalysisEvent *ev, TBits *fillMC, PairAnalysisTrack *track, Bool_t trackIsLeg, Double_t *values)
AnalysisFilter fPairFilter
THashList * fCutStepHistos
AnalysisFilter fPairFilterMC
void FillTrackArrays(PairAnalysisEvent *const ev)
AnalysisFilter fPairPreFilter
void ClearArrays()
virtual ~PairAnalysis()
Int_t GetPairIndex(Int_t arr1, Int_t arr2) const
Bool_t FillMCHistograms(Int_t label1, Int_t label2, Int_t nSignal)
TBits * fUsedVars
Bool_t fStoreRotatedPairs
PairAnalysis()
pair prefilter leg cut logic
void FillHistogramsMC(const PairAnalysisEvent *ev, PairAnalysisEvent *ev1)
Bool_t fPreFilterAllSigns
ECutType fCutType
void AddSignalMC(PairAnalysisSignalMC *signal)
static constexpr Int_t fLegTypes
TObjArray * PairArray(Int_t i)
Bool_t fNoPairing
Bool_t fProcessLS
TObjArray fTracks[4]
Bool_t fEventProcess
AnalysisFilter fTrackFilter
void PairPreFilter(Int_t arr1, Int_t arr2, TObjArray &arrTracks1, TObjArray &arrTracks2)
PairAnalysisHistos * fHistos
void FillHistogramsTracks(TObjArray **tracks)
void FillHistogramsFromPairArray(Bool_t pairInfoOnly=kFALSE)
PairAnalysisCutQa * fQAmonitor
PairAnalysisMixingHandler * fMixing
void FillCutStepHistogramsMC(AnalysisFilter *filter, UInt_t cutmask, Int_t label, const Double_t *values)
Bool_t fRefitMassAssump
AnalysisFilter fTrackFilterMC
static const char * fgkTrackClassNames[2]
Bool_t fDontClearArrays
void FillPairArrays(Int_t arr1, Int_t arr2)
AnalysisFilter fEventFilter
static constexpr Int_t fNTypes
void FilterTrackArrays(TObjArray &arrTracks1, TObjArray &arrTracks2)
void InitPairCandidateArrays()
void Process(TObjArray *arr)
void FillHistogramsPair(PairAnalysisPair *pair, Bool_t fromPreFilter=kFALSE)
TString GetDetName(ECbmModuleId det)