MpdEmcClusterCreation.cxx 9.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295
  1. ////////////////////////////////////////////////////////////////
  2. // //
  3. // MpdEmcClusterCreation //
  4. // EMC clsuters in MpdEmcCluster, v02 //
  5. // Author List : Martemianov M., 2019 //
  6. // //
  7. ////////////////////////////////////////////////////////////////
  8. #include "TClonesArray.h"
  9. #include "FairRootManager.h"
  10. #include "FairRunAna.h"
  11. #include "FairRunSim.h"
  12. #include "FairEventHeader.h"
  13. #include "MpdEmcHit.h"
  14. #include "MpdEmcCluster.h"
  15. #include "MpdEmcClusterCreation.h"
  16. #include "MpdMCTrack.h"
  17. #include "TGeoManager.h"
  18. namespace EMC {
  19. UInt_t nRows = 300; // number of rows (xy - plane)
  20. UInt_t nMods = 128; // number of modules in row
  21. Float_t Distance3DCut = 20.0; // cm
  22. Float_t HitThreshold = 5.0; // MeV
  23. UInt_t algorithmNum = 2; // algorithm number
  24. UInt_t frameRow = 5; // row frame
  25. UInt_t frameMod = 5; // module frame
  26. }
  27. using namespace std;
  28. using namespace TMath;
  29. // ----- Default constructor -------------------------------------------
  30. MpdEmcClusterCreation::MpdEmcClusterCreation() :
  31. FairTask("EMC cluster creation"), fEnergyThreshold(EMC::HitThreshold), algoIndex(EMC::algorithmNum),
  32. fMaxClusterRadius(EMC::Distance3DCut), rowFrame(EMC::frameRow), modFrame(EMC::frameMod){
  33. }
  34. // ----- Destructor ----------------------------------------------------
  35. MpdEmcClusterCreation::~MpdEmcClusterCreation() {
  36. }
  37. // -------------------------------------------------------------------------
  38. // ----- Public method Init --------------------------------------------
  39. InitStatus MpdEmcClusterCreation::Init() {
  40. cout << "******************* EMC INIT *********************" << endl;
  41. // Get RootManager
  42. FairRootManager* ioman = FairRootManager::Instance();
  43. if (!ioman) {
  44. cout << "-E- MpdEmcClusterCreation::Init: "
  45. << "RootManager not instantiated!" << endl;
  46. return kFATAL;
  47. }
  48. // Get input array
  49. fHitArray = (TClonesArray*) ioman->GetObject("MpdEmcHit");
  50. if (!fHitArray) {
  51. cout << "-W- MpdEmcHitCreation::Init: " << "No MpdEmcHit array!" << endl;
  52. return kERROR;
  53. }
  54. fMcTrackArray = (TClonesArray*) ioman->GetObject("MCTrack");
  55. if (!fMcTrackArray) {
  56. cout << "-W- MpdEmcClusterCreation::Init: " << "No MCTrack array!" << endl;
  57. return kERROR;
  58. }
  59. // Create and register output array
  60. fClusterArray = new TClonesArray("MpdEmcCluster", 100);
  61. ioman->Register("MpdEmcCluster", "EMC", fClusterArray, kTRUE);
  62. ioman->Register("MCTrack","EMC",fMcTrackArray, kTRUE);
  63. //// ioman->Register("MpdEmcHit","EMC",fHitArray, kTRUE);
  64. cout << "-I- MpdEmcClusterCreation: Intialization successfull" << endl;
  65. return kSUCCESS;
  66. }
  67. void MpdEmcClusterCreation::Finish() {
  68. cout << "\n-I- MpdEmcClusterCreation: Finish" << endl;
  69. }
  70. void MpdEmcClusterCreation::Exec(Option_t* opt) {
  71. cout << "\n-I- MpdEmcClusterCreation: Event No. " << FairRun::Instance()->GetEventHeader()->GetMCEntryNumber() << endl;
  72. // Reset output Array
  73. if (!fClusterArray) Fatal("MpdEmcClusterCreation::Exec)", "No array of clusters");
  74. fClusterArray->Delete();
  75. UInt_t nHits = fHitArray->GetEntriesFast();
  76. MpdEmcCluster *cluster = NULL;
  77. Bool_t addHit, addCut;
  78. vector<UInt_t> modId, relHits, supHits, outHits;
  79. vector<vector<UInt_t> > relCluster;
  80. vector<Float_t> xHit, yHit, zHit;
  81. vector<Float_t> energy, time;
  82. if (nHits > 0) {
  83. for (UInt_t iHit = 0; iHit < nHits; ++iHit) {
  84. MpdEmcHit* hit = (MpdEmcHit*)fHitArray->At(iHit);
  85. if (hit->GetE()*1000. > fEnergyThreshold) {
  86. modId.push_back(1000*(hit->GetRow()+1)+(hit->GetMod()+1));
  87. xHit.push_back(hit->GetRhoCenter()*cos(hit->GetPhiCenter()*DegToRad()));
  88. yHit.push_back(hit->GetRhoCenter()*sin(hit->GetPhiCenter()*DegToRad()));
  89. zHit.push_back(hit->GetZCenter());
  90. energy.push_back(hit->GetE());
  91. time.push_back(hit->GetTime());
  92. }
  93. }
  94. UInt_t row, mod;
  95. vector<Float_t> iterEnergy = energy;
  96. while (iterEnergy.size()) {
  97. relHits.clear(); outHits.clear();
  98. UInt_t indexMax = distance(iterEnergy.begin(), max_element(iterEnergy.begin(), iterEnergy.end()));
  99. if (iterEnergy.size() == 1) indexMax = 0;
  100. indexMax = distance(energy.begin(),find(energy.begin(),energy.end(),iterEnergy[indexMax]));
  101. relHits.push_back(indexMax);
  102. row = modId[indexMax]/1000; mod = modId[indexMax] - row*1000;
  103. SearchFrameHits(row, mod, outHits);
  104. for (UInt_t iHit = 0; iHit < energy.size(); iHit++) {
  105. Float_t dist = sqrt(pow(xHit[indexMax] - xHit[iHit],2) + pow(yHit[indexMax] - yHit[iHit],2) +
  106. pow(zHit[indexMax] - zHit[iHit],2));
  107. addHit = kTRUE; addCut = kFALSE;
  108. for (UInt_t iAdd = 0; iAdd < supHits.size(); iAdd++)
  109. if (supHits[iAdd] == iHit) addHit = kFALSE;
  110. if (algoIndex == 1) addCut = (dist < fMaxClusterRadius);
  111. if (algoIndex == 2) {
  112. for (UInt_t iRel = 0; iRel < outHits.size(); iRel++)
  113. if (modId[iHit] == outHits[iRel]) addCut = kTRUE;
  114. }
  115. if ( (addCut) && (addHit) ) {
  116. UInt_t pos = distance(iterEnergy.begin(),find(iterEnergy.begin(),iterEnergy.end(),energy[iHit]));
  117. if (iterEnergy.size() == 1) pos = 0;
  118. iterEnergy.erase(iterEnergy.begin()+pos);
  119. if (indexMax != iHit) relHits.insert(relHits.end(), iHit);
  120. }
  121. }
  122. for (UInt_t iRel = 0; iRel < relHits.size(); iRel++)
  123. supHits.insert(supHits.end(), relHits[iRel]);
  124. relCluster.push_back(relHits);
  125. }
  126. }
  127. // Fill cluster information (0,0 - central hit)
  128. TVector3 fPos;
  129. UInt_t numHits;
  130. Float_t rho, phi, theta, rhoCenter{0}, phiCenter{0}; //thetaCenter,
  131. Float_t thetaPos, phiPos{0}, rhoPos, zPos, fEnergy, fTime;
  132. vector<Float_t> eH, xH, yH, zH;
  133. for (UInt_t iCluster = 0; iCluster < relCluster.size(); iCluster++) {
  134. fEnergy = 0.; fTime = 0.; phiPos = 0.; thetaPos = 0.;
  135. numHits = relCluster[iCluster].size();
  136. eH.clear(); xH.clear(); yH.clear(); zH.clear();
  137. for (UInt_t iHit = 0; iHit < relCluster[iCluster].size(); iHit++) {
  138. eH.insert(eH.end(), energy[relCluster[iCluster][iHit]]);
  139. xH.insert(xH.end(),xHit[relCluster[iCluster][iHit]]);
  140. yH.insert(yH.end(),yHit[relCluster[iCluster][iHit]]);
  141. zH.insert(zH.end(),zHit[relCluster[iCluster][iHit]]);
  142. rho = sqrt(xH[iHit]*xH[iHit]+yH[iHit]*yH[iHit]+zH[iHit]*zH[iHit]);
  143. theta = ACos(zH[iHit]/rho);
  144. phi = ATan2(yH[iHit],xH[iHit])*RadToDeg();
  145. if (phi < 0) phi += 360;
  146. if (iHit == 0) {
  147. phiCenter = phi; rhoCenter = rho; //thetaCenter = theta;
  148. }
  149. fEnergy += eH[iHit]; thetaPos += theta*eH[iHit];
  150. phiPos += MpdEmcMath::GetPhiDiff(phi, phiCenter)*eH[iHit];
  151. fTime += time[relCluster[iCluster][iHit]]*eH[iHit];
  152. }
  153. thetaPos = thetaPos/fEnergy;
  154. phiPos = phiCenter + phiPos/fEnergy;
  155. fTime = fTime/fEnergy;
  156. rhoPos = rhoCenter*sin(thetaPos);
  157. zPos = rhoCenter*cos(thetaPos);
  158. fPos.SetXYZ(phiPos, rhoPos, zPos);
  159. cluster = new((*fClusterArray)[fClusterArray->GetEntriesFast()]) MpdEmcCluster(fEnergy, fTime, fPos);
  160. cluster->SetNHits(numHits);
  161. cluster->SetRad(cluster->ComputeClusterRadius(fPos, xH, yH, zH, eH));
  162. }
  163. }
  164. // Search relative hits in the definite frame
  165. void MpdEmcClusterCreation::SearchFrameHits(UInt_t row, UInt_t mod, vector<UInt_t> &outMod) {
  166. UInt_t ind, curRow;
  167. UInt_t nRow = EMC::nRows;
  168. UInt_t nMod = EMC::nMods;
  169. vector<Int_t> geoMod;
  170. Int_t dRowMin = row - rowFrame, dRowMax = row + rowFrame;
  171. Int_t dModMin = mod - modFrame, dModMax = mod + modFrame;
  172. if (dRowMin <= 0) dRowMin = nRow - fabs(dRowMin);
  173. if ((UInt_t)dRowMax >= nRow) dRowMax = dRowMax - nRow;
  174. if (dModMin <= 0) dModMin = 1;
  175. if ((UInt_t)dModMax >= nMod) dModMax = nMod;
  176. curRow = dRowMin;
  177. for (UInt_t iGeo1 = 0; iGeo1 < 2*rowFrame + 1; iGeo1++) {
  178. if (curRow == nRow + 1) { curRow = 1;}
  179. for (UInt_t iGeo2 = dModMin; iGeo2 < (UInt_t)dModMax + 1; iGeo2++)
  180. geoMod.insert(geoMod.end(), 1000*curRow + iGeo2);
  181. curRow++;
  182. }
  183. outMod.clear();
  184. for (UInt_t iHit = 0; iHit < (UInt_t)fHitArray->GetEntriesFast(); iHit++) {
  185. MpdEmcHit* hit = (MpdEmcHit*)fHitArray->At(iHit);
  186. ind = 1000*(hit->GetRow()+1) + (hit->GetMod()+1);
  187. for (UInt_t iGeo = 0; iGeo < geoMod.size(); iGeo++)
  188. if (ind == (UInt_t) geoMod[iGeo]) outMod.push_back(ind);
  189. }
  190. }
  191. // Search relative hits
  192. void MpdEmcClusterCreation::SearchRelativeHits(UInt_t row, UInt_t mod, vector<Int_t> &outMod) {
  193. UInt_t ind;
  194. UInt_t nRow = EMC::nRows;
  195. UInt_t nMod = EMC::nMods;
  196. vector<Int_t> geoMod;
  197. row = row + 1; mod = mod + 1;
  198. UInt_t rowPrev = row-1, rowNext = row+1;
  199. if (row == 1) rowPrev = nRow;
  200. if (row == nRow) rowNext = 1;
  201. geoMod = MpdEmcMath::make_vector<int>() <<
  202. 1000*rowPrev + (mod-1) << 1000*rowPrev + (mod) << 1000*rowPrev + (mod+1) <<
  203. 1000*row + (mod-1) << 1000*row + (mod+1) << 1000*rowNext + (mod-1) <<
  204. 1000*rowNext + (mod) << 1000*rowNext + (mod+1);
  205. if (mod - 1 == 0) geoMod.erase(geoMod.begin()+1,geoMod.begin()+3);
  206. if (mod + 1 > nMod) geoMod.erase(geoMod.begin()+5,geoMod.begin()+7);
  207. outMod.clear();
  208. outMod.insert(outMod.begin(), 1000*(row) + mod);
  209. for (UInt_t iHit = 0; iHit < (UInt_t) fHitArray->GetEntriesFast(); iHit++) {
  210. MpdEmcHit* hit = (MpdEmcHit*)fHitArray->At(iHit);
  211. ind = 1000*(hit->GetRow()+1) + (hit->GetMod()+1);
  212. for (UInt_t iGeo = 0; iGeo < geoMod.size(); iGeo++)
  213. if (ind == (UInt_t) geoMod[iGeo]) outMod.push_back(ind);
  214. }
  215. }
  216. // -------------------------------------------------------------------------
  217. ClassImp(MpdEmcClusterCreation)