58 if (f.IsOpen() && !f.IsZombie()) {
66 cout <<
"Num. of events to be used# " <<
dstChain->GetEntries() << endl;
79 UInt_t runId =
fHeader->GetRunId();
80 cout <<
"runId = " << runId << endl;
84 if (runId && runId < 10000) {
108 for (
Int_t iEvent = 0; iEvent < 2000; iEvent++) {
112 for (
Int_t iHit = 0; iHit <
fGemHits->GetEntriesFast(); iHit++) {
115 new ((*GEM)[
GEM->GetEntriesFast()])
123 new ((*SILICON)[SILICON->GetEntriesFast()])
141 nStats =
gem->GetNStations();
142 nMods =
gem->GetStation(0)->GetNModules();
144 nStats =
silicon->GetNStations();
145 nMods =
silicon->GetStation(2)->GetNModules();
148 set<Double_t>*** x =
new set<Double_t>**[nStats];
149 set<Double_t>*** y =
new set<Double_t>**[nStats];
150 set<Double_t>*** z =
new set<Double_t>**[nStats];
152 for (
Int_t iStat = 0; iStat < nStats; iStat++) {
153 x[iStat] =
new set<Double_t>*[nMods];
154 y[iStat] =
new set<Double_t>*[nMods];
155 z[iStat] =
new set<Double_t>*[nMods];
157 for (
Int_t iMod = 0; iMod < nMods; iMod++) {
158 x[iStat][iMod] =
new set<Double_t>();
159 y[iStat][iMod] =
new set<Double_t>();
160 z[iStat][iMod] =
new set<Double_t>();
166 for (
Int_t iStat = 0; iStat < nStats; iStat++)
167 for (
Int_t iMod = 0; iMod < nMods; iMod++)
171 x[iStat][iMod]->insert(g->
xyz().
X());
172 y[iStat][iMod]->insert(g->
xyz().
Y());
173 z[iStat][iMod]->insert(g->
xyz().
Z());
179 for (
Int_t iStat = 0; iStat < nStats; iStat++)
180 for (
Int_t iMod = 0; iMod < nMods; iMod++) {
182 Double_t minX = *min_element(x[iStat][iMod]->begin(), x[iStat][iMod]->end());
183 Double_t maxX = *max_element(x[iStat][iMod]->begin(), x[iStat][iMod]->end());
185 Double_t minY = *min_element(y[iStat][iMod]->begin(), y[iStat][iMod]->end());
186 Double_t maxY = *max_element(y[iStat][iMod]->begin(), y[iStat][iMod]->end());
188 vector<Double_t> borders{minX, maxX, minY, maxY};
190 if (z[iStat][iMod]->size()) {
191 fStatZ[make_pair(iStat + shift, iMod)] = *z[iStat][iMod]->begin();
197 for (
Int_t iStat = 0; iStat < nStats; iStat++)
198 for (
Int_t iMod = 0; iMod < nMods; iMod++) {
199 delete x[iStat][iMod];
200 delete y[iStat][iMod];
201 delete z[iStat][iMod];
235 if (f.IsOpen() && !f.IsZombie()) {
243 cout <<
"Num. of events to be used# " <<
dstChain->GetEntries() << endl;
256 UInt_t runId =
fHeader->GetRunId();
260 if (runId && runId < 10000) {
273 fField = fAna->GetField();
280 map<
Int_t, vector<pair<Double_t, Double_t>>> gYr,
281 map<
Int_t, vector<pair<Double_t, Double_t>>> sYr)
287 auto silFirst =
fStatZ.begin();
288 auto gemLast =
fStatZ.rbegin();
291 Double_t zBegin = silFirst->second - 5.;
292 Double_t zEnd = gemLast->second + 5.;
294 Int_t nModsSilicon = 0;
295 for (
Int_t iStat = 0; iStat <
silicon->GetNStations(); iStat++)
296 for (
Int_t iMod = 0; iMod <
silicon->GetStation(iStat)->GetNModules(); iMod++)
299 auto gemFirst = next(
fStatZ.begin(), nModsSilicon);
300 auto silLast = next(
fStatZ.begin(), nModsSilicon - 1);
302 Double_t borderSilGem = .5 * (gemFirst->second + silLast->second);
309 if (iEvent % 100000 == 0)
310 cout <<
"Event# " << iEvent << endl;
327 vector<BmnHit*> hits;
330 if (track->GetSilTrackIndex() != -1) {
348 for (
BmnHit* hit : hits) {
350 if (Z > zBegin && Z < zEnd)
358 Int_t nHitsPerSilicon = 0;
359 Int_t nHitsPerGem = 0;
361 for (
BmnHit* hit : hits) {
362 if (hit->GetZ() < borderSilGem)
380 FairTrackParam par = *first;
381 for (
BmnHit* hit : hits) {
382 BmnStatus status =
fKalman->TGeoTrackPropagate(&par, hit->GetZ(), last->GetQp() > 0. ? 2212 : -211,
383 nullptr,
nullptr,
kTRUE);
388 fKalman->Update(&par, hit, chi2);
393 fKalman->TGeoTrackPropagate(&par, z, last->GetQp() > 0. ? 2212 : -211,
nullptr,
nullptr,
kTRUE);
397 Double_t xPredictedAtVp = par.GetX();
398 Double_t yPredictedAtVp = par.GetY();
402 TMath::Sqrt((xPredictedAtVp - x) * (xPredictedAtVp - x) + (yPredictedAtVp - y) * (yPredictedAtVp - y));
404 if (distXY > distCut)
408 pair<FairTrackParam, FairTrackParam> trackParams = make_pair(parPredicted, *track->
GetParamLast());
411 Bool_t hitsPerStation[nHitsMax];
412 for (
Int_t iHit = 0; iHit < nHitsMax; iHit++)
413 hitsPerStation[iHit] =
kFALSE;
416 for (
BmnHit* hit : hits) {
420 Int_t stat = hit->GetStation();
421 hitsPerStation[stat + shift] =
kTRUE;
429 v.assign(hitsPerStation, hitsPerStation + nHitsMax);
447 input1.
det =
"SILICON";
458 vector<BmnHit*> hits,
459 pair<FairTrackParam, FairTrackParam> params,
460 pair<Double_t, Double_t>& xyPred)
465 vector<BmnHit*> hitsBeforeCurrentZ;
467 for (
BmnHit* hit : hits) {
473 hitsBeforeCurrentZ.push_back(hit);
476 for (
BmnHit* hit : hitsBeforeCurrentZ) {
478 &(params.first), hit->GetZ(), params.second.GetQp() > 0. ? 2212 : -211,
nullptr,
nullptr,
kTRUE);
483 fKalman->Update(&(params.first), hit, chi2);
486 BmnStatus status =
fKalman->TGeoTrackPropagate(&(params.first), Z, params.second.GetQp() > 0. ? 2212 : -211,
487 nullptr,
nullptr,
kTRUE);
492 Double_t xPred = params.first.GetX();
493 Double_t yPred = params.first.GetY();
495 xyPred.first = xPred;
496 xyPred.second = yPred;
499 Int_t stat = it.first.first;
500 Int_t mod = it.first.second;
507 if (xPred > xMin && xPred < xMax && yPred > yMin && yPred < yMax && (stat == iHit)) {
518 const Int_t hitCut = (input.det.Contains(
"GEM") ? input.nGemCut : input.nSilCut);
519 const Int_t shift = (input.det.Contains(
"GEM") ? -
silicon->GetNStations() : 0);
522 Int_t nHitsPerDetector = 0;
524 auto itBegin = input.usedStats.begin();
525 auto itEnd = input.usedStats.end();
526 auto itMiddle = next(input.usedStats.begin(),
silicon->GetNStations());
528 if (input.det.Contains(
"SILICON"))
529 for (
auto it = itBegin; it < itMiddle; it++) {
534 for (
auto it = itMiddle; it != itEnd; it++) {
539 if (!nHitsPerDetector)
542 for (
Int_t iHit = input.sfIndices.first; iHit < input.sfIndices.second; iHit++) {
544 if (nHitsPerDetector == hitCut && input.usedStats.at(iHit))
547 if (input.usedStats.at(iHit) && nHitsPerDetector != hitCut) {
552 for (
BmnHit* hit : input.hits) {
557 if (stat == (iHit + shift)) {
565 for (
Int_t iCont = 0; iCont < input.eContainers->GetEntriesFast(); iCont++) {
568 if (effCont->
Detector() == input.det && effCont->
Station() == (iHit + shift)) {
576 eff->
efficiency(((input.flParams.second.GetQp() > 0.) ? 0 : 1),
"", input.det)
577 ->
Fill(
kTRUE, hit0->GetX(), hit0->GetY());
579 if (input.det.Contains(
"SILICON"))
580 eff->
efficiency(((input.flParams.second.GetQp() > 0.) ? 0 : 1),
"xP", input.det)
584 FindBin(hit0->GetY(), stat, ((input.det.Contains(
"GEM")) ? input.gPairMap : input.sPairMap));
586 eff->
efficiency4range(yBin, ((input.flParams.second.GetQp() > 0.) ? 0 : 1), input.det)
591 pair<Double_t, Double_t> xyPred;
595 ((input.det.Contains(
"GEM")) ? input.gPairMap : input.sPairMap));
600 for (
Int_t iCont = 0; iCont < input.eContainers->GetEntriesFast(); iCont++) {
603 if (effCont->
Detector() == input.det && effCont->
Station() == (iHit + shift)) {
611 eff->
efficiency(((input.flParams.second.GetQp() > 0.) ? 0 : 1),
"", input.det)
614 if (input.det.Contains(
"SILICON"))
615 eff->
efficiency(((input.flParams.second.GetQp() > 0.) ? 0 : 1),
"xP", input.det)
618 eff->
efficiency4range(yBin, ((input.flParams.second.GetQp() > 0.) ? 0 : 1), input.det)
630 if (hit->GetZ() <
gem->GetStation(0)->GetZPosition() - 5.)
639 const Int_t nLayers =
gem->GetStation(stat)->GetModule(mod)->GetNStripLayers();
640 Bool_t isHitInside[nLayers];
642 for (
Int_t iLayer = 0; iLayer < nLayers; iLayer++) {
647 if (isHitInside[0] && isHitInside[1])
649 else if (isHitInside[2] && isHitInside[3])
R__EXTERN TStyle * gStyle
BmnStatus
Define status code.
void fillEfficiency(inputForEfficiency &&in)
BmnKalmanFilter * fKalman
BmnNewFieldMap * fMagField
TClonesArray * fInnerHits
void Efficiency(TClonesArray *, TClonesArray *effSilicon, map< Int_t, vector< pair< Double_t, Double_t > > >, map< Int_t, vector< pair< Double_t, Double_t > > >)
BmnGemStripStationSet * gem
TClonesArray * fSiliconTracks
map< pair< Int_t, Int_t >, Double_t > fStatZ
TString GetDetector(BmnHit *hit)
BmnSiliconStationSet * silicon
TClonesArray * fGemTracks
Double_t FindZ(Int_t stat)
Int_t FindBin(Double_t y, Int_t stat, map< Int_t, vector< pair< Double_t, Double_t > > > pairMap)
TClonesArray * fSiliconHits
void FillAcceptanceMaps(TClonesArray *, TClonesArray *SILs=nullptr)
map< pair< Int_t, Int_t >, vector< Double_t > > fStatAcceptance
TString GetGemZone(BmnHit *)
Bool_t isVirtualHitInAcceptance(Int_t, vector< BmnHit * >, pair< FairTrackParam, FairTrackParam >, pair< Double_t, Double_t > &)
TClonesArray * fGlobTracks
Bool_t IsPointInsideStripLayer(Double_t x, Double_t y)
Int_t GetGemTrackIndex() const
Accessors.
Short_t GetStation() const
BmnGemStripStationSet * GetGemGeometry()
BmnSiliconStationSet * GetSiliconGeometry()
Field Map does not have any symmetries w.r.t coordinate axises.
FairTrackParam * GetParamFirst()
Int_t GetNHits() const
Accessors.
Int_t GetHitIndex(Int_t iHit) const
FairTrackParam * GetParamLast()
Double_t GetX() const
Accessors.
TEfficiency * efficiency(Int_t qpBin=0, TString hType="", TString det="GEM")
TEfficiency * efficiency4range(Int_t range, Int_t qpBin=0, TString det="GEM")
void Fill(Bool_t bPassed, Double_t x, Double_t y=0, Double_t z=0)
Int_t GetEntriesFast() const
TObject * UncheckedAt(Int_t i) const
friend friend class TClonesArray
const char * Data() const
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
UniRun.h Class for the table: run_.
static UniRun * GetRun(int period_number, int run_number)
get run from the database
double * GetFieldVoltage()
get field voltage of the current run
Double_t Sqrt(Double_t x)