MpdEmcClusterKI.cxx 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494
  1. ////////////////////////////////////////////////////////////////
  2. // //
  3. // MpdEmcHitClusterKI //
  4. // Cluster production for EMC //
  5. // Author List : D.Peresunko., RRCKI, 2019 //
  6. // //
  7. ////////////////////////////////////////////////////////////////
  8. #include "MpdEmcClusterKI.h"
  9. #include <iostream>
  10. #include <numeric>
  11. #include "FairLogger.h"
  12. #include "MpdEmcClusterizerKI.h"
  13. #include "MpdEmcDigitKI.h"
  14. #include "MpdEmcGeoUtils.h"
  15. #include "MpdEmcSimParams.h"
  16. #include "TLorentzVector.h"
  17. #include "TMath.h"
  18. #include "TVector3.h"
  19. using namespace std;
  20. using namespace TMath;
  21. // Function used for sorting primaries
  22. bool compE(const std::pair<int, float>& a, const std::pair<int, float>& b) { return a.second > b.second; }
  23. // ----- Default constructor -------------------------------------------
  24. MpdEmcClusterKI::MpdEmcClusterKI()
  25. : fNDigits(0),
  26. fDitisId(nullptr),
  27. fDigitsE(nullptr),
  28. fNPrimaries(0),
  29. fPrimId(nullptr),
  30. fPrimE(nullptr),
  31. fE(0),
  32. fEcore(0),
  33. fEcore1p(0),
  34. fEcore2p(0),
  35. fTime(0),
  36. fX(0),
  37. fY(0),
  38. fZ(0),
  39. fdPhi(9999),
  40. fdZ(9999),
  41. fTrackId(-1),
  42. fDisp(0),
  43. fChi2(0),
  44. fLambda1(0),
  45. fLambda2(0),
  46. fNExLM(0)
  47. {
  48. }
  49. // ----- Standard constructor ------------------------------------------
  50. MpdEmcClusterKI::MpdEmcClusterKI(const MpdEmcDigitKI* digit)
  51. : fNDigits(0),
  52. fDitisId(nullptr),
  53. fDigitsE(nullptr),
  54. fNPrimaries(0),
  55. fPrimId(nullptr),
  56. fPrimE(nullptr),
  57. fE(0),
  58. fEcore(0),
  59. fEcore1p(0),
  60. fEcore2p(0),
  61. fTime(0),
  62. fX(0),
  63. fY(0),
  64. fZ(0),
  65. fdPhi(9999),
  66. fdZ(9999),
  67. fTrackId(-1),
  68. fDisp(0),
  69. fChi2(0),
  70. fLambda1(0),
  71. fLambda2(0),
  72. fNExLM(0)
  73. {
  74. AddDigit(digit);
  75. }
  76. // ----- Destructor ----------------------------------------------------
  77. MpdEmcClusterKI::~MpdEmcClusterKI()
  78. {
  79. if (fDitisId) {
  80. delete[] fDitisId;
  81. }
  82. if (fDigitsE) {
  83. delete[] fDigitsE;
  84. }
  85. if (fPrimId) {
  86. delete[] fPrimId;
  87. }
  88. if (fPrimE) {
  89. delete[] fPrimE;
  90. }
  91. }
  92. void MpdEmcClusterKI::GetMomentum(TLorentzVector& p, const TVector3* vertex) const
  93. {
  94. // Returm momentum of photon assuming it came from the provided vertex
  95. if (fX == 0. && fY == 0.) { // position not defined
  96. p.SetPxPyPzE(0., 0., 0., 0.);
  97. return;
  98. }
  99. double dx = fX, dy = fY, dz = fZ;
  100. if (vertex) { // calculate direction from vertex
  101. dx -= vertex->X();
  102. dy -= vertex->Y();
  103. dz -= vertex->Z();
  104. }
  105. Double_t r = TMath::Sqrt(dx * dx + dy * dy + dz * dz);
  106. if (r > 0) {
  107. p.SetPxPyPzE(fE * dx / r, fE * dy / r, fE * dz / r, fE);
  108. } else {
  109. p.SetPxPyPzE(0., 0., 0., 0.);
  110. }
  111. }
  112. // ----- Print -----------------------------------------------------------
  113. void MpdEmcClusterKI::Print(const Option_t* opt) const
  114. {
  115. cout << "MpdEmcClusterKI: " << endl;
  116. cout << "\tDeposited energy: " << fE << "\tMean time: " << fTime << "\tDigits: " << fNDigits
  117. << " Rho cluster: " << GetRho() << " Phi cluster: " << GetPhi() << " Z cluster: " << fZ
  118. << " dPhi cluster: " << fdPhi << " dZ cluster: " << fdZ << " index cluster: " << fTrackId << endl;
  119. cout << "\tNumber of tracks in module: " << fNPrimaries << endl;
  120. for (int i = 0; i < fNPrimaries; i++)
  121. cout << " " << fPrimId[i] << ", ";
  122. cout << endl;
  123. }
  124. //____________________________________________________________________________
  125. void MpdEmcClusterKI::AddDigit(const MpdEmcDigitKI* digit, double energy)
  126. {
  127. // Add another digit to cluster:
  128. // increase full energy, re-calculate time if necessary, add MC info, add digit to list of digits
  129. if (energy == 0) { // use full energy of digit, otherwise method is called from Unfolding.
  130. energy = digit->GetE();
  131. }
  132. fE += energy;
  133. // check if this is contribution with highest energy
  134. std::vector<pair<int, float>>::iterator elist = fDigitIDEnergy.begin();
  135. bool found = false;
  136. while (elist != fDigitIDEnergy.end()) {
  137. if ((*elist).second > energy) {
  138. found = true;
  139. break;
  140. }
  141. elist++;
  142. }
  143. if (!found) {
  144. fTime = digit->GetTime();
  145. }
  146. fDigitIDEnergy.emplace_back(digit->GetDetId(), energy);
  147. // now add new primaries
  148. int nPrim = digit->GetNPrimaries();
  149. for (int iprim = 0; iprim < nPrim; iprim++) {
  150. Int_t primId;
  151. Float_t primEdep;
  152. digit->GetPrimary(iprim, primId, primEdep);
  153. //Scale primary energy accosging to proportion
  154. if(digit->GetE()>0){
  155. primEdep*=energy/digit->GetE() ;
  156. std::map<Int_t, Float_t>::iterator it =
  157. fMCTracks.lower_bound(primId); // pointer, where new pair should appear, not necessarilly existing
  158. if ((it == fMCTracks.end()) || ((*it).first != primId)) { // do not exist yet
  159. fMCTracks.insert(it, { primId, primEdep });
  160. } else { // Add energy deposition
  161. (*it).second += primEdep;
  162. }
  163. }
  164. }
  165. }
  166. //____________________________________________________________________________
  167. void MpdEmcClusterKI::EvalAll()
  168. {
  169. // Calculate
  170. // Cluster coordinates
  171. // the shower dispersion, moments and Chi2
  172. // Fills final arrays with accociated digits and energy deposition
  173. if (fE <= 0) { // Energy can not be zero in good cluster, do nothing
  174. return;
  175. }
  176. //Eval position and dispersion
  177. double wtot = 0.;
  178. double x = 0.;
  179. double y = 0.;
  180. double z = 0.;
  181. double phi = 0.;
  182. double dphiphi = 0.;
  183. double dphiz = 0.;
  184. double dzz = 0.;
  185. fDisp = 0.;
  186. fLambda1 = 0;
  187. fLambda2 = 0;
  188. MpdEmcGeoUtils* geom = MpdEmcGeoUtils::GetInstance();
  189. MpdEmcSimParams* simParams = MpdEmcSimParams::GetInstance();
  190. double logWeight = simParams->LogWeight();
  191. // 1) Find covariance matrix elements:
  192. // || dxx dxz ||
  193. // || dxz dzz ||
  194. std::vector<std::pair<int, float>>::iterator digIterator;
  195. digIterator = fDigitIDEnergy.begin();
  196. while (digIterator != fDigitIDEnergy.end()) {
  197. int detID = (*digIterator).first;
  198. float energy = (*digIterator).second;
  199. digIterator++;
  200. double xi, yi, zi, w, phii;
  201. geom->DetIdToGlobalPosition(detID, xi, yi, zi);
  202. if (energy > 0) {
  203. w = TMath::Max(0., logWeight + TMath::Log(energy / fE));
  204. x += w * xi;
  205. y += w * yi;
  206. z += w * zi;
  207. phii = ATan2(yi, xi);
  208. // in case of cluster around 2pi: avoid averaging of delta and 2pi+delta
  209. if (wtot > 0) { // not first digit
  210. double phiCurrent = phi / wtot;
  211. if (phii > phiCurrent + TMath::Pi()) {
  212. phii -= TMath::TwoPi();
  213. }
  214. if (phii < phiCurrent - TMath::Pi()) {
  215. phii += TMath::TwoPi();
  216. }
  217. }
  218. phi += w * phii;
  219. dphiphi += w * phii * phii;
  220. dphiz += w * phii * zi;
  221. dzz += w * zi * zi;
  222. wtot += w;
  223. }
  224. }
  225. if (wtot > 0) {
  226. x /= wtot;
  227. y /= wtot;
  228. z /= wtot;
  229. double r = TMath::Sqrt(x * x + y * y);
  230. phi = phi * r / wtot;
  231. dphiphi = dphiphi * r * r / wtot;
  232. dphiz = dphiz * r / wtot;
  233. dphiphi -= phi * phi;
  234. dphiz -= phi * z;
  235. dzz /= wtot;
  236. dzz -= z * z;
  237. fLambda2 = 0.5 * (dphiphi + dzz) + TMath::Sqrt(0.25 * (dphiphi - dzz) * (dphiphi - dzz) + dphiz * dphiz);
  238. fLambda1 = 0.5 * (dphiphi + dzz) - TMath::Sqrt(0.25 * (dphiphi - dzz) * (dphiphi - dzz) + dphiz * dphiz);
  239. } else {
  240. fLambda2 = fLambda1 = 0.;
  241. }
  242. fX = x;
  243. fY = y;
  244. fZ = z;
  245. // calculates agreement of the cluster shape with parameterization in MpdEmcClusterizer::ShowerShape
  246. // and core energies
  247. fChi2 = 0;
  248. int ndf = 0;
  249. fEcore1p = 0.;
  250. fEcore2p = 0.;
  251. digIterator = fDigitIDEnergy.begin();
  252. while (digIterator != fDigitIDEnergy.end()) {
  253. int detID = (*digIterator).first;
  254. float energy = (*digIterator).second;
  255. digIterator++;
  256. double xi, yi, zi;
  257. geom->DetIdToGlobalPosition(detID, xi, yi, zi);
  258. double dz = zi - fZ;
  259. double dphi = TMath::Sqrt((xi - fX) * (xi - fX) + (yi - fY) * (yi - fY));
  260. double ss = MpdEmcClusterizerKI::ShowerShape(dphi, dz);
  261. if (ss > simParams->EcoreCut1()) {
  262. fEcore1p += energy;
  263. }
  264. if (ss > simParams->EcoreCut2()) {
  265. fEcore2p += energy;
  266. }
  267. if (ss > simParams->Chi2radiusCut()) {
  268. double sigma = 0.008 * ss * fE * (1 - ss); // TODO: What is the meaning of 0.008???
  269. if (sigma != 0) {
  270. fChi2 += pow(ss * fE - energy, 2) / sigma;
  271. ndf++;
  272. }
  273. }
  274. }
  275. if (ndf > 0) {
  276. fChi2 /= ndf;
  277. } else {
  278. fChi2 = -1;
  279. }
  280. //Correct full energy
  281. fE= simParams->ENonLinCorrection(0)+simParams->ENonLinCorrection(1)*fE ;
  282. // Prepare for writing
  283. FillArrays();
  284. }
  285. //____________________________________________________________________________
  286. int MpdEmcClusterKI::GetNumberOfLocalMax(int* maxAt, float* maxAtEnergy) const
  287. {
  288. // Calculates the number of local maxima in the cluster using LocalMaxCut as the minimum
  289. // energy difference between maximum and surrounding digits
  290. int n = GetMultiplicity();
  291. MpdEmcGeoUtils* geom = MpdEmcGeoUtils::GetInstance();
  292. MpdEmcSimParams* simParams = MpdEmcSimParams::GetInstance();
  293. float locMaxCut = simParams->LocalMaximumCut();
  294. bool* isLocalMax = new bool[n];
  295. for (int i = 0; i < n; i++)
  296. //isLocalMax[i] = true;
  297. {
  298. isLocalMax[i] = false;
  299. float en1 = fDigitIDEnergy[i].second;
  300. if (en1 > simParams->ClusteringThreshold()) isLocalMax[i] = true;
  301. }
  302. for (int i = 0; i < n; i++) {
  303. int detId1 = fDigitIDEnergy[i].first;
  304. float en1 = fDigitIDEnergy[i].second;
  305. for (int j = i + 1; j < n; j++) {
  306. int detId2 = fDigitIDEnergy[j].first;
  307. float en2 = fDigitIDEnergy[j].second;
  308. if (geom->AreNeighboursVertex(detId1, detId2) == 1) {
  309. if (en1 > en2) {
  310. isLocalMax[j] = false;
  311. // but may be digit too is not local max ?
  312. if (en2 > en1 - locMaxCut) {
  313. isLocalMax[i] = false;
  314. }
  315. } else {
  316. isLocalMax[i] = false;
  317. // but may be digitN is not local max too?
  318. if (en1 > en2 - locMaxCut) {
  319. isLocalMax[j] = false;
  320. }
  321. }
  322. } // if Areneighbours
  323. } // digit j
  324. } // digit i
  325. int iDigitN = 0;
  326. for (int i = 0; i < n; i++) {
  327. if (isLocalMax[i]) {
  328. maxAt[iDigitN] = i;
  329. maxAtEnergy[iDigitN] = fDigitIDEnergy[i].second;
  330. iDigitN++;
  331. if (iDigitN >= simParams->NLMMax()) { // Note that size of output arrays is limited:
  332. LOG(ERROR) << "Too many local maxima, cluster multiplicity " << n;
  333. return 0;
  334. }
  335. }
  336. }
  337. delete[] isLocalMax;
  338. return iDigitN;
  339. }
  340. //____________________________________________________________________________
  341. void MpdEmcClusterKI::Purify()
  342. {
  343. // Removes digits below threshold
  344. // If after purifying isolated digits remain, remove them too
  345. /*
  346. Int_t tempo[fMaxDigit];
  347. Float_t tempoE[fMaxDigit];
  348. Int_t mult = 0 ;
  349. for(Int_t iDigit=0;iDigit< fMulDigit ;iDigit++){
  350. if(fEnergyList[iDigit] > threshold){
  351. tempo[mult] = fDigitsList[iDigit] ;
  352. tempoE[mult] = fEnergyList[iDigit] ;
  353. mult++ ;
  354. }
  355. }
  356. if(mult==0){ //too soft cluster
  357. fMulDigit =0 ;
  358. fAmp = 0.; //Recalculate total energy
  359. }
  360. //Remove non-connected cells
  361. Int_t index[mult] ;
  362. Bool_t used[mult] ;
  363. for(Int_t i=0; i<mult; i++) used[i]=0 ;
  364. Int_t inClu=0 ;
  365. Double_t eMax=0. ;
  366. //find maximum
  367. for(Int_t iDigit=0; iDigit<mult; iDigit++) {
  368. if(eMax<tempoE[iDigit]){
  369. eMax=tempoE[iDigit];
  370. index[0]=iDigit ;
  371. inClu=1 ;
  372. }
  373. }
  374. if(mult>0)
  375. used[index[0]]=kTRUE ; //mark as used
  376. for(Int_t i=0; i<inClu; i++){
  377. AliPHOSDigit * digit = (AliPHOSDigit *) digits->At(tempo[index[i]]) ;
  378. for(Int_t iDigit=0 ;iDigit<mult; iDigit++){
  379. if(used[iDigit]) //already used
  380. continue ;
  381. AliPHOSDigit * digitN = (AliPHOSDigit *) digits->At(tempo[iDigit]) ;
  382. if(AreNeighbours(digit,digitN)){
  383. index[inClu]= iDigit ;
  384. inClu++ ;
  385. used[iDigit]=kTRUE ;
  386. }
  387. }
  388. }
  389. fMulDigit = inClu ;
  390. delete [] fDigitsList ;
  391. delete [] fEnergyList ;
  392. fDigitsList = new Int_t[fMulDigit];
  393. fEnergyList = new Float_t[fMulDigit];
  394. fAmp = 0.; //Recalculate total energy
  395. for(Int_t iDigit=0;iDigit< fMulDigit ;iDigit++){
  396. fDigitsList[iDigit] = tempo[index[iDigit]];
  397. fEnergyList[iDigit] = tempoE[index[iDigit]] ;
  398. fAmp+=tempoE[index[iDigit]];
  399. }
  400. */
  401. }
  402. void MpdEmcClusterKI::FillArrays()
  403. {
  404. // As root is not able to handle std::map and std::vector, Fill arrays to store to disk
  405. fNDigits = fDigitIDEnergy.size();
  406. fDitisId = new Int_t[fNDigits];
  407. fDigitsE = new Float_t[fNDigits];
  408. std::vector<std::pair<int, float>>::iterator digIterator = fDigitIDEnergy.begin();
  409. int i = 0;
  410. while (digIterator != fDigitIDEnergy.end()) {
  411. fDitisId[i] = (*digIterator).first;
  412. fDigitsE[i] = (*digIterator).second;
  413. digIterator++;
  414. i++;
  415. }
  416. // Sort map accordig to deposited energy
  417. std::vector<std::pair<int, float>> vec(fMCTracks.begin(), fMCTracks.end());
  418. std::sort(vec.begin(), vec.end(), compE);
  419. fNPrimaries = fMCTracks.size();
  420. fNPrimaries = TMath::Min(fNPrimaries, MpdEmcSimParams::GetInstance()->NPrimMax());
  421. fPrimId = new Int_t[fNPrimaries];
  422. fPrimE = new Float_t[fNPrimaries];
  423. i = 0;
  424. std::vector<std::pair<int, float>>::iterator primIterator = vec.begin();
  425. while (primIterator != vec.end() && i < fNPrimaries) {
  426. fPrimId[i] = (*primIterator).first;
  427. fPrimE[i] = (*primIterator).second;
  428. i++;
  429. primIterator++;
  430. }
  431. }
  432. void MpdEmcClusterKI::CorrectVertex(double zVtx)
  433. {
  434. // correct cluster position for Z ccordinate of vertex
  435. // function is sA*sin(z/sW)+a+b*z, where sA,sW,a,b may depend on Zvtx and E
  436. MpdEmcSimParams* simParams = MpdEmcSimParams::GetInstance();
  437. double logE = TMath::Log(fE); // NB: here we assume that energy already corrected for 1/3
  438. double sA = simParams->ZcorrSinA(0) + simParams->ZcorrSinA(1) * logE; //-5.59366e-01 -2.67599e-02*logE
  439. double sW = simParams->ZcorrSinW(0) + simParams->ZcorrSinW(1) * logE;
  440. double a = (simParams->ZcorrA(0) + simParams->ZcorrA(1) * logE) * zVtx;
  441. double b = simParams->ZcorrB(0) + simParams->ZcorrB(1) * logE;
  442. fZ += (sA * TMath::Sin(fZ / sW) + a + b * fZ);
  443. }
  444. ClassImp(MpdEmcClusterKI)