13 #include "FairRootManager.h"
25 #include "FairLogger.h"
26 #include "TDatabasePDG.h"
27 #include "TParticlePDG.h"
37 #include "TObjArray.h"
42 #include "FairRuntimeDb.h"
56 : FairTask(name, verbose)
72 , fPointInfos(new TClonesArray(
"CbmMuchPointInfo", 10))
75 , fhChargeEnergyLog(NULL)
76 , fhChargeEnergyLogPi(NULL)
77 , fhChargeEnergyLogPr(NULL)
78 , fhChargeEnergyLogEl(NULL)
79 , fhChargeTrackLength(NULL)
80 , fhChargeTrackLengthPi(NULL)
81 , fhChargeTrackLengthPr(NULL)
82 , fhChargeTrackLengthEl(NULL)
84 , fhChargePr_1GeV_3mm(NULL)
105 , fhPointsInCluster(NULL)
106 , fhDigisInCluster(NULL)
107 , fhHitsInCluster(NULL)
117 , fPointsUnderCounted(0)
118 , fPointsOverCounted(0)
119 , fOccupancyQaOn(kTRUE)
121 , fDigitizerQaOn(kTRUE)
122 , fStatisticsQaOn(kTRUE)
123 , fClusterDeconvQaOn(kTRUE)
124 , fPrintToFileOn(kTRUE)
143 FairRootManager* fManager = FairRootManager::Instance();
144 fMCTracks = (TClonesArray*) fManager->GetObject(
"MCTrack");
145 fPoints = (TClonesArray*) fManager->GetObject(
"MuchPoint");
146 fHits = (TClonesArray*) fManager->GetObject(
"MuchPixelHit");
149 fClusters = (TClonesArray*) fManager->GetObject(
"MuchCluster");
163 TObjArray* stations = (TObjArray*)
f->Get(
"stations");
166 printf(
"Init: fNstations = %i\n",
fNstations);
193 #define BINNING_CHARGE 100, 0, 3.0
194 #define BINNING_LENGTH 100, 0, 2.5
195 #define BINNING_CHARGE_LOG 100, 4, 8
196 #define BINNING_ENERGY_LOG 100, -0.5, 4.5
197 #define BINNING_ENERGY_LOG_EL 100, -0.5, 4.5
207 "Charge vs energy (log.) for all tracks",
212 "Charge vs energy (log.) for pions",
217 "Charge vs energy (log.) for protons",
222 "Charge vs energy (log.) for electrons",
227 "Charge vs length for all tracks",
232 "Charge vs length for pions",
236 "Charge vs length for proton",
241 "Charge vs length for electrons",
246 new TH1D(
"hChargePr_1GeV_3mm",
"Charge for 1 GeV protons",
BINNING_CHARGE);
262 fhChargeLog->GetXaxis()->SetTitle(
"Lg (Charge) [Number of electrons]");
278 for (Int_t
i = 0;
i < 8;
i++) {
281 histo->GetYaxis()->SetDecimals(kTRUE);
282 histo->GetYaxis()->SetTitleOffset(1.4);
283 histo->GetYaxis()->CenterTitle();
284 histo->GetYaxis()->SetTitle(
"Charge [10^{6} electrons]");
286 histo->GetXaxis()->SetTitle(
"Lg E_{kin} [MeV]");
288 histo->GetXaxis()->SetTitle(
"Track length [cm]");
302 Double_t rMax = station->
GetRmax();
303 Double_t rMin = station->
GetRmin();
305 new TH1D(Form(
"hPadsTotal%i",
i + 1),
306 Form(
"Number of pads vs radius: station %i;Radius [cm]",
i + 1),
311 Form(
"hPadsFired%i",
i + 1),
312 Form(
"Number of fired pads vs radius: station %i;Radius [cm]",
i + 1),
317 Form(
"hOccupancy%i",
i + 1),
318 Form(
"Occupancy vs radius: station %i;Radius [cm];Occupancy, %%",
i + 1),
324 new TH1D(Form(
"hCharge%i",
i + 1),
325 Form(
"Charge : station %i; Charge(1e4); Count",
i + 1),
330 new TH1I(Form(
"hPointsInCluster%i",
i + 1),
331 Form(
"Points in Cluster : Station %i ",
i + 1),
336 new TH1I(Form(
"hDigisInCluster%i",
i + 1),
337 Form(
"Digis in Cluster : Station %i ",
i + 1),
342 Form(
"Hits in Cluster : Station %i ",
i + 1),
349 for (vector<CbmMuchModule*>::iterator im = modules.begin();
356 vector<CbmMuchPad*> pads = module->
GetPads();
357 for (UInt_t ip = 0; ip < pads.size(); ip++) {
360 Double_t x0 = pad->
GetX();
361 Double_t y0 = pad->
GetY();
362 Double_t r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
379 Int_t nTotChannels = 0;
380 for (Int_t iStation = 0; iStation <
fNstations; iStation++) {
397 printf(
"%i\t\t| %i\t\t\n", iStation + 1, nChannels);
398 nTotChannels += nChannels;
400 printf(
"---------------------------------------------------------------------"
401 "--------------------\n");
402 printf(
" Total:\t\t| %i\t\t\n", nTotChannels);
403 printf(
"====================================================================="
404 "====================\n");
437 "hPullX",
"Pull distribution X;(x_{RC} - x_{MC}) / dx_{RC}", 500, -5, 5);
439 "hPullY",
"Pull distribution Y;(y_{RC} - y_{MC}) / dy_{RC}", 500, -5, 5);
441 "hPullT",
"Pull distribution T;(t_{RC} - t_{MC}) / dt_{RC}", 120, -10, 50);
445 "hResidualX",
"Residual distribution X;(x_{RC} - x_{MC})(cm)", 500, -5, 5);
447 "hResidualY",
"Residual distribution Y;(y_{RC} - y_{MC})(cm)", 500, -5, 5);
449 "Residual distribution T;(t_{RC} - t_{MC})(ns)",
455 "Number of fired pads vs pad area:area:n pads",
483 LOG(info) <<
"Event: " <<
fEvent;
494 for (
int i = 0;
i <
fPoints->GetEntriesFast();
i++) {
497 if (stId != 0)
continue;
499 if (layerId != 0)
continue;
500 printf(
"point %4i xin=%6.2f yin=%6.2f xout=%6.2f yout=%6.2f zin=%6.2f\n",
508 "%7.3f %7.3f %7.3f %7.3f\n",
522 if (stId != 0)
continue;
524 if (layerId != 0)
continue;
527 if (!module)
continue;
529 Double_t x0 = pad->
GetX();
530 Double_t y0 = pad->
GetY();
531 UInt_t charge = digi->
GetAdc();
532 printf(
"digi %4i x0=%5.1f y0=%5.1f charge=%3i\n",
i, x0, y0, charge);
533 fprintf(
padsFile,
"%5.1f %5.1f %3i\n", x0, y0, charge);
541 printf(
"FinishTask\n");
545 gStyle->SetPaperSize(20, 20);
558 Double_t errors[100];
559 for (Int_t
i = 0;
i < 100;
i++)
576 printf(
"===================================\n");
577 printf(
"PullsQa:\n");
592 fhPullX->GetFunction(
"gaus")->SetLineWidth(3);
593 fhPullX->GetFunction(
"gaus")->SetLineColor(kRed);
601 fhPullY->GetFunction(
"gaus")->SetLineWidth(3);
602 fhPullY->GetFunction(
"gaus")->SetLineColor(kRed);
610 fhPullT->GetFunction(
"gaus")->SetLineWidth(3);
611 fhPullT->GetFunction(
"gaus")->SetLineColor(kRed);
624 fhResidualX->GetFunction(
"gaus")->SetLineColor(kRed);
633 fhResidualY->GetFunction(
"gaus")->SetLineColor(kRed);
642 fhResidualT->GetFunction(
"gaus")->SetLineColor(kRed);
719 printf(
"===================================\n");
720 printf(
"OccupancyQa:\n");
788 printf(
"===================================\n");
789 printf(
"DigitizerQa:\n");
791 TF1* fit_el =
new TF1(
"fit_el",
LandauMPV, -0.5, 4.5, 1);
792 fit_el->SetParameter(0, 0.511);
793 fit_el->SetLineWidth(2);
794 fit_el->SetLineColor(kBlack);
796 TF1* fit_pi =
new TF1(
"fit_pi",
LandauMPV, -0.5, 4.5, 1);
797 fit_pi->SetParameter(0, 139.57);
798 fit_pi->SetLineWidth(2);
799 fit_pi->SetLineColor(kBlack);
801 TF1* fit_pr =
new TF1(
"fit_pr",
LandauMPV, -0.5, 4.5, 1);
802 fit_pr->SetParameter(0, 938.27);
803 fit_pr->SetLineWidth(2);
804 fit_pr->SetLineColor(kBlack);
813 hLength->SetTitle(
"Track length for all tracks");
814 hLengthEl->SetTitle(
"Track length for electrons");
815 hLengthPi->SetTitle(
"Track length for pions");
816 hLengthPr->SetTitle(
"Track length for protons");
825 for (Int_t
i = 0;
i < 8;
i++) {
875 for (Int_t
i = 10;
i < 14;
i++) {
927 for (Int_t iS = 1; iS <= 10; iS++) {
931 for (Int_t iN = 1; iN <= 10; iN++) {
932 nMean[iS] += iN *
fhNpadsVsS->GetBinContent(iS, iN);
935 if (total > 0) nMean[iS] /= total;
943 printf(
"All tracks: ;");
947 printf(
"------------;");
951 printf(
"Primary: ;");
955 printf(
"Secondary: ;");
959 printf(
"-------------");
963 printf(
"Protons: ;");
965 printf(
"%8i;",
fNpr[
i]);
969 printf(
"%8i;",
fNpi[
i]);
971 printf(
"Electrons: ;");
973 printf(
"%8i;",
fNel[
i]);
977 printf(
"%8i;",
fNmu[
i]);
981 printf(
"%8i;",
fNka[
i]);
987 printf(
"===================================\n");
988 printf(
"StatisticsQa:\n");
1048 Double_t gaz_gain_mean = 1.7e+4;
1049 Double_t scale = 1.e+6;
1050 gaz_gain_mean /= scale;
1051 Double_t mass = par[0];
1052 Double_t
x = TMath::Power(10, lg_x[0]);
1053 return gaz_gain_mean *
MPV_n_e(
x, mass);
1067 for (Int_t
i = 0;
i <
fPoints->GetEntriesFast();
i++) {
1072 Int_t trackId = point->GetTrackID();
1087 TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(pdgCode);
1088 if (!particle || pdgCode == 22 ||
1098 if (pdgCode == 2212)
1100 else if (pdgCode == -211 || pdgCode == 211)
1102 else if (pdgCode == -11 || pdgCode == 11)
1104 else if (pdgCode == -13 || pdgCode == 13)
1106 else if (pdgCode == -321 || pdgCode == 321)
1115 Double_t mass = particle->Mass();
1120 Double_t length = (vOut - vIn).Mag();
1121 Double_t kine =
sqrt(p.Mag2() + mass * mass) - mass;
1123 Int_t motherPdg = 0;
1126 if (motherTrack) motherPdg = motherTrack->
GetPdgCode();
1127 new ((*fPointInfos)[
i])
1149 if (!module)
continue;
1153 LOG(debug) << GetName() <<
" Processing MuchDigi " <<
i <<
" at address "
1163 for (Int_t pt = 0; pt < nofLinks; pt++) {
1181 Double_t kine = 1000 * (info->
GetKine());
1184 if (pdg == 0)
continue;
1185 if (charge <= 0)
continue;
1186 Double_t log_kine = TMath::Log10(kine);
1187 Double_t log_charge = TMath::Log10(charge);
1188 charge = charge / 1.e+4;
1191 Double_t area = info->
GetArea() / nPads;
1201 else if (pdg == 211 || pdg == -211)
1203 else if (pdg == 11 || pdg == -11)
1208 else if (pdg == 211 || pdg == -211)
1210 else if (pdg == 11 || pdg == -11)
1213 if (pdg == 2212 && length > 0.3 && length < 0.32 && kine > 1000
1233 if (!module)
continue;
1241 Double_t x0 = pad->
GetX();
1242 Double_t y0 = pad->
GetY();
1243 r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
1253 Int_t nClusters =
fClusters->GetEntriesFast();
1254 TArrayI hitsInCluster;
1255 TArrayI pointsInCluster;
1256 hitsInCluster.Set(nClusters);
1257 pointsInCluster.Set(nClusters);
1258 cout <<
" start Stat QA " << endl;
1259 for (Int_t
i = 0;
i < nClusters;
i++)
1260 hitsInCluster[
i] = 0;
1261 for (Int_t
i = 0;
i < nClusters;
i++)
1262 pointsInCluster[
i] = 0;
1264 for (Int_t
i = 0;
i <
fHits->GetEntriesFast();
i++) {
1270 hitsInCluster[clusterId]++;
1273 for (Int_t
i = 0;
i < nClusters;
i++) {
1276 map<Int_t, Int_t> map_points;
1285 for (Int_t digiId = 0; digiId < nDigis; digiId++) {
1286 Int_t index = cluster->
GetDigi(digiId);
1293 for (Int_t iRefPoint = 0; iRefPoint < nPoints; iRefPoint++) {
1295 map_points[pointId] = 1;
1298 pointsInCluster[
i] = map_points.size();
1303 for (Int_t
i = 0;
i < nClusters;
i++) {
1311 Int_t nPts = pointsInCluster[
i];
1312 Int_t nHits = hitsInCluster[
i];
1328 for (Int_t
i = 0;
i <
fHits->GetEntriesFast();
i++) {
1340 printf(
" Hit %i, station %i, layer %i ",
i, iStation, iLayer);
1343 Bool_t hit_unique = 1;
1344 for (Int_t j =
i + 1; j <
fHits->GetEntriesFast(); j++) {
1347 if (hit1->
GetRefId() == clusterId) {
1352 if (verbose) printf(
"hit_unique=%i", hit_unique);
1354 if (verbose) printf(
"\n");
1360 Bool_t point_unique = 1;
1372 for (Int_t digiId = 0; digiId < cluster->
GetNofDigis(); digiId++) {
1373 Int_t index = cluster->
GetDigi(digiId);
1378 if (!digi)
continue;
1379 if (index < 0)
continue;
1383 if (!match)
continue;
1387 if (verbose) printf(
" n=%i", match->
GetNofLinks());
1389 printf(
" noise hit");
1400 if (!module)
continue;
1413 pointId = currentPointId;
1417 if (currentPointId != pointId) {
1424 if (verbose) printf(
" point_unique=%i", point_unique);
1425 if (!point_unique) {
1426 if (verbose) printf(
"\n");
1434 if (verbose) printf(
" pointId=%i", pointId);
1439 Double_t tMC = point->GetTime();
1441 Double_t xRC = hit->
GetX();
1442 Double_t yRC = hit->
GetY();
1443 Double_t dx = hit->
GetDx();
1444 Double_t dy = hit->
GetDy();
1446 Double_t tRC = hit->
GetTime();
1451 printf(
"Anomalously small dx\n");
1455 printf(
"Anomalously small dy\n");
1458 fhPullX->Fill((xRC - xMC) / dx);
1459 fhPullY->Fill((yRC - yMC) / dy);
1460 fhPullT->Fill((tRC - tMC) / dt);
1467 if (verbose) printf(
"\n");
1492 Int_t nPoints =
fPoints->GetEntriesFast();
1494 Int_t nClusters =
fClusters->GetEntriesFast();
1495 vector<Int_t> pIndices;
1496 vector<Int_t>::iterator it;
1498 for (Int_t iPoint = 0; iPoint < nPoints; ++iPoint) {
1502 for (Int_t iCluster = 0; iCluster < nClusters; ++iCluster) {
1504 if (!cluster)
continue;
1506 for (Int_t
id = 0;
id < nDigis; ++id) {
1507 Int_t iDigi = cluster->
GetDigi(
id);
1514 if (!match)
continue;
1515 for (Int_t ip = 0; ip < match->
GetNofLinks(); ++ip) {
1517 it = find(pIndices.begin(), pIndices.end(), iPoint);
1518 if (it != pIndices.end())
continue;
1519 pIndices.push_back(iPoint);
1528 Int_t nPoints =
fPoints->GetEntriesFast();
1529 Int_t nTracks =
fMCTracks->GetEntriesFast();
1530 CbmMuchPoint* point = (iPoint < 0 || iPoint > nPoints - 1)
1533 if (!point)
return kFALSE;
1534 Int_t iTrack = point->GetTrackID();
1535 CbmMCTrack* track = (iTrack < 0 || iTrack > nTracks - 1)
1538 if (!track)
return kFALSE;
1539 if (iTrack != 0 && iTrack != 1)
1544 if (TMath::Abs(pdgCode) == 13) {
1553 Int_t nChannels = 0;
1555 for (UInt_t im = 0; im < modules.size(); im++) {
1559 if (!module)
continue;
1569 for (UInt_t im = 0; im < modules.size(); im++) {
1573 if (!module)
continue;
1584 for (UInt_t im = 0; im < modules.size(); im++) {
1588 vector<CbmMuchPad*> pads = module->
GetPads();
1589 for (UInt_t ip = 0; ip < pads.size(); ip++) {
1591 if (pad->
GetDx() < padMinLx) padMinLx = pad->
GetDx();
1592 if (pad->
GetDy() < padMinLy) padMinLy = pad->
GetDy();
1595 return TVector2(padMinLx, padMinLy);
1604 for (UInt_t im = 0; im < modules.size(); im++) {
1608 vector<CbmMuchPad*> pads = module->
GetPads();
1609 for (UInt_t ip = 0; ip < pads.size(); ip++) {
1611 if (pad->
GetDx() > padMaxLx) padMaxLx = pad->
GetDx();
1612 if (pad->
GetDy() > padMaxLy) padMaxLy = pad->
GetDy();
1615 return TVector2(padMaxLx, padMaxLy);
1622 TF1 fPol6(
"fPol6",
"pol6", -5, 10);
1624 logT =
log(Tkin * 0.511 / mass);
1625 if (logT > 9.21034) logT = 9.21034;
1627 return fPol6.EvalPar(&logT,
mpv_e);
1628 }
else if (mass >= 0.1 && mass < 0.2) {
1629 logT =
log(Tkin * 105.658 / mass);
1630 if (logT > 9.21034) logT = 9.21034;
1632 return fPol6.EvalPar(&logT,
mpv_mu);
1634 logT =
log(Tkin * 938.272 / mass);
1635 if (logT > 9.21034) logT = 9.21034;
1637 return fPol6.EvalPar(&logT,
mpv_p);