MpdEmcClusterizerKI.cxx 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481
  1. //-----------------------------------------------------------
  2. //
  3. // Author List:
  4. // D.Peresunko, KI, 2019
  5. //
  6. //-----------------------------------------------------------
  7. // This Class' Header ------------------
  8. #include "MpdEmcClusterizerKI.h"
  9. #include "MpdEmcClusterKI.h"
  10. #include "MpdEmcDigitKI.h"
  11. #include "MpdEmcGeoUtils.h"
  12. #include "MpdEmcSimParams.h"
  13. #include "MpdEmcCalibParams.h"
  14. #include "TFile.h"
  15. #include "TClonesArray.h"
  16. #include "TMath.h"
  17. #include <iostream>
  18. using namespace std;
  19. //__________________________________________________________________________
  20. MpdEmcClusterizerKI::MpdEmcClusterizerKI()
  21. : FairTask("EMC Clusterizer"),
  22. fNumberOfClusters(0),
  23. fDigitsArray(nullptr),
  24. fClustersArray(nullptr),
  25. fSimParams(nullptr),
  26. fGeom(nullptr),
  27. fCalibData(nullptr)
  28. {
  29. }
  30. //__________________________________________________________________________
  31. InitStatus MpdEmcClusterizerKI::Init()
  32. {
  33. // Get ROOT Manager
  34. FairRootManager* ioman = FairRootManager::Instance();
  35. if (ioman == 0) {
  36. LOG(ERROR) << "RootManager not instantiated!";
  37. return kERROR;
  38. }
  39. //Temporary solution: get calibration from file TODO!!!!
  40. if(!fCalibData){
  41. TString CalFile = "$VMCWORKDIR/input/MpdEmcCalib.root";
  42. TFile * f = new TFile(CalFile) ;
  43. fCalibData =(MpdEmcCalibParams*)f->Get("CalibrationDef") ;
  44. cout<<"ECAL: Read out EMC calib data"<<endl;
  45. }
  46. // Get input collection
  47. fDigitsArray = (TClonesArray*)ioman->GetObject("EmcDigit");
  48. if (fDigitsArray == 0) {
  49. LOG(ERROR) << "Array of digits not found!";
  50. return kERROR;
  51. }
  52. // Create and register output array
  53. fClustersArray = new TObjArray(); // Clusters may contain different number of digits=>TObjArray
  54. ioman->Register("EmcCluster", "Emc", fClustersArray, true);
  55. // Class with list of parameters
  56. if (!fSimParams) {
  57. fSimParams = MpdEmcSimParams::GetInstance();
  58. }
  59. if (!fGeom) {
  60. fGeom = MpdEmcGeoUtils::GetInstance();
  61. }
  62. return kSUCCESS;
  63. }
  64. //__________________________________________________________________________
  65. void MpdEmcClusterizerKI::Exec(Option_t* opt)
  66. {
  67. fClustersArray->Delete(); // No possibility to Clear for TObjArray
  68. fNumberOfClusters = 0;
  69. // Prepare Digits: Calibration, bad map, cleanup
  70. PrepareDigits();
  71. // Collect list of clusters
  72. MakeClusters();
  73. // Split clusters with several local maxima if necessary
  74. MakeUnfoldings();
  75. // Evaluate cluster position, dispersion etc.
  76. EvalClusters();
  77. }
  78. //__________________________________________________________________________
  79. void MpdEmcClusterizerKI::PrepareDigits()
  80. {
  81. // Apply (de)calibration
  82. // TODO: Implement Energy and time calibration for data and de-calibrations for MC
  83. // Apply BadMap
  84. // TODO: Implement bad map
  85. // Remove digits below clustering threshold
  86. int n = fDigitsArray->GetEntriesFast();
  87. for (int i = 0; i < n; i++) {
  88. MpdEmcDigitKI* digit = static_cast<MpdEmcDigitKI*>(fDigitsArray->UncheckedAt(i));
  89. if (!digit) { // already removed e.g. by bad map selection
  90. continue;
  91. }
  92. //Calibrate energy
  93. digit->SetE(digit->GetE()*fCalibData->GetGain(digit->GetDetId()));//V
  94. //Calibrate Time TODO!!!
  95. if (digit->GetE() < fSimParams->DigitMinEnergy()) {
  96. fDigitsArray->RemoveAt(i);
  97. }
  98. }
  99. fDigitsArray->Compress();
  100. // Alredy sorted by construction, but ROOT requires this to allow search
  101. fDigitsArray->Sort();
  102. LOG(INFO) << "Digits used for clusters:" << fDigitsArray->GetEntriesFast() << endl;
  103. }
  104. //__________________________________________________________________________
  105. void MpdEmcClusterizerKI::MakeClusters()
  106. {
  107. // Combine digits into cluster according to definition of neighbours in
  108. // MpdEmcGeoUtils::AreNighbours()
  109. // Note, that only digits in same sector contribute to the cluster
  110. // Mark all digits as unused yet
  111. int nDigits = fDigitsArray->GetEntriesFast();
  112. if (nDigits < 1) { // nothing to do
  113. return;
  114. }
  115. bool* digitsUsed = new bool[nDigits]{ false };
  116. int iFirst = 0; // first index of digit which potentially can be a part of cluster
  117. // e.g. first digit in this sector
  118. for (int i = 0; i < nDigits; i++) {
  119. if (digitsUsed[i])
  120. continue;
  121. const MpdEmcDigitKI* digitSeed = static_cast<MpdEmcDigitKI*>(fDigitsArray->UncheckedAt(i));
  122. // is this digit so energetic that start cluster?
  123. MpdEmcClusterKI* clu = nullptr;
  124. int iDigitInCluster = 0;
  125. if (digitSeed->GetE() > fSimParams->ClusteringThreshold()) {
  126. // start a new EMC RecPoint
  127. clu = new MpdEmcClusterKI(digitSeed);
  128. fClustersArray->AddAtAndExpand(clu, fNumberOfClusters);
  129. fNumberOfClusters++;
  130. digitsUsed[i] = true;
  131. iDigitInCluster = 1;
  132. } else {
  133. continue;
  134. }
  135. // Now scan remaining digits in list to find neigbours of our seed
  136. int index = 0;
  137. while (index < iDigitInCluster) { // scan over digits already in cluster
  138. int digitInCluTowerId = clu->GetDigitTowerId(index);
  139. index++;
  140. for (Int_t j = iFirst; j < nDigits;
  141. j++) { // upper limit not really matters, AreNeighbour stops this loop for too far digits
  142. if (digitsUsed[j]) {
  143. continue; // look through remaining digits
  144. }
  145. const MpdEmcDigitKI* digitN = static_cast<MpdEmcDigitKI*>(fDigitsArray->UncheckedAt(j));
  146. // call (digit,digitN) in THAT oder !!!!!
  147. Int_t ineb = fGeom->AreNeighbours(digitInCluTowerId, digitN->GetDetId());
  148. switch (ineb) {
  149. case -1: // too early (e.g. previous sector), do not look before j at subsequent passes
  150. iFirst = j + 1;
  151. continue;
  152. case 0: // not a neighbour
  153. continue;
  154. case 1: // are neighbours
  155. // check if digits have consistent times
  156. // Be VERY careful with this cut!!! Too strict may strongly modify your spectrum
  157. if (std::fabs(digitN->GetTime() - clu->GetTime()) < fSimParams->ClusteringTimeGate()) {
  158. clu->AddDigit(digitN);
  159. iDigitInCluster++;
  160. digitsUsed[j] = true;
  161. }
  162. continue;
  163. case 2: // too far from each other, stop loop
  164. default:
  165. goto nextDigit;
  166. } // switch
  167. }
  168. nextDigit:;
  169. } // loop over cluster
  170. } // energy theshold
  171. delete[] digitsUsed;
  172. }
  173. //__________________________________________________________________________
  174. void MpdEmcClusterizerKI::MakeUnfoldings()
  175. {
  176. //Split cluster if several local maxima are found
  177. if (!fSimParams->UnfoldClusters()) {
  178. return;
  179. }
  180. int* maxAt = new int[fSimParams->NLMMax()]; // NLMMax:Maximal number of local maxima
  181. float* maxAtEnergy = new float[fSimParams->NLMMax()];
  182. float localMaxCut = fSimParams->LocalMaximumCut();
  183. int numberOfNotUnfolded = fNumberOfClusters;
  184. for (int index = 0; index < numberOfNotUnfolded; index++) {
  185. MpdEmcClusterKI* clu = static_cast<MpdEmcClusterKI*>(fClustersArray->At(index));
  186. Int_t nMultipl = clu->GetMultiplicity();
  187. int nMax = clu->GetNumberOfLocalMax(maxAt, maxAtEnergy);
  188. if (nMax > 1) {
  189. UnfoldOneCluster(clu, nMax, maxAt, maxAtEnergy);
  190. fClustersArray->Remove(clu);
  191. fClustersArray->Compress();
  192. index--;
  193. fNumberOfClusters--;
  194. numberOfNotUnfolded--;
  195. } else {
  196. clu->SetNLM(1); // Only one local maximum
  197. }
  198. }
  199. delete[] maxAt;
  200. delete[] maxAtEnergy;
  201. }
  202. //____________________________________________________________________________
  203. void MpdEmcClusterizerKI::UnfoldOneCluster(MpdEmcClusterKI* iniClu, Int_t nMax, int* digitId, float* maxAtEnergy)
  204. {
  205. // Performs the unfolding of a cluster with nMax overlapping showers
  206. // Parameters: iniClu cluster to be unfolded
  207. // nMax number of local maxima found (this is the number of new clusters)
  208. // digitId: index of digits, corresponding to local maxima
  209. // maxAtEnergy: energies of digits, corresponding to local maxima
  210. // Take initial cluster and calculate local coordinates of digits
  211. // To avoid multiple re-calculation of same parameters
  212. int mult = iniClu->GetMultiplicity();
  213. std::vector<double> x(mult);
  214. std::vector<double> y(mult);
  215. std::vector<double> z(mult);
  216. std::vector<double> e(mult);
  217. std::vector<std::vector<double>> eInClusters(mult, std::vector<double>(nMax));
  218. for (int idig = 0; idig < mult; idig++) {
  219. Int_t detID;
  220. Float_t eDigit;
  221. ;
  222. iniClu->GetTransientDigitParams(idig, detID, eDigit);
  223. e[idig] = eDigit;
  224. double lx, ly, lz;
  225. fGeom->DetIdToGlobalPosition(detID, lx, ly, lz);
  226. x[idig] = lx;
  227. y[idig] = ly;
  228. z[idig] = lz;
  229. }
  230. // Coordinates of centers of clusters
  231. std::vector<double> xMax(nMax);
  232. std::vector<double> yMax(nMax);
  233. std::vector<double> zMax(nMax);
  234. std::vector<double> eMax(nMax);
  235. for (int iclu = 0; iclu < nMax; iclu++) {
  236. xMax[iclu] = x[digitId[iclu]];
  237. yMax[iclu] = y[digitId[iclu]];
  238. zMax[iclu] = z[digitId[iclu]];
  239. eMax[iclu] = e[digitId[iclu]];
  240. }
  241. std::vector<double> prop(nMax); // proportion of clusters in the current digit
  242. // Try to decompose cluster to contributions
  243. int nIterations = 0;
  244. bool insuficientAccuracy = true;
  245. while (insuficientAccuracy && nIterations < fSimParams->NMaxIterations()) {
  246. // Loop over all digits of parent cluster and split their energies between daughter clusters
  247. // according to shower shape
  248. for (int idig = 0; idig < mult; idig++) {
  249. double eEstimated = 0;
  250. for (int iclu = 0; iclu < nMax; iclu++) {
  251. prop[iclu] = eMax[iclu] * ShowerShape(sqrt((x[idig] - xMax[iclu]) * (x[idig] - xMax[iclu]) +
  252. (y[idig] - yMax[iclu]) * (y[idig] - yMax[iclu])),
  253. z[idig] - zMax[iclu]);
  254. eEstimated += prop[iclu];
  255. }
  256. if (eEstimated == 0.) { // numerical accuracy
  257. continue;
  258. }
  259. // Split energy of digit according to contributions
  260. for (int iclu = 0; iclu < nMax; iclu++) {
  261. eInClusters[idig][iclu] = e[idig] * prop[iclu] / eEstimated;
  262. }
  263. }
  264. // Recalculate parameters of clusters and check relative variation of energy and absolute of position
  265. insuficientAccuracy = false; // will be true if at least one parameter changed too much
  266. for (int iclu = 0; iclu < nMax; iclu++) {
  267. double oldX = xMax[iclu];
  268. double oldY = yMax[iclu];
  269. double oldZ = zMax[iclu];
  270. double oldE = eMax[iclu];
  271. // new energy, need for weight
  272. eMax[iclu] = 0;
  273. for (int idig = 0; idig < mult; idig++) {
  274. eMax[iclu] += eInClusters[idig][iclu];
  275. }
  276. xMax[iclu] = 0;
  277. yMax[iclu] = 0;
  278. zMax[iclu] = 0.;
  279. double wtot = 0.;
  280. for (int idig = 0; idig < mult; idig++) {
  281. double w = std::max(std::log(eInClusters[idig][iclu] / eMax[iclu]) + fSimParams->LogWeight(), 0.);
  282. xMax[iclu] += x[idig] * w;
  283. yMax[iclu] += y[idig] * w;
  284. zMax[iclu] += z[idig] * w;
  285. wtot += w;
  286. }
  287. if (wtot > 0.) {
  288. xMax[iclu] /= wtot;
  289. yMax[iclu] /= wtot;
  290. zMax[iclu] /= wtot;
  291. }
  292. // Compare variation of parameters
  293. insuficientAccuracy += (std::abs(eMax[iclu] - oldE) > fSimParams->UnfogingEAccuracy());
  294. insuficientAccuracy += (std::abs(xMax[iclu] - oldX) > fSimParams->UnfogingXZAccuracy());
  295. insuficientAccuracy += (std::abs(yMax[iclu] - oldY) > fSimParams->UnfogingXZAccuracy());
  296. insuficientAccuracy += (std::abs(zMax[iclu] - oldZ) > fSimParams->UnfogingXZAccuracy());
  297. }
  298. nIterations++ ;
  299. }
  300. // Iterations finished, add new clusters
  301. for (int iclu = 0; iclu < nMax; iclu++) {
  302. MpdEmcClusterKI* clu = new MpdEmcClusterKI();
  303. fClustersArray->AddAtAndExpand(clu, fNumberOfClusters + iclu);
  304. clu->SetNLM(nMax);
  305. }
  306. for (int idig = 0; idig < mult; idig++) {
  307. Int_t detID;
  308. Float_t eDigit;
  309. ;
  310. iniClu->GetTransientDigitParams(idig, detID, eDigit);
  311. MpdEmcDigitKI testdigit(detID, 0, 0, 0); // test digit
  312. int jdigit = fDigitsArray->BinarySearch(&testdigit); // Look for digit with same detID
  313. if (jdigit == -1) {
  314. LOG(ERROR) << "MpdEmcClusterizerKI::UnfoldOneCluster: Can not find Digit with detID=" << detID;
  315. continue;
  316. }
  317. MpdEmcDigitKI* digit = static_cast<MpdEmcDigitKI*>(fDigitsArray->At(jdigit));
  318. for (int iclu = 0; iclu < nMax; iclu++) {
  319. MpdEmcClusterKI* clu = static_cast<MpdEmcClusterKI*>(fClustersArray->UncheckedAt(fNumberOfClusters + iclu));
  320. clu->AddDigit(digit, eInClusters[idig][iclu]); // Fills geometry and MC infor from Digit,+ correct energy
  321. }
  322. }
  323. fNumberOfClusters += nMax;
  324. }
  325. //__________________________________________________________________________
  326. void MpdEmcClusterizerKI::EvalClusters()
  327. {
  328. // Calculate cluster properties
  329. int n = fClustersArray->GetEntriesFast();
  330. LOG(DEBUG) << "EvalCluProperties: nclu=" << n;
  331. for (int i = 0; i < n; i++) {
  332. MpdEmcClusterKI* clu = static_cast<MpdEmcClusterKI*>(fClustersArray->UncheckedAt(i));
  333. // Remove remaining digits below threshold, e.g. after unfolding
  334. clu->Purify();
  335. // Eval all variables: Energy, CoreEnergy, position, Dispersion,...
  336. clu->EvalAll();
  337. }
  338. }
  339. //__________________________________________________________________________
  340. double MpdEmcClusterizerKI::ShowerShape(double dx, double dz)
  341. {
  342. // Shower shape of EM cluster: proportion of energy in cell as a distancs dx,dz (phi,z directions)
  343. // from the center of gravity of cluster.
  344. // Due to projective geometry may also depend on Z coordinate. TODO: explore Z-dependence
  345. // Parameterization from V.Riabov. TODO: verify with beam-test
  346. double frac = 0;
  347. double x = std::sqrt(dx * dx + dz * dz);
  348. if (x < 0.25) return 0.73;
  349. if (x < 4.5)
  350. {
  351. frac = 7.55666e+000*exp(-2.35773e+000 + 1.23342e-001*x -2.53958e-001*x*x +1.41214e-002*x*x*x );
  352. }
  353. if ( x >= 4.5 && x < 8.33 )
  354. {
  355. frac = 4.26832e-002*exp(1.47109e+000 -2.06712e-001*x -6.60220e-002*x*x +4.88207e-003*x*x*x );
  356. }
  357. if ( x >= 8.33 && x < 2.65*4 )
  358. {
  359. frac = 4.46227e-002*exp( 1.56817e+000 -7.41051e-001*x/4 -4.80563e-001*x/4*x/4 );
  360. }
  361. if ( x >= 2.65*4 )
  362. {
  363. frac = 1.24899e-002*exp( 3.60075e-001 -8.15748e-001*x/4 -3.74305e-002*x/4*x/4*x/4 );
  364. }
  365. if (frac < 1e-24) frac = 1e-24;
  366. return frac;
  367. }
  368. /*
  369. //__________________________________________________________________________
  370. void MpdEmcClusterizerKI::GetPhiTheta(Double_t &phi, Double_t &theta)
  371. {
  372. // Convert COG in units of bins to angles
  373. static Int_t first = 1, offset = 0;
  374. static TSpline3 *phiS, *theS;
  375. if (first) {
  376. // Get phi and theta angles of the tower centers at their inner face
  377. first = 0;
  378. const vector<Double_t> &phis = fEmcGeo->GetPhiRow();
  379. const vector<Double_t> &thes = fEmcGeo->GetThetaBox();
  380. const vector<Double_t> &rhos = fEmcGeo->GetRhoCenterBox();
  381. const vector<Double_t> &zs = fEmcGeo->GetZCenterBox();
  382. Int_t nphi = phis.size();
  383. // Offset due to the fact that the 1'st sector starts at phi = -Phi_sec/2;
  384. offset = nphi / (fEmcGeo->GetNsec()/2 - 1) / 2;
  385. Double_t *phia = new Double_t [nphi];
  386. Double_t *ind = new Double_t [nphi];
  387. for (Int_t j = 0; j < nphi; ++j) {
  388. phia[j] = phis[j];
  389. ind[j] = j;
  390. }
  391. phiS = new TSpline3("grs",ind,phia,nphi); // phi vs ind
  392. delete [] phia;
  393. delete [] ind;
  394. Int_t nthe = thes.size();
  395. Double_t *the = new Double_t [nthe];
  396. Double_t *ind1 = new Double_t [nthe];
  397. Double_t height = fEmcGeo->GetLengthBox(); // tower half-height
  398. for (Int_t j = nthe-1; j >= 0; --j) {
  399. Double_t rho = rhos[j];
  400. Double_t z = zs[j];
  401. Double_t theta1 = thes[j];
  402. if (j < nthe-1 && thes[j] <= thes[j+1]+0.1) theta1 = 180 - theta1;
  403. Double_t costhe = TMath::Cos(theta1*TMath::DegToRad());
  404. Double_t sinthe = TMath::Sin(theta1*TMath::DegToRad());
  405. rho -= height * sinthe;
  406. z -= height * costhe;
  407. the[j] = TMath::ATan2(rho,z) * TMath::RadToDeg();
  408. ind1[j] = j; // - nthe/2;
  409. }
  410. theS = new TSpline3("grs1",ind1,the,nthe); // theta vs ind
  411. delete [] the;
  412. delete [] ind1;
  413. }
  414. phi = phiS->Eval(phi-offset);
  415. if (phi > 180) phi -= 360;
  416. theta = theS->Eval(theta);
  417. }
  418. */
  419. //__________________________________________________________________________
  420. ClassImp(MpdEmcClusterizerKI)