MpdFemtoLikeSignAnalysis.cxx 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319
  1. //
  2. // The class is a base class for like-sign analysis with z-binning
  3. //
  4. // MpdFemtoMaker headers
  5. // Base
  6. #include "MpdFemtoBaseTrackCut.h"
  7. #include "MpdFemtoBaseV0Cut.h"
  8. // Infrastructure
  9. #include "MpdFemtoLikeSignAnalysis.h"
  10. #include "MpdFemtoParticleCollection.h"
  11. #include "MpdFemtoPicoEventCollectionVector.h"
  12. #include "MpdFemtoPicoEventCollectionVectorHideAway.h"
  13. ClassImp(MpdFemtoLikeSignAnalysis)
  14. /// This little function used to apply ParticleCuts (TrackCuts or V0Cuts)
  15. /// and fill ParticleCollections of picoEvent it is called
  16. /// from MpdFemtoAnalysis::ProcessEvent()
  17. extern void fillHbtParticleCollection(MpdFemtoBaseParticleCut* partCut,
  18. MpdFemtoEvent* hbtEvent,
  19. MpdFemtoParticleCollection* partCollection, char doReshuffle);
  20. //_________________
  21. MpdFemtoLikeSignAnalysis::MpdFemtoLikeSignAnalysis(const unsigned int& bins, const double& min,
  22. const double& max) : MpdFemtoAnalysis() {
  23. // Constructor
  24. mVertexBins = bins;
  25. mVertexZ[0] = min;
  26. mVertexZ[1] = max;
  27. mUnderFlow = 0;
  28. mOverFlow = 0;
  29. if (mMixingBuffer) delete mMixingBuffer;
  30. mPicoEventCollectionVectorHideAway =
  31. new MpdFemtoPicoEventCollectionVectorHideAway(mVertexBins, mVertexZ[0], mVertexZ[1]);
  32. }
  33. //_________________
  34. MpdFemtoLikeSignAnalysis::MpdFemtoLikeSignAnalysis(const MpdFemtoLikeSignAnalysis& a) : MpdFemtoAnalysis(a) {
  35. // Copy constructor
  36. mVertexBins = a.mVertexBins;
  37. mVertexZ[0] = a.mVertexZ[0];
  38. mVertexZ[1] = a.mVertexZ[1];
  39. mUnderFlow = 0;
  40. mOverFlow = 0;
  41. if (mMixingBuffer) delete mMixingBuffer;
  42. mPicoEventCollectionVectorHideAway =
  43. new MpdFemtoPicoEventCollectionVectorHideAway(mVertexBins, mVertexZ[0], mVertexZ[1]);
  44. }
  45. //_________________
  46. MpdFemtoLikeSignAnalysis& MpdFemtoLikeSignAnalysis::operator=(const MpdFemtoLikeSignAnalysis& a) {
  47. // Assignement operator
  48. if (this != &a) {
  49. mVertexBins = a.mVertexBins;
  50. mVertexZ[0] = a.mVertexZ[0];
  51. mVertexZ[1] = a.mVertexZ[1];
  52. mUnderFlow = 0;
  53. mOverFlow = 0;
  54. if (mMixingBuffer) delete mMixingBuffer;
  55. if (mPicoEventCollectionVectorHideAway) delete mPicoEventCollectionVectorHideAway;
  56. mPicoEventCollectionVectorHideAway =
  57. new MpdFemtoPicoEventCollectionVectorHideAway(mVertexBins, mVertexZ[0], mVertexZ[1]);
  58. }
  59. return *this;
  60. }
  61. //_________________
  62. MpdFemtoLikeSignAnalysis::~MpdFemtoLikeSignAnalysis() {
  63. delete mPicoEventCollectionVectorHideAway;
  64. mPicoEventCollectionVectorHideAway = nullptr;
  65. }
  66. //_________________
  67. MpdFemtoString MpdFemtoLikeSignAnalysis::report() {
  68. char Ctemp[200];
  69. std::cout << "MpdFemtoLikeSignAnalysis - constructing Report..." << std::endl;
  70. MpdFemtoString temp = "-----------\nHbt Analysis Report:\n";
  71. sprintf(Ctemp, "Events are mixed in %d bins in the range %E cm to %E cm.\n", mVertexBins, mVertexZ[0], mVertexZ[1]);
  72. temp += Ctemp;
  73. sprintf(Ctemp, "Events underflowing: %d\n", mUnderFlow);
  74. temp += Ctemp;
  75. sprintf(Ctemp, "Events overflowing: %d\n", mOverFlow);
  76. temp += Ctemp;
  77. sprintf(Ctemp, "Now adding MpdFemtoAnalysis(base) Report\n");
  78. temp += Ctemp;
  79. temp += "Adding MpdFemtoAnalysis(base) Report now:\n";
  80. temp += MpdFemtoAnalysis::report();
  81. temp += "-------------\n";
  82. MpdFemtoString returnThis = temp;
  83. return returnThis;
  84. }
  85. //_________________________
  86. void MpdFemtoLikeSignAnalysis::processEvent(const MpdFemtoEvent* hbtEvent) {
  87. // Perform all the analysis tasks for a single event
  88. // get right mixing buffer
  89. double vertexZ = hbtEvent->primaryVertex().Z();
  90. mMixingBuffer = mPicoEventCollectionVectorHideAway->picoEventCollection(vertexZ);
  91. if (!mMixingBuffer) {
  92. if (vertexZ < mVertexZ[0]) mUnderFlow++;
  93. if (vertexZ > mVertexZ[1]) mOverFlow++;
  94. return;
  95. }
  96. // Startup for EbyE
  97. eventBegin(hbtEvent);
  98. // Event cut and event cut monitor
  99. bool tmpPassEvent = mEventCut->pass(hbtEvent);
  100. mEventCut->fillCutMonitor(hbtEvent, tmpPassEvent);
  101. if (tmpPassEvent) {
  102. mNeventsProcessed++;
  103. std::cout << "MpdFemtoLikeSignAnalysis::processEvent() - " << hbtEvent->trackCollection()->size();
  104. std::cout << " #track=" << hbtEvent->trackCollection()->size();
  105. // OK, analysis likes the event-- build a pico event from it,
  106. // using tracks the analysis likes...
  107. // This is what we will make pairs from and put in Mixing Buffer
  108. MpdFemtoPicoEvent* picoEvent = new MpdFemtoPicoEvent;
  109. fillHbtParticleCollection(mFirstParticleCut, (MpdFemtoEvent*) hbtEvent, picoEvent->firstParticleCollection(), mReshuffle);
  110. if (!(analyzeIdenticalParticles())) {
  111. fillHbtParticleCollection(mSecondParticleCut, (MpdFemtoEvent*) hbtEvent, picoEvent->secondParticleCollection(), mReshuffle);
  112. }
  113. std::cout << " #particles in First, Second Collections: "
  114. << picoEvent->firstParticleCollection()->size() << " "
  115. << picoEvent->secondParticleCollection()->size() << std::endl;
  116. if (picoEvent->secondParticleCollection()->size() *
  117. picoEvent->firstParticleCollection()->size() == 0) {
  118. delete picoEvent;
  119. std::cout << "MpdFemtoLikeSignAnalysis - picoEvent deleted due to empty collection " << std::endl;
  120. return;
  121. }
  122. // OK, pico event is built
  123. // make real pairs...
  124. // Fabrice points out that we do not need to keep creating/deleting pairs all the time
  125. // We only ever need ONE pair, and we can just keep changing internal pointers
  126. // this should help speed things up
  127. MpdFemtoPair* ThePair = new MpdFemtoPair;
  128. MpdFemtoParticleIterator PartIter1;
  129. MpdFemtoParticleIterator PartIter2;
  130. MpdFemtoCorrFctnIterator CorrFctnIter;
  131. // Always
  132. MpdFemtoParticleIterator StartOuterLoop = picoEvent->firstParticleCollection()->begin();
  133. // Will be one less if identical
  134. MpdFemtoParticleIterator EndOuterLoop = picoEvent->firstParticleCollection()->end();
  135. MpdFemtoParticleIterator StartInnerLoop;
  136. MpdFemtoParticleIterator EndInnerLoop;
  137. // Only use First collection
  138. if (analyzeIdenticalParticles()) {
  139. // Iuter loop goes to next-to-last particle in First collection
  140. EndOuterLoop--;
  141. // Inner loop goes to last particle in First collection
  142. EndInnerLoop = picoEvent->firstParticleCollection()->end();
  143. } else {
  144. // Non-identical - loop over First and Second collections
  145. // Inner loop starts at first particle in Second collection
  146. StartInnerLoop = picoEvent->secondParticleCollection()->begin();
  147. // Inner loop goes to last particle in Second collection
  148. EndInnerLoop = picoEvent->secondParticleCollection()->end();
  149. }
  150. // real pairs
  151. for (PartIter1 = StartOuterLoop; PartIter1 != EndOuterLoop; PartIter1++) {
  152. if (analyzeIdenticalParticles()) {
  153. StartInnerLoop = PartIter1;
  154. StartInnerLoop++;
  155. }
  156. ThePair->setTrack1(*PartIter1);
  157. for (PartIter2 = StartInnerLoop; PartIter2 != EndInnerLoop; PartIter2++) {
  158. ThePair->setTrack2(*PartIter2);
  159. // The following lines have to be uncommented if you want pairCutMonitors
  160. // they are not in for speed reasons
  161. // bool tmpPassPair = mPairCut->Pass(ThePair);
  162. // mPairCut->FillCutMonitor(ThePair, tmpPassPair);
  163. // if ( tmpPassPair ) {
  164. if (mPairCut->pass(ThePair)) {
  165. for (CorrFctnIter = mCorrFctnCollection->begin();
  166. CorrFctnIter != mCorrFctnCollection->end(); CorrFctnIter++) {
  167. MpdFemtoBaseLikeSignCorrFctn* CorrFctn = dynamic_cast<MpdFemtoBaseLikeSignCorrFctn*> (*CorrFctnIter);
  168. if (CorrFctn) {
  169. CorrFctn->addRealPair(ThePair);
  170. }
  171. }
  172. } //if ( mPairCut->pass(ThePair) )
  173. } //for (PartIter2 = StartInnerLoop; PartIter2!=EndInnerLoop; PartIter2++)
  174. } //for ( PartIter1=StartOuterLoop; PartIter1!=EndOuterLoop; PartIter1++ )
  175. #ifdef STHBTDEBUG
  176. std::cout << "MpdFemtoLikeSignAnalysis::ProcessEvent() - reals done" << std::endl;
  177. #endif
  178. MpdFemtoParticleIterator nextIter;
  179. MpdFemtoParticleIterator prevIter;
  180. // Like-sign first partilce collection pairs
  181. prevIter = EndOuterLoop;
  182. prevIter--;
  183. for (PartIter1 = StartOuterLoop; PartIter1 != prevIter; PartIter1++) {
  184. ThePair->setTrack1(*PartIter1);
  185. nextIter = PartIter1;
  186. nextIter++;
  187. for (PartIter2 = nextIter; PartIter2 != EndOuterLoop; PartIter2++) {
  188. ThePair->setTrack2(*PartIter2);
  189. // The following lines have to be uncommented if you want pairCutMonitors
  190. // they are not in for speed reasons
  191. // bool tmpPassPair = mPairCut->Pass(ThePair);
  192. // mPairCut->FillCutMonitor(ThePair, tmpPassPair);
  193. // if ( tmpPassPair ) {
  194. if (mPairCut->pass(ThePair)) {
  195. for (CorrFctnIter = mCorrFctnCollection->begin();
  196. CorrFctnIter != mCorrFctnCollection->end(); CorrFctnIter++) {
  197. MpdFemtoBaseLikeSignCorrFctn* CorrFctn = dynamic_cast<MpdFemtoBaseLikeSignCorrFctn*> (*CorrFctnIter);
  198. if (CorrFctn) {
  199. CorrFctn->addLikeSignPositivePair(ThePair);
  200. }
  201. }
  202. } //if ( mPairCut->pass(ThePair) )
  203. } //for ( PartIter2 = nextIter; PartIter2!=EndOuterLoop; PartIter2++)
  204. } //for ( PartIter1=StartOuterLoop; PartIter1!=prevIter; PartIter1++)
  205. #ifdef STHBTDEBUG
  206. std::cout << "MpdFemtoLikeSignAnalysis::processEvent() - like sign first collection done" << std::endl;
  207. #endif
  208. // Like-sign second partilce collection pairs
  209. prevIter = EndInnerLoop;
  210. prevIter--;
  211. for (PartIter1 = StartInnerLoop; PartIter1 != prevIter; PartIter1++) {
  212. ThePair->setTrack1(*PartIter1);
  213. nextIter = PartIter1;
  214. nextIter++;
  215. for (PartIter2 = nextIter; PartIter2 != EndInnerLoop; PartIter2++) {
  216. ThePair->setTrack2(*PartIter2);
  217. // The following lines have to be uncommented if you want pairCutMonitors
  218. // they are not in for speed reasons
  219. // bool tmpPassPair = mPairCut->Pass(ThePair);
  220. // mPairCut->FillCutMonitor(ThePair, tmpPassPair);
  221. // if ( tmpPassPair ) {
  222. if (mPairCut->pass(ThePair)) {
  223. for (CorrFctnIter = mCorrFctnCollection->begin();
  224. CorrFctnIter != mCorrFctnCollection->end(); CorrFctnIter++) {
  225. MpdFemtoBaseLikeSignCorrFctn* CorrFctn = dynamic_cast<MpdFemtoBaseLikeSignCorrFctn*> (*CorrFctnIter);
  226. if (CorrFctn) {
  227. CorrFctn->addLikeSignNegativePair(ThePair);
  228. }
  229. }
  230. } //if ( mPairCut->pass(ThePair) )
  231. } //for (PartIter2 = nextIter; PartIter2!=EndInnerLoop; PartIter2++)
  232. //for ( PartIter1=StartInnerLoop; PartIter1!=prevIter; PartIter1++)
  233. #ifdef STHBTDEBUG
  234. std::cout << "MpdFemtoLikeSignAnalysis::processEvent() - like sign second collection done" << std::endl;
  235. #endif
  236. if (mixingBufferFull()) {
  237. #ifdef STHBTDEBUG
  238. std::cout << "Mixing Buffer is full - lets rock and roll" << std::endl;
  239. #endif
  240. } else {
  241. std::cout << "Mixing Buffer not full -gotta wait " << mixingBuffer()->size() << std::endl;
  242. }
  243. if (mixingBufferFull()) {
  244. StartOuterLoop = picoEvent->firstParticleCollection()->begin();
  245. EndOuterLoop = picoEvent->firstParticleCollection()->end();
  246. MpdFemtoPicoEvent* storedEvent;
  247. MpdFemtoPicoEventIterator picoEventIter;
  248. for (picoEventIter = mixingBuffer()->begin();
  249. picoEventIter != mixingBuffer()->end(); picoEventIter++) {
  250. storedEvent = *picoEventIter;
  251. if (analyzeIdenticalParticles()) {
  252. StartInnerLoop = storedEvent->firstParticleCollection()->begin();
  253. EndInnerLoop = storedEvent->firstParticleCollection()->end();
  254. } else {
  255. StartInnerLoop = storedEvent->secondParticleCollection()->begin();
  256. EndInnerLoop = storedEvent->secondParticleCollection()->end();
  257. }
  258. for (PartIter1 = StartOuterLoop; PartIter1 != EndOuterLoop; PartIter1++) {
  259. ThePair->setTrack1(*PartIter1);
  260. for (PartIter2 = StartInnerLoop; PartIter2 != EndInnerLoop; PartIter2++) {
  261. ThePair->setTrack2(*PartIter2);
  262. // testing... cout << "ThePair defined... going to pair cut... ";
  263. if (mPairCut->pass(ThePair)) {
  264. // testing... cout << " ThePair passed PairCut... ";
  265. for (CorrFctnIter = mCorrFctnCollection->begin();
  266. CorrFctnIter != mCorrFctnCollection->end(); CorrFctnIter++) {
  267. MpdFemtoBaseLikeSignCorrFctn* CorrFctn = dynamic_cast<MpdFemtoBaseLikeSignCorrFctn*> (*CorrFctnIter);
  268. if (CorrFctn) {
  269. CorrFctn->addMixedPair(ThePair);
  270. //cout << " ThePair has been added to MixedPair method " << endl;
  271. } //if (CorrFctn)
  272. } // for ( CorrFctnIter=mCorrFctnCollection->begin();
  273. } //if ( mPairCut->pass(ThePair) )
  274. } //for (PartIter2=StartInnerLoop; PartIter2!=EndInnerLoop; PartIter2++)
  275. } //for ( PartIter1=StartOuterLoop; PartIter1!=EndOuterLoop; PartIter1++)
  276. } //for ( picoEventIter=mixingBuffer()->begin(); picoEventIter!=mixingBuffer()->end(); picoEventIter++)
  277. // Now get rid of oldest stored pico-event in buffer.
  278. // This means (1) delete the event from memory, (2) "pop" the pointer to it from the MixingBuffer
  279. delete mixingBuffer()->back();
  280. mixingBuffer()->pop_back();
  281. } //if ( mixingBufferFull() )
  282. delete ThePair;
  283. // Store the current pico-event in buffer
  284. mixingBuffer()->push_front(picoEvent);
  285. } //if (tmpPassEvent)
  286. // Cleanup for EbyE
  287. eventEnd(hbtEvent);
  288. // cout << "MpdFemtoLikeSignAnalysis::ProcessEvent() - return to caller ... " << endl;
  289. }
  290. }