MpdMiniDstFillTask.cxx 60 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523
  1. // MpdRoot headers
  2. #include "MpdEvent.h"
  3. #include "MpdVertex.h"
  4. #include "MpdTrack.h"
  5. #include "MpdTpcKalmanTrack.h"
  6. #include "MpdKalmanFilter.h"
  7. #include "MpdKalmanGeoScheme.h"
  8. #include "MpdTofHit.h"
  9. #include "MpdEmcDigitKI.h"
  10. #include "MpdEmcClusterKI.h"
  11. #include "MpdTofMatchingData.h"
  12. #include "MpdZdcDigi.h"
  13. #include "MpdMCTrack.h"
  14. #include "MpdGenTrack.h"
  15. // MiniDst headers
  16. #include "MpdMiniDstFillTask.h"
  17. #include "MpdMiniDst.h"
  18. #include "MpdMiniEvent.h"
  19. #include "MpdMiniTrack.h"
  20. #include "MpdMiniBTofHit.h"
  21. #include "MpdMiniBTofPidTraits.h"
  22. #include "MpdMiniBECalCluster.h"
  23. #include "MpdMiniTrackCovMatrix.h"
  24. #include "MpdMiniFHCalHit.h"
  25. #include "MpdMiniMcEvent.h"
  26. #include "MpdMiniMcTrack.h"
  27. #include "MpdMiniMessMgr.h"
  28. #include "MpdMiniArrays.cxx"
  29. //_________________
  30. MpdMiniDstFillTask::MpdMiniDstFillTask() {
  31. // Default constructor
  32. /* empty */
  33. }
  34. //_________________
  35. MpdMiniDstFillTask::MpdMiniDstFillTask(TString name) :
  36. fIsUseCovMatrix(kTRUE),
  37. fIsUseECal(kTRUE),
  38. fEventHeaders(nullptr),
  39. fEvents(nullptr),
  40. fVertices(nullptr),
  41. fTpcTracks(nullptr),
  42. fTofHits(nullptr),
  43. fTofMatching(nullptr),
  44. fMCTracks(nullptr),
  45. fGenTracks(nullptr),
  46. fEmcClusters(nullptr),
  47. fZdcDigits(nullptr),
  48. fMiniDst(new MpdMiniDst()),
  49. fBField(0.5), // mag. field in T (MUST be changed to the real magnetic field
  50. fOutputFile(nullptr),
  51. fTTree(nullptr),
  52. fSplit(99),
  53. fCompression(9),
  54. fBufferSize(65536 * 4),
  55. fMiniArrays(nullptr),
  56. fNSigmaDedxEstimator(0) {
  57. // Standard constructor
  58. TString truncatedInFileName = name;
  59. Ssiz_t lastOccurence = name.Last('/');
  60. if (lastOccurence != kNPOS) {
  61. truncatedInFileName = name((lastOccurence + 1), name.Length());
  62. }
  63. fOutputFileName = truncatedInFileName.ReplaceAll(".root", ".MiniDst.root");
  64. if (!fMcTrk2MiniMcTrk.empty()) fMcTrk2MiniMcTrk.clear();
  65. if (!fMcTrk2EcalCluster.empty()) fMcTrk2EcalCluster.clear();
  66. streamerOff();
  67. createArrays();
  68. }
  69. //_________________
  70. MpdMiniDstFillTask::~MpdMiniDstFillTask() {
  71. // Destructor
  72. if (!fMcTrk2MiniMcTrk.empty()) fMcTrk2MiniMcTrk.clear();
  73. if (!fMcTrk2EcalCluster.empty()) fMcTrk2EcalCluster.clear();
  74. delete fMiniDst;
  75. }
  76. //_________________
  77. InitStatus MpdMiniDstFillTask::Init() {
  78. // Output name must exist
  79. if (fOutputFileName.IsNull()) {
  80. return kERROR;
  81. }
  82. // Creating output tree with miniDST
  83. fOutputFile = new TFile(fOutputFileName.Data(), "RECREATE");
  84. // Inform about the creation
  85. LOG_INFO << " Output file: " << fOutputFileName.Data() << " created." << endm;
  86. // Set compression level
  87. fOutputFile->SetCompressionLevel(fCompression);
  88. int bufsize = fBufferSize;
  89. if (fSplit) {
  90. bufsize /= 4;
  91. }
  92. // Create TTree
  93. fTTree = new TTree("MiniDst", "MpdMiniDst", fSplit);
  94. fTTree->SetAutoSave(1000000);
  95. // Create arrays
  96. for (Int_t i = 0; i < MpdMiniArrays::NAllMiniArrays; ++i) {
  97. fTTree->Branch(MpdMiniArrays::miniArrayNames[i], &fMiniArrays[i], bufsize, fSplit);
  98. }
  99. // Initialize FairRoot manager
  100. FairRootManager* ioman = FairRootManager::Instance();
  101. if (!ioman) {
  102. std::cout << "[ERROR] MpdMiniDstFillTask::Init - Not FairRootManager instance has been found"
  103. << std::endl;
  104. exit(0);
  105. }
  106. // Reading all possible input branches ...
  107. // Get MC event header
  108. fEventHeaders = (FairMCEventHeader*) ioman->GetObject("MCEventHeader.");
  109. if (!fEventHeaders) {
  110. std::cout << "[WARNING] MpdMiniDstFillTask::Init - No MCEventHeader has been found"
  111. << std::endl;
  112. }
  113. // Get information on event
  114. fEvents = (MpdEvent*) ioman->GetObject("MPDEvent.");
  115. if (!fEvents) {
  116. std::cout << "[WARNING] MpdMiniDstFillTask::Init - No MpdEvent has been found"
  117. << std::endl;
  118. }
  119. // Get information on reconstructed vertex
  120. fVertices = (TClonesArray*) ioman->GetObject("Vertex");
  121. if (!fVertices) {
  122. std::cout << "[WARNING] MpdMiniDstFillTask::Init - No Vertex has been found"
  123. << std::endl;
  124. }
  125. // Retrieve information about tracks in TPC
  126. fTpcTracks = (TClonesArray*) ioman->GetObject("TpcKalmanTrack");
  127. if (!fTpcTracks) {
  128. std::cout << "[WARNING] MpdMiniDstFillTask::Init - No TpcKalmanTrack has been found"
  129. << std::endl;
  130. }
  131. // Retrieve information about TOF hits
  132. fTofHits = (TClonesArray*) ioman->GetObject("TOFHit");
  133. if (!fTofHits) {
  134. std::cout << "[WARNING] MpdMiniDstFillTask::Init - No TOFHit has been found"
  135. << std::endl;
  136. }
  137. // Get information about TPC tracks that match TOF
  138. fTofMatching = (TClonesArray*) ioman->GetObject("TOFMatching");
  139. if (!fTofMatching) {
  140. std::cout << "[WARNING] MpdMiniDstFillTask::Init - No TOFMatching has been found"
  141. << std::endl;
  142. }
  143. // Retrieve MC Tracks
  144. fMCTracks = (TClonesArray*) ioman->GetObject("MCTrack");
  145. if (!fMCTracks) {
  146. std::cout << "[WARNING] MpdMiniDstFillTask::Init - No MCTracks have been found"
  147. << std::endl;
  148. }
  149. // Retrieve generator (primary) tracks
  150. fGenTracks = (TClonesArray*) ioman->GetObject("GenTracks");
  151. if (!fGenTracks) {
  152. std::cout << "[WARNING] MpdMiniDstFillTask::Init - No GenTracks have been found"
  153. << std::endl;
  154. }
  155. // Retrieve information about barrel ECal digits
  156. // fEmcDigits = (TClonesArray*) ioman->GetObject("EmcDigit");
  157. // if (!fEmcDigits) {
  158. // std::cout << "[WARNING] MpdMiniDstFillTask::Init - No EmcDigits have been found"
  159. // << std::endl;
  160. // }
  161. if (fIsUseECal) {
  162. // Retrieve information about barrel ECal clusters
  163. fEmcClusters = (TClonesArray*) ioman->GetObject("EmcCluster");
  164. if (!fEmcClusters) {
  165. std::cout << "[WARNING] MpdMiniDstFillTask::Init - No EmcClusters have been found"
  166. << std::endl;
  167. }
  168. } // if ( fIsUseECal )
  169. // Retrieve information about FHCal
  170. fZdcDigits = (TClonesArray*) ioman->GetObject("ZdcDigi");
  171. if (!fZdcDigits) {
  172. std::cout << "[WARNING] MpdMiniDstFillTask::Init - No ZdcDigits have been found"
  173. << std::endl;
  174. }
  175. return kSUCCESS;
  176. }
  177. //_________________
  178. void MpdMiniDstFillTask::Exec(Option_t* option) {
  179. // Clear container at the beginning of each event
  180. if (!fMcTrk2MiniMcTrk.empty()) fMcTrk2MiniMcTrk.clear();
  181. if (!fMcTrk2EcalCluster.empty()) fMcTrk2EcalCluster.clear();
  182. Bool_t isGood = isGoodEvent();
  183. // Convert only events that are normally reconstructed
  184. if (isGood) {
  185. fillEvent();
  186. fillFHCalHits();
  187. fillBTofHits();
  188. fillMcTracks();
  189. if (fIsUseECal) {
  190. fillECalClusters();
  191. }
  192. fillTracks();
  193. if (fVerbose != 0) {
  194. fMiniDst->printTracks();
  195. }
  196. // Fill TTree
  197. fTTree->Fill();
  198. }// if (isGood)
  199. else {
  200. std::cout << "[WARNING] MpdMiniDstFillTask::Exec - Bad event reconstuction. Skipping event"
  201. << std::endl;
  202. }
  203. }
  204. //_________________
  205. Bool_t MpdMiniDstFillTask::isGoodEvent() {
  206. // Bool_t isGoodAmount = kFALSE;
  207. // Bool_t isGoodOrder = kTRUE;
  208. TClonesArray* glTracks = (TClonesArray*) fEvents->GetGlobalTracks();
  209. Int_t nGlobalTracks = glTracks->GetEntriesFast();
  210. Int_t nKalmanTracks = fTpcTracks->GetEntriesFast();
  211. if (nGlobalTracks != nKalmanTracks)
  212. return kFALSE;
  213. // isGoodAmount = kTRUE;
  214. // Loop over kalman (aka TPC) tracks
  215. for (Int_t iKTrack = 0; iKTrack < nKalmanTracks; iKTrack++) {
  216. // Retrieve kalman track
  217. MpdTpcKalmanTrack* kTrack = (MpdTpcKalmanTrack*) fTpcTracks->UncheckedAt(iKTrack);
  218. if (!kTrack) continue;
  219. // Retrieve global track
  220. MpdTrack* gTrack = (MpdTrack*) glTracks->UncheckedAt(iKTrack);
  221. if (!glTracks || kTrack->GetTrackID() != gTrack->GetID()) continue;
  222. } // for (Int_t iKTracks=0; iKTracks<nKalmanTracks; iKTracks++)
  223. return kTRUE;
  224. }
  225. //_________________
  226. void MpdMiniDstFillTask::Finish() {
  227. // Write data to the output file and close it
  228. if (fOutputFile) {
  229. fOutputFile->Write();
  230. fOutputFile->Close();
  231. }
  232. }
  233. //_________________
  234. void MpdMiniDstFillTask::streamerOff() {
  235. MpdMiniEvent::Class()->IgnoreTObjectStreamer();
  236. MpdMiniTrack::Class()->IgnoreTObjectStreamer();
  237. MpdMiniBTofHit::Class()->IgnoreTObjectStreamer();
  238. MpdMiniBECalCluster::Class()->IgnoreTObjectStreamer();
  239. MpdMiniBTofPidTraits::Class()->IgnoreTObjectStreamer();
  240. MpdMiniTrackCovMatrix::Class()->IgnoreTObjectStreamer();
  241. MpdMiniFHCalHit::Class()->IgnoreTObjectStreamer();
  242. MpdMiniMcEvent::Class()->IgnoreTObjectStreamer();
  243. MpdMiniMcTrack::Class()->IgnoreTObjectStreamer();
  244. }
  245. //_________________
  246. void MpdMiniDstFillTask::createArrays() {
  247. // Create MiniDst arrays
  248. fMiniArrays = new TClonesArray*[MpdMiniArrays::NAllMiniArrays];
  249. for (Int_t iArr = 0; iArr < MpdMiniArrays::NAllMiniArrays; iArr++) {
  250. fMiniArrays[iArr] = new TClonesArray(MpdMiniArrays::miniArrayTypes[iArr],
  251. MpdMiniArrays::miniArraySizes[iArr]);
  252. }
  253. // Set pointers to the arrays
  254. fMiniDst->set(fMiniArrays);
  255. }
  256. //_________________
  257. void MpdMiniDstFillTask::fillEvent() {
  258. // Fill event information
  259. TClonesArray* miniEventHeaders = fMiniArrays[MpdMiniArrays::Event];
  260. miniEventHeaders->Delete();
  261. // Create new event information
  262. MpdMiniEvent* miniEvent = new ((*miniEventHeaders)[miniEventHeaders->GetEntriesFast()]) MpdMiniEvent();
  263. if (!miniEvent) {
  264. std::cout << "[ERROR] MpdMiniDstFillTask::fillEvent - No MiniEvent has been created"
  265. << std::endl;
  266. // In this case we need to quit (probably completely)
  267. return;
  268. }
  269. // Retrieve information about reconstructed primary vertex
  270. MpdVertex* vertex = (MpdVertex*) fVertices->UncheckedAt(0);
  271. // Primary vertex position
  272. if (vertex) {
  273. miniEvent->setPrimaryVertexPosition(vertex->GetX(),
  274. vertex->GetY(),
  275. vertex->GetZ());
  276. } else {
  277. // In case that vertex has not been reconstructed
  278. miniEvent->setPrimaryVertexPosition(-999., -999., -999.);
  279. }
  280. // Number of global tracks in the current event
  281. if (!fEvents) {
  282. std::cout << "[ERROR] MpdMiniDstFillTask::fillEvent - No event information has been found"
  283. << std::endl;
  284. return;
  285. }
  286. // Retrieve global tracks
  287. TClonesArray* glTracks = (TClonesArray*) fEvents->GetGlobalTracks();
  288. miniEvent->setNumberOfGlobalTracks(glTracks->GetEntriesFast());
  289. // Fill McEvent info
  290. TClonesArray* miniMcEventHeaders = fMiniArrays[MpdMiniArrays::McEvent];
  291. miniMcEventHeaders->Delete();
  292. // Retrieve number of primary tracks
  293. Int_t nPart = -1;
  294. Int_t nColl = -1;
  295. Double_t rp = -1;
  296. // Impact parameter
  297. Double_t b = fEventHeaders->GetB();
  298. // Run number
  299. UInt_t runId = fEventHeaders->GetRunID();
  300. // Event number
  301. UInt_t eventId = fEventHeaders->GetEventID();
  302. // Time
  303. Double_t time = fEventHeaders->GetT();
  304. // Fill MC vertex
  305. TVector3 vtx;
  306. fEventHeaders->GetVertex(vtx);
  307. // Create and fill MiniMcEvent
  308. MpdMiniMcEvent* mcEvent = new ((*miniMcEventHeaders)[miniMcEventHeaders->GetEntriesFast()])
  309. MpdMiniMcEvent(runId, eventId, rp, b, nPart, nColl, vtx, time);
  310. mcEvent->setReactionPlaneAngle(fEventHeaders->GetRotZ());
  311. }
  312. //_________________
  313. void MpdMiniDstFillTask::fillFHCalHits() {
  314. TClonesArray* miniFHCal = fMiniArrays[MpdMiniArrays::FHCalHit];
  315. miniFHCal->Delete();
  316. // Retrieve number of hits (aka channels = detId(1,2)*modId(1-45)*ch(1-42) )
  317. Int_t nFHCalHits = fZdcDigits->GetEntries();
  318. // Total number of FHCal modules (aka towers)
  319. Int_t nFHCalModules = 45;
  320. // Energy deposition in the modules (towers)
  321. Float_t eLoss[90] = {};
  322. Int_t det = 0;
  323. Int_t mod = 0;
  324. // Loop over hits
  325. for (Int_t iHit = 0; iHit < nFHCalHits; iHit++) {
  326. // Retrieve i-th tower information
  327. MpdZdcDigi *zdcHit = (MpdZdcDigi*) fZdcDigits->At(iHit);
  328. if (!zdcHit) continue;
  329. // FHCal miniHit
  330. MpdMiniFHCalHit* hit = new ((*miniFHCal)[miniFHCal->GetEntriesFast()]) MpdMiniFHCalHit();
  331. det = zdcHit->GetDetectorID();
  332. mod = zdcHit->GetModuleID();
  333. // Set hit information
  334. hit->setId(det, mod, zdcHit->GetChannelID());
  335. hit->setEDep(zdcHit->GetELoss());
  336. //hit->setAdc( zdcHit->ADC( zdcHit->GetELoss() ) );
  337. // Measure energy deposition for each module
  338. eLoss[ (det - 1) * nFHCalModules + mod - 1 ] += zdcHit->GetELoss();
  339. } // for (Int_t iTower=0; iTower<nFHCalTower; iTower++)
  340. MpdMiniEvent* miniEvent = (MpdMiniEvent*) fMiniArrays[MpdMiniArrays::Event]->At(0);
  341. // Fill energy deposition in each tower
  342. for (Int_t iTow = 0; iTow < (2 * nFHCalModules); iTow++) {
  343. miniEvent->setFHCalEnergyDepositInModule(iTow, eLoss[iTow]);
  344. }
  345. }
  346. //_________________
  347. void MpdMiniDstFillTask::fillBTofHits() {
  348. // Instantiate MpdMiniBTofHit array
  349. TClonesArray* miniToF = fMiniArrays[MpdMiniArrays::BTofHit];
  350. miniToF->Delete();
  351. if (!fTofHits)
  352. return;
  353. // Loop over TOF hits
  354. for (Int_t iHit = 0; iHit < fTofHits->GetEntriesFast(); iHit++) {
  355. // Retrieve TOF hit
  356. MpdTofHit* tofHit = (MpdTofHit*) fTofHits->UncheckedAt(iHit);
  357. if (!tofHit) continue;
  358. MpdMiniBTofHit* miniTofHit =
  359. new ((*miniToF)[miniToF->GetEntriesFast()]) MpdMiniBTofHit();
  360. miniTofHit->setId(tofHit->GetDetectorID());
  361. miniTofHit->setHitPositionXYZ(tofHit->GetX(), tofHit->GetY(), tofHit->GetZ());
  362. miniTofHit->setTime(tofHit->GetTime());
  363. } // for (Int_t iHit = 0; iHit < fTofHits->GetEntriesFast(); iHit++)
  364. }
  365. //_________________
  366. void MpdMiniDstFillTask::fillMcTracks() {
  367. //
  368. // Prepare mapping of MC tracks to reconstructed ones
  369. // and those used in barrel ECal clusters
  370. //
  371. TClonesArray* miniMcTracks = fMiniArrays[MpdMiniArrays::McTrack];
  372. miniMcTracks->Delete();
  373. // Fill McTracks if exist
  374. Bool_t isEmcTrack = kFALSE;
  375. Bool_t isGenLevelTrack = kFALSE;
  376. if (!fMCTracks)
  377. return;
  378. // Loop over MC tracks
  379. for (Int_t iMcTrk = 0; iMcTrk < fMCTracks->GetEntriesFast(); iMcTrk++) {
  380. // Retrieve MCTrack
  381. MpdMCTrack* mcTrack = (MpdMCTrack*) fMCTracks->UncheckedAt(iMcTrk);
  382. // MC track must exist
  383. if (!mcTrack) continue;
  384. // Clean variables
  385. isEmcTrack = kFALSE;
  386. isGenLevelTrack = kFALSE;
  387. // Check if MC track is a generator level track
  388. if (mcTrack->GetMotherId() == -1) isGenLevelTrack = kTRUE;
  389. // {
  390. // std::cout << "McTrk: " << iMcTrk << " pdgCode: " << mcTrack->GetPdgCode()
  391. // << " px/py/pz: " << Form("%4.2f/%4.2f/%4.2f",
  392. // mcTrack->GetPx(),
  393. // mcTrack->GetPy(),
  394. // mcTrack->GetPz())
  395. // << " MothId: " << mcTrack->GetMotherId()
  396. // << std::endl;
  397. // }
  398. if (fIsUseECal && fEmcClusters) {
  399. // Check if MC track that was used in ECal clusters
  400. for (Int_t iCluster = 0; iCluster < fEmcClusters->GetEntriesFast(); iCluster++) {
  401. // Retrieve barrel ECal cluster
  402. MpdEmcClusterKI* cluster =
  403. (MpdEmcClusterKI*) fEmcClusters->UncheckedAt(iCluster);
  404. if (!cluster) continue;
  405. // Loop over tracks in the cluster
  406. for (Int_t iTrk = 0; iTrk < cluster->GetNumberOfTracks(); iTrk++) {
  407. Int_t id = -1;
  408. Float_t eDep = -1.;
  409. cluster->GetMCTrack(iTrk, id, eDep);
  410. if (id == iMcTrk) {
  411. isEmcTrack = kTRUE;
  412. fMcTrk2EcalCluster[iMcTrk] = iCluster;
  413. break;
  414. }
  415. } // for (Int_t iTrk=0; iTrk < cluster->GetNumberOfTracks(); iTrk++)
  416. } // for (Int_t iCluster = 0; iCluster < fEmcClusters->GetEntriesFast(); iCluster++)
  417. } // if ( fIsUseECal )
  418. // Check if generator level or ECal track
  419. if (isGenLevelTrack || isEmcTrack) {
  420. // Create new MiniMcTrack
  421. MpdMiniMcTrack* miniMcTrack =
  422. new ((*miniMcTracks)[miniMcTracks->GetEntriesFast()]) MpdMiniMcTrack();
  423. if (!miniMcTrack) {
  424. std::cout << "[WARNING] MpdMiniDstFillTask::fillTracks - No miniMcTrack has been found"
  425. << std::endl;
  426. continue;
  427. }
  428. // Set McTrack information
  429. miniMcTrack->setId(iMcTrk);
  430. miniMcTrack->setPdgId(mcTrack->GetPdgCode());
  431. miniMcTrack->setPx(mcTrack->GetPx());
  432. miniMcTrack->setPy(mcTrack->GetPy());
  433. miniMcTrack->setPz(mcTrack->GetPz());
  434. miniMcTrack->setEnergy(mcTrack->GetEnergy());
  435. // Assume that there is no branch with GenTracks ...
  436. miniMcTrack->setX(-1.);
  437. miniMcTrack->setY(-1.);
  438. miniMcTrack->setZ(-1.);
  439. miniMcTrack->setT(-1.);
  440. if (isGenLevelTrack) {
  441. miniMcTrack->setIsFromGenerator(kTRUE);
  442. }
  443. if (fGenTracks && isGenLevelTrack) {
  444. for (Int_t iGenTrack = 0; iGenTrack < fGenTracks->GetEntriesFast(); iGenTrack++) {
  445. MpdGenTrack* genTrack = (MpdGenTrack*) fGenTracks->UncheckedAt(iGenTrack);
  446. if (genTrack->GetIsUsed()) continue;
  447. Double_t absMomDiff = TMath::Abs(genTrack->GetMomentum().Mag() - mcTrack->GetP());
  448. if (absMomDiff < DBL_EPSILON) {
  449. genTrack->SetIsUsed(kTRUE);
  450. TLorentzVector spaceTime = genTrack->GetCoordinates();
  451. miniMcTrack->setX(spaceTime.X());
  452. miniMcTrack->setY(spaceTime.Y());
  453. miniMcTrack->setZ(spaceTime.Z());
  454. miniMcTrack->setT(spaceTime.T());
  455. } // if (absMomDiff < DBL_EPSILON)
  456. } // for (Int_t iGenTrack = 0; iGenTrack < fGenTracks->GetEntriesFast(); iGenTrack++)
  457. } // if (fGenTracks)
  458. // Store indices
  459. fMcTrk2MiniMcTrk.push_back(std::make_pair(iMcTrk,
  460. (UShort_t) miniMcTracks->GetEntriesFast() - 1));
  461. } // if ( isGenLevelTrack || isEmcTrack )
  462. } // for (Int_t iMcTrk = 0; iMcTrk < fMCTracks->GetEntriesFast(); iMcTrk++)
  463. }
  464. //_________________
  465. void MpdMiniDstFillTask::fillECalClusters() {
  466. //
  467. // Fill barrel ECal cluster information. Need to do it here
  468. // because of the mapping
  469. //
  470. // Instantiate MpdMiniBECalCluster array
  471. TClonesArray* miniEmcClusters = fMiniArrays[MpdMiniArrays::BECalCluster];
  472. miniEmcClusters->Delete();
  473. // Check if barrel ECal clusters exist
  474. if (fEmcClusters) {
  475. // Variables for smaller and larger dispertion axes
  476. Float_t lambda1 = 0;
  477. Float_t lambda2 = 0;
  478. // Loop over barrel ECal clusters
  479. for (Int_t iCluster = 0; iCluster < fEmcClusters->GetEntriesFast(); iCluster++) {
  480. // Retrieve barrel ECal cluster
  481. MpdEmcClusterKI* cluster =
  482. (MpdEmcClusterKI*) fEmcClusters->UncheckedAt(iCluster);
  483. // Cluster must exist
  484. if (!cluster) continue;
  485. // Clear
  486. lambda1 = 0;
  487. lambda2 = 0;
  488. // Create miniCluster
  489. MpdMiniBECalCluster* miniCluster =
  490. new ((*miniEmcClusters)[miniEmcClusters->GetEntriesFast()]) MpdMiniBECalCluster();
  491. // Fill miniCluster info
  492. miniCluster->setEnergy(cluster->GetE());
  493. miniCluster->setECore(cluster->GetEcore());
  494. miniCluster->setECore1p(cluster-> GetEcore_1p());
  495. miniCluster->setECore2p(cluster-> GetEcore_2p());
  496. miniCluster->setTime(cluster->GetTime());
  497. miniCluster->setXYZ(cluster->GetX(), cluster->GetY(), cluster->GetZ());
  498. miniCluster->setDPhi(cluster->GetDPhi());
  499. miniCluster->setDz(cluster->GetDZ());
  500. miniCluster->setTrackId(miniMcIdxFromMcIdx(cluster->GetTrackIndex()));
  501. cluster->GetLambdas(lambda1, lambda2);
  502. miniCluster->setLambdas(lambda1, lambda2);
  503. miniCluster->setChi2(cluster->GetChi2());
  504. miniCluster->setNLM(cluster->GetNLM());
  505. // Fill digit info
  506. Int_t id = -1;
  507. Float_t eDep = -1.;
  508. for (Int_t iDigi = 0; iDigi < cluster->GetMultiplicity(); iDigi++) {
  509. id = -1;
  510. eDep = -1.;
  511. cluster->GetDigitParams(iDigi, id, eDep);
  512. miniCluster->addDigit(id, eDep);
  513. } // for (Int_t iDigi = 0; iDigi < cluster->GetMultiplicity(); iDigi++)
  514. // Fill MC track info
  515. for (Int_t iTrk = 0; iTrk < cluster->GetNumberOfTracks(); iTrk++) {
  516. id = -1;
  517. eDep = -1.;
  518. cluster->GetMCTrack(iTrk, id, eDep);
  519. miniCluster->addMcTrack(miniMcIdxFromMcIdx(id), eDep);
  520. } // for (Int_t iTrk=0; iTrk < cluster->GetNumberOfTracks(); iTrk++)
  521. } // for (Int_t iCluster = 0; iCluster < fEmcClusters->GetEntriesFast(); iCluster++)
  522. } // if (fEmcClusters)
  523. }
  524. //_________________
  525. void MpdMiniDstFillTask::fillTracks() {
  526. //
  527. // Fill miniTrack, corresponding covariance matrix
  528. // and BTof matching information
  529. //
  530. // Reconstructed tracks
  531. TClonesArray* miniTracks = fMiniArrays[MpdMiniArrays::Track];
  532. miniTracks->Delete();
  533. // Convariance matrices of reconstructed tracks
  534. TClonesArray* miniTrackCovMatrices = fMiniArrays[MpdMiniArrays::TrackCovMatrix];
  535. miniTrackCovMatrices->Delete();
  536. // Create TOF-matching information
  537. TClonesArray* miniBTofTraits = fMiniArrays[MpdMiniArrays::BTofPidTraits];
  538. miniBTofTraits->Delete();
  539. // Retrieve global tracks from event
  540. TClonesArray* glTracks = (TClonesArray*) fEvents->GetGlobalTracks();
  541. // Reconstructed primary vertex in event
  542. MpdVertex* vtx = (MpdVertex*) fVertices->UncheckedAt(0);
  543. // Retrieve primary track indices in the fTpcTracks array that were used
  544. // for the primary vertex reconstruction
  545. TArrayI* ind = vtx->GetIndices();
  546. std::vector< Int_t > indices;
  547. for (Int_t iEle = 0; iEle < ind->GetSize(); iEle++) {
  548. indices.push_back(ind->At(iEle));
  549. }
  550. // Reference multiplicities
  551. Int_t grefMult = 0;
  552. Int_t refMultPos = 0;
  553. Int_t refMultNeg = 0;
  554. Int_t refMultHalfPosEast = 0;
  555. Int_t refMultHalfNegEast = 0;
  556. Int_t refMultHalfPosWest = 0;
  557. Int_t refMultHalfNegWest = 0;
  558. Int_t refMult2PosEast = 0;
  559. Int_t refMult2NegEast = 0;
  560. Int_t refMult2PosWest = 0;
  561. Int_t refMult2NegWest = 0;
  562. Int_t nTofMatched = 0;
  563. Int_t nBECalMatched = 0;
  564. std::vector< Double_t> nSigma;
  565. // Fill global and primary track information. All tracks originate
  566. // from kalman track, so loop over them and copy information
  567. for (Int_t iKalmanTrk = 0; iKalmanTrk < fTpcTracks->GetEntriesFast(); iKalmanTrk++) {
  568. // Retrieve i-th TpcKalmanTrack
  569. MpdTpcKalmanTrack* kalTrack = (MpdTpcKalmanTrack*) fTpcTracks->UncheckedAt(iKalmanTrk);
  570. // Skip non-existing tracks
  571. if (!kalTrack) continue;
  572. // Create miniTrack
  573. MpdMiniTrack* miniTrack =
  574. new ((*miniTracks)[miniTracks->GetEntriesFast()]) MpdMiniTrack();
  575. // Clean vector for each track
  576. if (!nSigma.empty()) nSigma.clear();
  577. // Set track parameters
  578. miniTrack->setId(kalTrack->GetTrackID());
  579. miniTrack->setChi2(kalTrack->GetChi2());
  580. miniTrack->setNHits(kalTrack->GetNofHits() * kalTrack->Charge());
  581. // Fill track covariance matrix if needed
  582. if (fIsUseCovMatrix) {
  583. // Create miniTrack covariance matrix
  584. MpdMiniTrackCovMatrix* miniTrackCovMatrix =
  585. new ((*miniTrackCovMatrices)[miniTrackCovMatrices->GetEntriesFast()]) MpdMiniTrackCovMatrix();
  586. fillCovMatrix(kalTrack, miniTrackCovMatrix);
  587. }
  588. // Let's find primary tracks from the whole set of reconstructed
  589. // tracks from the current event
  590. Bool_t isPrimary = kFALSE;
  591. for (auto it : indices) {
  592. if (TMath::Abs(it - iKalmanTrk) == 0) {
  593. isPrimary = kTRUE;
  594. break;
  595. }
  596. }
  597. // If track is primary then fill primary momentum
  598. if (isPrimary) {
  599. // Refit track to primary vertex
  600. refit2Vp(miniTrack, iKalmanTrk, vtx);
  601. //
  602. // TODO: Here we should check TOF-matching of the primary track
  603. //
  604. // Estimate reference multiplicities
  605. if (miniTrack->nHits() > 15) {
  606. TVector3 pMom = miniTrack->pMom();
  607. Int_t charge = (miniTrack->charge() > 0) ? 1 : -1;
  608. // Central region
  609. if (TMath::Abs(pMom.Eta()) < 0.5) {
  610. (charge > 0) ? refMultPos++ : refMultNeg++;
  611. }
  612. // Halfs of the TPC
  613. if (pMom.Eta() > 0) {
  614. (charge > 0) ? refMultHalfPosWest++ : refMultHalfNegWest++;
  615. }
  616. if (pMom.Eta() < 0) {
  617. (charge > 0) ? refMultHalfPosEast++ : refMultHalfNegEast++;
  618. }
  619. if (0.5 < pMom.Eta() && pMom.Eta() < 1.) {
  620. (charge > 0) ? refMult2PosWest++ : refMult2NegWest++;
  621. }
  622. if (-1. < pMom.Eta() && pMom.Eta() < -0.5) {
  623. (charge > 0) ? refMult2PosEast++ : refMult2NegEast++;
  624. }
  625. } // if ( miniTrack->nHits() > 15 )
  626. }// if (isPrimary)
  627. else {
  628. miniTrack->setPrimaryMomentum(TVector3(0., 0., 0.));
  629. }
  630. // Retrieve global track to get more parameters that are not available
  631. // directly from kalman track
  632. MpdTrack* glTrack = (MpdTrack*) glTracks->UncheckedAt(iKalmanTrk);
  633. // Track must exist
  634. if (!glTrack || (kalTrack->GetTrackID() != glTrack->GetID())) {
  635. std::cout << "[WARNING] MpdMiniDstFillTask::fillTracks -"
  636. << " No gtrk that corresponds to ktrk"
  637. << std::endl;
  638. continue;
  639. }
  640. // Setting global track momentum at DCA to primary vertex
  641. TVector3 globMom(glTrack->GetPx(), glTrack->GetPy(), glTrack->GetPz());
  642. miniTrack->setGlobalMomentum(globMom.X(), globMom.Y(), globMom.Z());
  643. // Set nSigma and dE/dx info
  644. if (fNSigmaDedxEstimator == 0) {
  645. miniTrack->setNSigmaElectron(glTrack->GetNSigmaElectron());
  646. miniTrack->setNSigmaPion(glTrack->GetNSigmaPion());
  647. miniTrack->setNSigmaKaon(glTrack->GetNSigmaKaon());
  648. miniTrack->setNSigmaProton(glTrack->GetNSigmaProton());
  649. } else if (fNSigmaDedxEstimator == 1) {
  650. nSigma = nSigmaDedx(globMom.Mag(), glTrack->GetdEdXTPC());
  651. miniTrack->setNSigmaElectron(nSigma.at(0));
  652. miniTrack->setNSigmaPion(nSigma.at(1));
  653. miniTrack->setNSigmaKaon(nSigma.at(2));
  654. miniTrack->setNSigmaProton(nSigma.at(3));
  655. } else {
  656. std::cout << "[WARNING] Wrong fNSigmaDedxEstimator: " << fNSigmaDedxEstimator << std::endl;
  657. }
  658. miniTrack->setDedx(glTrack->GetdEdXTPC());
  659. miniTrack->setHitMap(glTrack->GetLayerHitMap());
  660. // Getting primary vertex the first primary vertex
  661. MpdVertex* vertex = (MpdVertex*) fVertices->UncheckedAt(0);
  662. // Get primary vertex position
  663. TVector3 primVertex(vertex->GetX(), vertex->GetY(), vertex->GetZ());
  664. TVector3 firstPoint(glTrack->GetFirstPointX(),
  665. glTrack->GetFirstPointY(),
  666. glTrack->GetFirstPointZ());
  667. // Physical helix instantiation
  668. MpdMiniPhysicalHelix helix(globMom, firstPoint,
  669. fBField * kilogauss,
  670. glTrack->GetCharge());
  671. double pathLength = helix.pathLength(primVertex);
  672. TVector3 dcaPosition = helix.at(pathLength);
  673. miniTrack->setOrigin(dcaPosition.X(),
  674. dcaPosition.Y(),
  675. dcaPosition.Z());
  676. // Estimate refMult of global tracks
  677. if (TMath::Abs(miniTrack->gMom().Eta()) < 0.5 &&
  678. miniTrack->charge() != 0 && miniTrack->nHits() > 15) {
  679. grefMult++;
  680. }
  681. // Store index to miniMcTrack
  682. miniTrack->setMcTrackIndex(miniMcIdxFromMcIdx(kalTrack->GetTrackID()));
  683. // Update miniMcTrack global track ID information
  684. // in case of match
  685. if (miniTrack->mcTrackIndex() >= 0) {
  686. MpdMiniMcTrack* miniMcTrack =
  687. (MpdMiniMcTrack*) fMiniArrays[MpdMiniArrays::McTrack]->At(miniMcIdxFromMcIdx(kalTrack->GetTrackID()));
  688. miniMcTrack->addGlobalTrackId((UShort_t) miniTracks->GetEntriesFast() - 1);
  689. }
  690. //
  691. //
  692. // ECal-matching information
  693. //
  694. //
  695. auto itr2 = fMcTrk2EcalCluster.find(kalTrack->GetTrackID());
  696. if (itr2 != fMcTrk2EcalCluster.end()) {
  697. miniTrack->setBECalClusterIndex(itr2->second);
  698. } else {
  699. miniTrack->setBECalClusterIndex(-1);
  700. }
  701. // Check number of BECal-matched tracks
  702. if (miniTrack->isBECalTrack()) nBECalMatched++;
  703. //
  704. //
  705. // TOF-matching information
  706. //
  707. // TODO: must recalculate the information for primary tracks
  708. // because it does not make any sense for kalman or global ones
  709. //
  710. //
  711. // Loop over TOF-matching information
  712. if (fTofMatching)
  713. for (Int_t iMatch = 0; iMatch < fTofMatching->GetEntriesFast(); iMatch++) {
  714. // Retrieve TOF-matching information
  715. MpdTofMatchingData* dataMatch = (MpdTofMatchingData*) fTofMatching->UncheckedAt(iMatch);
  716. // TOF matching information should exist
  717. if (!dataMatch || (dataMatch->GetKFTrackIndex() != iKalmanTrk)) continue;
  718. // Create barrel TOF traits
  719. MpdMiniBTofPidTraits* bTofPidTraits =
  720. new ((*miniBTofTraits)[miniBTofTraits->GetEntriesFast()]) MpdMiniBTofPidTraits();
  721. miniTrack->setBTofPidTraitsIndex(miniBTofTraits->GetEntriesFast() - 1);
  722. bTofPidTraits->setTrackIndex(miniTracks->GetEntriesFast() - 1);
  723. bTofPidTraits->setHitIndex(dataMatch->GetTofHitIndex());
  724. bTofPidTraits->setBeta(dataMatch->GetBeta());
  725. bTofPidTraits->setMomentum(dataMatch->GetMomentum());
  726. // Increment number of TOF-matched tracks
  727. nTofMatched++;
  728. break;
  729. } // for (Int_t iMatch = 0; iMatch < fTofMatching->GetEntriesFast(); iMatch++)
  730. } // for (Int_t iKalmanTrk = 0; iKalmanTrk < fTpcTracks->GetEntriesFast(); iKalmanTrk++)
  731. //
  732. // Fill MpdMiniEvent with reference multiplicity values
  733. //
  734. MpdMiniEvent* miniEvent = (MpdMiniEvent*) fMiniArrays[MpdMiniArrays::Event]->At(0);
  735. miniEvent->setRefMultNeg(refMultNeg);
  736. miniEvent->setRefMultPos(refMultPos);
  737. miniEvent->setRefMultHalfNegEast(refMultHalfNegEast);
  738. miniEvent->setRefMultHalfPosEast(refMultHalfPosEast);
  739. miniEvent->setRefMultHalfNegWest(refMultHalfNegWest);
  740. miniEvent->setRefMultHalfPosWest(refMultHalfPosWest);
  741. miniEvent->setRefMult2NegEast(refMult2NegEast);
  742. miniEvent->setRefMult2PosEast(refMult2PosEast);
  743. miniEvent->setRefMult2NegWest(refMult2NegWest);
  744. miniEvent->setRefMult2PosWest(refMult2PosWest);
  745. miniEvent->setGRefMult(grefMult);
  746. miniEvent->setNumberOfBTOFMatch(nTofMatched);
  747. miniEvent->setNumberOfBECalMatch(nBECalMatched);
  748. }
  749. //_________________
  750. void MpdMiniDstFillTask::refit2Vp(MpdMiniTrack* miniTrack, Int_t iTpcKalmanTrack,
  751. MpdVertex* vtx) {
  752. // Get primary tracks from event
  753. // Done by smooth tracks from primary vertex (update momentum and track length -
  754. // covariance matrix is not updated !!!)
  755. // Got from MpdKfPrimaryVertexFinder::Smooth()
  756. MpdKalmanHit hit;
  757. TMatrixD c(3, 3), xk(3, 1), ck0(5, 1);
  758. TMatrixD a(5, 3), b(5, 3);
  759. TVector3 vert;
  760. vtx->Position(vert);
  761. xk(0, 0) = vert.X();
  762. xk(1, 0) = vert.Y();
  763. xk(2, 0) = vert.Z();
  764. Double_t rad = vert.Pt();
  765. MpdKalmanTrack* track =
  766. (MpdKalmanTrack*) fTpcTracks->UncheckedAt(iTpcKalmanTrack);
  767. MpdKalmanTrack* trVert = new MpdKalmanTrack(*track);
  768. MpdKalmanTrack track1 = *track;
  769. track1.SetParamNew(*track1.GetParam());
  770. track1.SetPos(track1.GetPosNew());
  771. track1.ReSetWeight();
  772. track1.SetLength(0.);
  773. TMatrixD g = *track1.GetWeight(); // track weight matrix
  774. if (track->GetNode() == "") {
  775. hit.SetType(MpdKalmanHit::kFixedR);
  776. hit.SetPos(track->GetPos());
  777. } else {
  778. hit.SetType(MpdKalmanHit::kFixedP);
  779. TString detName = track->GetNode();
  780. if (track->GetUniqueID()) {
  781. // ITS
  782. detName = detName(16, detName.Length());
  783. detName += "#0";
  784. }
  785. MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
  786. hit.SetDetectorID(geo->DetId(detName));
  787. // Find distance from the current track position to the last point (plane) -
  788. // to define direction (mainly for ITS)
  789. TVector3 pos = geo->GlobalPos(&hit);
  790. TVector3 norm = geo->Normal(&hit);
  791. Double_t v7[7] = {0.0};
  792. track1.SetNode("");
  793. MpdKalmanFilter::Instance()->SetGeantParamB(&track1, v7, 1);
  794. Double_t d = -(pos * norm); // Ax+By+Cz+D=0, A=nx, B=ny, C=nz
  795. TVector3 v3(v7[0], v7[1], v7[2]);
  796. d += v3 * norm;
  797. if (d < 0) {
  798. track1.SetDirection(MpdKalmanTrack::kOutward);
  799. }
  800. } // else
  801. MpdKalmanFilter::Instance()->PropagateToHit(&track1, &hit, kTRUE, kTRUE);
  802. //TMatrixD* par2 = track->GetParamNew();
  803. //cout << par2 << endl;
  804. computeAandB(xk, track, track1, a, b, ck0); // compute matrices of derivatives
  805. // W = (Bt*G*B)'
  806. TMatrixD tmp(g, TMatrixD::kMult, b);
  807. TMatrixD w(b, TMatrixD::kTransposeMult, tmp);
  808. w.Invert();
  809. TMatrixD m = *track1.GetParamNew();
  810. m -= ck0; // m-ck0
  811. // qk = W*Bt*G*(m-ck0-A*xk)
  812. TMatrixD tmp21(a, TMatrixD::kMult, xk);
  813. tmp21 *= -1;
  814. tmp21 += m; // m-ck0-A*xk
  815. TMatrixD tmp22(g, TMatrixD::kMult, tmp21);
  816. TMatrixD tmp23(b, TMatrixD::kTransposeMult, tmp22);
  817. TMatrixD qk(w, TMatrixD::kMult, tmp23);
  818. // Update momentum and last coordinate
  819. TMatrixD* parPointer = nullptr;
  820. if (track->GetParamNew())
  821. parPointer = track->GetParamNew();
  822. else
  823. parPointer = &m;
  824. TMatrixD par = *parPointer;
  825. for (Int_t i = 0; i < 3; ++i) par(i + 2, 0) = qk(i, 0);
  826. par(0, 0) = rad * vert.Phi();
  827. par(1, 0) = vert.Z();
  828. trVert->SetParam(par);
  829. trVert->SetPosNew(rad);
  830. Double_t pT = 1. / trVert->GetParam(4); // Signed pT
  831. Double_t phi = trVert->GetParam(2);
  832. Double_t theta = TMath::PiOver2() - trVert->GetParam(3);
  833. Double_t Px = TMath::Abs(pT) * TMath::Cos(phi);
  834. Double_t Py = TMath::Abs(pT) * TMath::Sin(phi);
  835. Double_t Pz = -1.;
  836. if (TMath::Sin(theta) != 0.) {
  837. Pz = TMath::Abs(pT) / TMath::Tan(theta);
  838. }
  839. miniTrack->setPrimaryMomentum(TVector3(Px, Py, Pz));
  840. }
  841. //_________________
  842. void MpdMiniDstFillTask::computeAandB(TMatrixD &xk0, const MpdKalmanTrack *track,
  843. const MpdKalmanTrack &trackM,
  844. TMatrixD &a, TMatrixD &b, TMatrixD &ck0) {
  845. // Compute matrices of derivatives w.r.t. vertex coordinates and track momentum
  846. Double_t vert0[3], *vert = xk0.GetMatrixArray(); //zero[3] = {0},
  847. for (Int_t i = 0; i < 3; ++i) vert0[i] = vert[i];
  848. MpdKalmanTrack trackk = *track;
  849. trackk.SetPos(trackk.GetPosNew());
  850. //trackk.GetParam()->Print();
  851. // Propagate track to PCA w.r.t. point xk0
  852. MpdKalmanFilter::Instance()->FindPca(&trackk, vert0);
  853. //MpdKalmanFilter::Instance()->FindPca(&trackk,zero); // just for test
  854. //std::cout << trackk.GetPosNew() << std::endl;
  855. trackk.SetParam(*trackk.GetParamNew());
  856. //trackk.GetParam()->Print();
  857. // Put track at xk0
  858. Double_t r = TMath::Sqrt(vert0[0] * vert0[0] + vert0[1] * vert0[1]);
  859. Double_t phi = trackk.GetParamNew(2); // track Phi
  860. if (r > 1.e-7) phi = TMath::ATan2(vert0[1], vert0[0]);
  861. trackk.SetPos(r);
  862. trackk.SetParam(0, r * phi);
  863. trackk.SetParam(1, vert0[2]);
  864. trackk.SetNode("");
  865. MpdKalmanTrack track0 = trackk;
  866. // Propagate track to chosen radius
  867. MpdKalmanHit hit;
  868. //hit.SetR(35.);
  869. //hit = *(MpdKalmanHitR*)track->GetTrHits()->Last();
  870. if (track->GetNode() == "") {
  871. hit.SetType(MpdKalmanHit::kFixedR);
  872. //hit.SetR(35.);
  873. //hit = *(MpdKalmanHitR*)track->GetTrHits()->Last();
  874. hit.SetPos(track->GetPos());
  875. MpdKalmanFilter::Instance()->PropagateParamR(&trackk, &hit, kFALSE);
  876. //trackk.GetParamNew()->Print();
  877. Proxim(trackM, trackk);
  878. } else {
  879. hit.SetType(MpdKalmanHit::kFixedP);
  880. TString detName = track->GetNode();
  881. if (track->GetUniqueID()) {
  882. // ITS
  883. detName = detName(16, detName.Length());
  884. detName += "#0";
  885. }
  886. MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
  887. hit.SetDetectorID(geo->DetId(detName));
  888. // Find distance from the current track position to the last point (plane) -
  889. // to define direction (mainly for ITS)
  890. TVector3 pos = geo->GlobalPos(&hit);
  891. TVector3 norm = geo->Normal(&hit);
  892. Double_t v7[7] = {0.0};
  893. MpdKalmanFilter::Instance()->SetGeantParamB(&trackk, v7, 1);
  894. Double_t d = -(pos * norm); // Ax+By+Cz+D=0, A=nx, B=ny, C=nz
  895. TVector3 v3(v7[0], v7[1], v7[2]);
  896. d += v3 * norm;
  897. if (d < 0) trackk.SetDirection(MpdKalmanTrack::kOutward);
  898. MpdKalmanFilter::Instance()->PropagateParamP(&trackk, &hit, kFALSE, kTRUE);
  899. track0.SetDirection(trackk.GetDirection());
  900. }
  901. //Double_t shift = 0.01; // 100 um coordinate shift
  902. Double_t shift = 0.1; // 1 mm coordinate shift
  903. for (Int_t i = 0; i < 3; ++i) {
  904. MpdKalmanTrack track1 = track0;
  905. vert0[i] += shift;
  906. if (i > 0) vert0[i - 1] -= shift;
  907. r = TMath::Sqrt(vert0[0] * vert0[0] + vert0[1] * vert0[1]);
  908. if (r > 1.e-7) phi = TMath::ATan2(vert0[1], vert0[0]);
  909. else phi = track0.GetParamNew(2); // track Phi
  910. track1.SetPos(r);
  911. track1.SetParam(0, r * phi);
  912. track1.SetParam(1, vert0[2]);
  913. if (track->GetNode() == "") {
  914. MpdKalmanFilter::Instance()->PropagateParamR(&track1, &hit, kFALSE);
  915. Proxim(trackk, track1);
  916. //proxim(track1,trackk);
  917. } else MpdKalmanFilter::Instance()->PropagateParamP(&track1, &hit, kFALSE, kTRUE);
  918. // Derivatives
  919. for (Int_t j = 0; j < 5; ++j) {
  920. a(j, i) = (track1.GetParamNew(j) - trackk.GetParamNew(j)) / shift;
  921. }
  922. } // for (Int_t i = 0; i < 3; ++i)
  923. for (Int_t i = 0; i < 3; ++i) {
  924. MpdKalmanTrack track1 = track0;
  925. Int_t j = i + 2;
  926. shift = (*track->GetCovariance())(j, j);
  927. shift = TMath::Sqrt(shift);
  928. if (j == 4) shift *= TMath::Sign(1., -track0.GetParamNew(j)); // 1/p
  929. track1.SetParam(j, track0.GetParamNew(j) + shift);
  930. //if (j == 2 && track1.GetParamNew(j)*TMath::Sign(1.,track1.GetParamNew(j)) > TMath::Pi())
  931. //track1.SetParam(j,track0.GetParamNew(j)-shift);
  932. if (track->GetNode() == "") {
  933. MpdKalmanFilter::Instance()->PropagateParamR(&track1, &hit, kFALSE);
  934. Proxim(trackk, track1);
  935. //proxim(track1,trackk);
  936. } else MpdKalmanFilter::Instance()->PropagateParamP(&track1, &hit, kFALSE, kTRUE);
  937. // Derivatives
  938. for (Int_t k = 0; k < 5; ++k) {
  939. b(k, i) = (track1.GetParamNew(k) - trackk.GetParamNew(k)) / shift;
  940. }
  941. } // for (Int_t i = 0; i < 3; ++i)
  942. TMatrixD qk0(3, 1);
  943. for (Int_t i = 0; i < 3; ++i) qk0(i, 0) = track0.GetParamNew(i + 2);
  944. //qk0.Print();
  945. ck0 = *trackk.GetParamNew();
  946. ck0 -= TMatrixD(a, TMatrixD::kMult, xk0);
  947. ck0 -= TMatrixD(b, TMatrixD::kMult, qk0);
  948. }
  949. //________________
  950. void MpdMiniDstFillTask::Proxim(const MpdKalmanTrack &track0, MpdKalmanTrack &track) {
  951. if (track0.GetType() != MpdKalmanTrack::kBarrel) {
  952. std::cout << " !!! Implemented only for kBarrel tracks !!!" << std::endl;
  953. exit(0);
  954. }
  955. //Double_t tmp = track.GetParamNew(0);
  956. Double_t phi0 = track0.GetParamNew(0) / track0.GetPosNew();
  957. Double_t phi = track.GetParamNew(0) / track.GetPosNew();
  958. phi = MpdKalmanFilter::Instance()->Proxim(phi0, phi);
  959. TMatrixD *par = track.GetParamNew();
  960. (*par)(0, 0) = phi * track.GetPosNew();
  961. phi0 = track0.GetParamNew(2);
  962. phi = track.GetParamNew(2);
  963. phi = MpdKalmanFilter::Instance()->Proxim(phi0, phi);
  964. (*par)(2, 0) = phi;
  965. track.SetParamNew(*par);
  966. }
  967. //_________________
  968. void MpdMiniDstFillTask::fillCovMatrix(MpdTpcKalmanTrack* tpcTrack, MpdMiniTrackCovMatrix* miniTrackCovMatrix) {
  969. // Covariance matrix is stored at DCA2beamline as being written in the source track !!!
  970. Double_t rPhi0 = tpcTrack->GetParam()->GetMatrixArray()[0];
  971. Double_t Z = tpcTrack->GetParam()->GetMatrixArray()[1];
  972. Double_t Phi = tpcTrack->GetParam()->GetMatrixArray()[2];
  973. Double_t Lambda = tpcTrack->GetParam()->GetMatrixArray()[3];
  974. Double_t MinusQPt = tpcTrack->GetParam()->GetMatrixArray()[4];
  975. Double_t R = tpcTrack->GetPosNew();
  976. Double_t R0 = tpcTrack->GetPos();
  977. miniTrackCovMatrix->SetStateVector(rPhi0, Z, Phi, Lambda, MinusQPt);
  978. miniTrackCovMatrix->SetR(R);
  979. miniTrackCovMatrix->SetR0(R0);
  980. vector <Double_t> sigmas;
  981. vector <Double_t> correlations;
  982. for (Int_t i = 0; i < 5; i++)
  983. for (Int_t j = 0; j < 5; j++) {
  984. if (i == j)
  985. sigmas.push_back((*tpcTrack->GetCovariance())[i][j]);
  986. else if (i > j)
  987. correlations.push_back((*tpcTrack->GetCovariance())[i][j]);
  988. }
  989. miniTrackCovMatrix->setSigmas(sigmas);
  990. miniTrackCovMatrix->setCorrelations(correlations);
  991. }
  992. //_________________
  993. Int_t MpdMiniDstFillTask::miniMcIdxFromMcIdx(Int_t mcIdx) {
  994. Int_t ret = -1; // Not found
  995. for (UInt_t i = 0; i < fMcTrk2MiniMcTrk.size(); i++) {
  996. if (fMcTrk2MiniMcTrk.at(i).first == mcIdx) {
  997. ret = fMcTrk2MiniMcTrk.at(i).second;
  998. break;
  999. }
  1000. }
  1001. return ret;
  1002. }
  1003. //________________
  1004. std::vector< Double_t > MpdMiniDstFillTask::nSigmaDedx(Double_t p, Double_t dEdx) {
  1005. std::vector<Float_t> electronMean{
  1006. 3757.34, 3770.90, 3784.94, 3799.03, 3808.90, 3819.55, 3828.35, 3839.03,
  1007. 3851.42, 3860.54, 3868.62, 3875.06, 3883.33, 3891.81, 3898.33, 3906.15,
  1008. 3912.76, 3920.30, 3924.64, 3929.60, 3935.65, 3939.86, 3942.09, 3949.23,
  1009. 3953.97, 3958.03, 3961.38, 3966.54, 3967.48, 3972.15, 3974.84, 3976.74,
  1010. 3982.19, 3983.85, 3989.49, 3990.05, 3990.49, 3996.64, 3996.57, 3998.25,
  1011. 4000.55, 4004.68, 4004.29, 4009.05, 4011.62, 4018.04, 4016.26, 4016.22,
  1012. 4019.27, 4017.20, 4020.18, 4024.86, 4023.67, 4025.07, 4028.15, 4031.01,
  1013. 4032.08, 4035.21, 4029.78, 4043.04, 4039.36, 4040.60, 4042.16, 4041.30,
  1014. 4034.07, 4046.29, 4043.97, 4052.67, 4050.12, 4055.28, 4053.21, 4047.72,
  1015. 4051.43, 4050.56, 4057.61, 4061.93, 4053.38, 4062.25, 4066.25, 4063.20,
  1016. 4064.63, 4066.04, 4066.93, 4074.67, 4080.85, 4067.16, 4059.44, 4065.39,
  1017. 4075.59, 4078.03, 4071.79, 4072.62, 4085.47, 4077.24, 4085.17, 4085.90,
  1018. 4067.84, 4085.41, 4094.88, 4084.99, 4079.28, 4090.99, 4085.75, 4109.16,
  1019. 4085.81};
  1020. std::vector<Float_t> pionMean{
  1021. 4262.84, 4036.03, 3845.60, 3719.01, 3577.86, 3465.71, 3393.91, 3310.76,
  1022. 3226.63, 3177.69, 3136.60, 3089.69, 3041.34, 3000.81, 2984.68, 2988.91,
  1023. 3006.90, 2988.47, 2963.48, 2946.04, 2939.60, 2928.94, 2917.17, 2905.05,
  1024. 2892.70, 2880.60, 2870.70, 2868.63, 2866.40, 2862.11, 2857.59, 2852.92,
  1025. 2848.26, 2843.56, 2839.01, 2835.19, 2834.63, 2836.07, 2835.99, 2835.48,
  1026. 2834.87, 2834.14, 2833.61, 2833.06, 2832.32, 2831.84, 2832.19, 2833.61,
  1027. 2834.33, 2835.11, 2835.17, 2835.92, 2836.62, 2837.92, 2839.48, 2840.95,
  1028. 2842.79, 2844.75, 2846.28, 2848.05, 2849.51, 2851.45, 2853.13, 2854.91,
  1029. 2856.73, 2858.61, 2859.91, 2861.68, 2863.49, 2865.84, 2867.34, 2869.98,
  1030. 2871.76, 2873.91, 2876.03, 2877.96, 2880.20, 2882.43, 2884.52, 2886.60,
  1031. 2888.32, 2890.91, 2892.97, 2895.19, 2897.32, 2899.75, 2901.64, 2903.70,
  1032. 2906.02, 2908.15, 2910.38, 2912.85, 2914.53, 2916.71, 2918.79, 2920.87,
  1033. 2922.91, 2925.07, 2928.05, 2929.70, 2931.79, 2933.62, 2936.53, 2938.11,
  1034. 2940.14, 2942.70, 2944.81, 2947.09, 2948.71, 2951.36, 2953.11, 2955.20,
  1035. 2957.73, 2959.67, 2962.50, 2963.49, 2965.71, 2967.79, 2969.70, 2971.68,
  1036. 2974.71, 2976.76, 2978.99, 2980.92, 2982.49, 2983.69, 2986.54, 2987.96,
  1037. 2989.94, 2992.44, 2994.25, 2996.08, 2996.84, 3000.75, 3002.46, 3002.93,
  1038. 3005.36, 3008.21, 3009.69, 3011.52, 3012.96, 3014.92, 3016.88, 3018.95,
  1039. 3021.46, 3023.63, 3025.01, 3027.19, 3027.37, 3030.46, 3032.71, 3033.61,
  1040. 3037.03, 3037.25, 3039.41, 3041.14, 3042.50, 3046.58, 3049.13, 3049.04,
  1041. 3051.06, 3052.68, 3053.95, 3054.28, 3057.88, 3057.84, 3059.71, 3060.80,
  1042. 3063.40, 3067.77, 3066.87, 3069.38, 3071.58, 3071.18, 3074.22, 3074.04,
  1043. 3077.39, 3077.13, 3079.85, 3081.78, 3082.59, 3085.07, 3086.51, 3088.73,
  1044. 3089.27, 3089.88, 3093.16, 3095.17, 3094.81, 3097.77, 3096.85, 3102.71,
  1045. 3100.34, 3103.17, 3105.11, 3108.30, 3107.99, 3110.87, 3112.43, 3112.21,
  1046. 3113.52, 3116.44, 3118.44, 3122.84, 3119.92, 3124.31, 3125.29, 3124.75,
  1047. 3125.22, 3128.23, 3128.21, 3126.40, 3128.76, 3131.73, 3131.05, 3135.59,
  1048. 3136.53, 3133.55, 3137.90, 3136.49, 3140.83, 3143.74, 3138.88, 3142.34,
  1049. 3147.02, 3140.30, 3138.59, 3139.91, 3139.16, 3143.03, 3146.28, 3137.80,
  1050. 3134.58, 3137.14, 3138.64};
  1051. std::vector<Float_t> kaonMean{
  1052. 24902.63, 22512.42, 20653.74, 18948.00, 17529.94, 16309.94, 15078.43, 14133.55,
  1053. 13229.39, 12313.20, 11611.54, 11000.31, 10330.84, 9745.07, 9322.95, 8911.44,
  1054. 8465.27, 8038.40, 7739.63, 7482.00, 7201.13, 6908.41, 6630.43, 6425.87,
  1055. 6264.29, 6089.32, 5904.01, 5713.26, 5533.83, 5404.25, 5301.87, 5191.32,
  1056. 5074.85, 4955.59, 4838.39, 4726.18, 4637.42, 4578.42, 4512.20, 4442.14,
  1057. 4369.90, 4295.83, 4220.43, 4145.17, 4078.87, 4030.65, 3993.99, 3955.50,
  1058. 3913.45, 3870.88, 3825.12, 3779.50, 3734.01, 3689.94, 3645.96, 3611.15,
  1059. 3585.50, 3565.63, 3543.24, 3518.78, 3493.96, 3468.58, 3441.34, 3415.68,
  1060. 3388.67, 3362.41, 3336.47, 3312.95, 3293.65, 3278.79, 3267.47, 3253.99,
  1061. 3240.52, 3227.13, 3214.76, 3198.35, 3184.53, 3170.68, 3156.23, 3142.34,
  1062. 3128.98, 3116.82, 3106.47, 3096.47, 3085.51, 3074.93, 3066.53, 3059.45,
  1063. 3054.90, 3053.29, 3059.15, 3063.03, 3066.18, 3068.04, 3066.32, 3062.43,
  1064. 3057.32, 3054.12, 3047.56, 3042.65, 3038.70, 3034.49, 3028.39, 3023.09,
  1065. 3016.85, 3013.23, 3006.28, 3000.60, 2996.01, 2990.71, 2986.44, 2981.70,
  1066. 2978.35, 2973.94, 2969.85, 2966.46, 2963.08, 2963.20, 2960.06, 2956.58,
  1067. 2952.46, 2949.20, 2949.07, 2944.48, 2943.52, 2941.92, 2939.51, 2935.96,
  1068. 2934.03, 2933.08, 2929.91, 2927.59, 2925.62, 2925.12, 2920.79, 2918.91,
  1069. 2917.29, 2916.21, 2914.11, 2909.00, 2910.72, 2906.60, 2907.18, 2906.63,
  1070. 2903.04, 2902.20, 2903.17, 2902.20, 2902.07, 2898.14, 2896.82, 2899.09,
  1071. 2897.26, 2898.67, 2895.51, 2896.22, 2895.20, 2894.95, 2896.60, 2892.48,
  1072. 2894.59, 2893.13, 2891.10, 2890.79, 2892.82, 2889.15, 2888.74, 2888.66,
  1073. 2885.85, 2883.53, 2886.15, 2885.93, 2885.98, 2885.17, 2883.49, 2884.08,
  1074. 2883.31, 2884.29, 2880.32, 2885.39, 2885.55, 2884.25, 2881.87, 2883.70,
  1075. 2883.74, 2884.78, 2884.68, 2887.26, 2888.48, 2881.29, 2886.50, 2882.68,
  1076. 2880.41, 2887.89, 2885.83, 2887.97, 2885.93, 2887.57, 2884.88, 2883.68,
  1077. 2889.82, 2882.39, 2884.60, 2883.05, 2880.41, 2887.52, 2885.94, 2883.64,
  1078. 2879.01, 2886.20, 2887.72, 2886.49, 2888.96, 2880.14, 2883.86, 2889.81,
  1079. 2890.89, 2888.92, 2900.42, 2890.46, 2892.30, 2893.44, 2885.62, 2895.89,
  1080. 2890.43, 2891.32, 2895.11, 2890.97, 2889.99, 2898.71, 2892.84, 2896.20,
  1081. 2889.65, 2897.49, 2887.01};
  1082. std::vector<Float_t> protonMean{
  1083. 38972.37, 37573.08, 35969.80, 34459.91, 33019.80, 31650.71, 30150.28, 28471.70,
  1084. 26478.74, 24943.29, 23634.93, 22394.18, 21419.11, 20571.20, 19735.88, 18887.98,
  1085. 18096.00, 17437.12, 16846.17, 16207.97, 15554.46, 14939.65, 14434.02, 14004.97,
  1086. 13563.70, 13095.23, 12612.15, 12169.34, 11803.16, 11504.49, 11193.70, 10861.69,
  1087. 10516.00, 10173.81, 9866.30, 9626.29, 9428.65, 9223.85, 9006.58, 8783.11,
  1088. 8547.25, 8314.05, 8099.92, 7924.99, 7789.16, 7658.45, 7516.56, 7369.65,
  1089. 7219.12, 7064.38, 6911.46, 6761.07, 6628.94, 6523.83, 6437.91, 6351.78,
  1090. 6261.61, 6166.89, 6070.43, 5970.76, 5870.26, 5769.79, 5671.57, 5581.27,
  1091. 5506.87, 5446.44, 5392.92, 5336.64, 5277.38, 5217.60, 5155.84, 5093.41,
  1092. 5030.03, 4966.69, 4905.62, 4842.56, 4783.99, 4734.08, 4691.10, 4654.32,
  1093. 4622.63, 4587.78, 4552.15, 4516.16, 4478.47, 4439.35, 4400.60, 4361.30,
  1094. 4321.57, 4282.49, 4243.61, 4204.60, 4167.69, 4133.43, 4104.35, 4078.59,
  1095. 4057.82, 4037.31, 4017.78, 3996.53, 3973.57, 3952.16, 3929.09, 3905.43,
  1096. 3882.12, 3858.92, 3836.02, 3812.96, 3790.78, 3769.17, 3749.01, 3728.11,
  1097. 3708.20, 3689.55, 3670.49, 3652.33, 3634.04, 3615.72, 3597.41, 3581.73,
  1098. 3564.05, 3548.06, 3533.77, 3519.34, 3505.96, 3494.08, 3481.90, 3471.18,
  1099. 3460.27, 3448.81, 3437.92, 3427.87, 3417.28, 3406.27, 3395.20, 3386.34,
  1100. 3375.71, 3365.49, 3355.56, 3344.80, 3334.77, 3323.59, 3314.08, 3303.63,
  1101. 3294.71, 3285.58, 3275.74, 3265.11, 3256.75, 3249.46, 3241.93, 3235.16,
  1102. 3227.95, 3219.91, 3215.56, 3208.66, 3202.83, 3197.98, 3190.52, 3185.80,
  1103. 3180.11, 3174.12, 3168.97, 3162.83, 3157.70, 3152.17, 3146.97, 3141.88,
  1104. 3137.22, 3131.89, 3126.64, 3122.98, 3117.68, 3114.53, 3109.62, 3106.55,
  1105. 3104.21, 3102.82, 3101.72, 3102.88, 3101.63, 3103.90, 3104.95, 3105.39,
  1106. 3106.98, 3107.46, 3107.23, 3107.80, 3107.03, 3105.42, 3105.41, 3105.67,
  1107. 3101.84, 3099.91, 3098.73, 3094.83, 3093.52, 3091.87, 3088.97, 3086.46,
  1108. 3082.45, 3080.88, 3079.28, 3075.34, 3072.60, 3071.95, 3068.22, 3066.21,
  1109. 3063.35, 3061.61, 3057.89, 3056.60, 3052.09, 3051.42, 3048.96, 3045.14,
  1110. 3040.91, 3040.11, 3035.33, 3035.78, 3032.39, 3030.79, 3028.90, 3028.71,
  1111. 3024.20, 3024.21, 3020.38, 3018.97, 3015.58, 3015.54, 3010.61};
  1112. std::vector<Float_t> electronSigma{
  1113. 258.53, 258.66, 258.60, 259.82, 260.60, 262.01, 262.63, 264.18,
  1114. 266.58, 265.97, 266.51, 266.67, 267.31, 268.45, 266.98, 269.08,
  1115. 269.89, 270.52, 268.64, 269.40, 271.10, 270.64, 271.42, 272.45,
  1116. 270.73, 273.04, 274.42, 275.20, 276.39, 276.95, 277.20, 278.11,
  1117. 278.62, 279.04, 279.23, 279.30, 279.41, 278.45, 280.47, 282.71,
  1118. 282.09, 282.52, 282.94, 283.31, 286.75, 285.68, 289.90, 285.61,
  1119. 286.83, 284.81, 288.44, 285.30, 286.75, 288.75, 289.52, 287.85,
  1120. 288.96, 293.38, 292.45, 293.87, 291.02, 293.36, 294.46, 290.53,
  1121. 291.99, 297.53, 297.21, 293.31, 291.58, 296.65, 301.40, 294.78,
  1122. 294.42, 299.74, 295.31, 297.58, 293.69, 300.13, 293.72, 298.41,
  1123. 289.15, 299.23, 311.42, 297.80, 305.38, 304.40, 301.76, 296.18,
  1124. 302.25, 290.06, 307.85, 302.16, 313.76, 303.03, 298.38, 307.94,
  1125. 297.13, 303.07, 296.57, 290.59, 303.75, 293.38, 300.15, 292.96,
  1126. 289.06};
  1127. std::vector<Float_t> pionSigma{
  1128. 343.49, 323.76, 291.77, 277.81, 262.79, 244.54, 237.33, 229.45,
  1129. 220.80, 213.92, 210.48, 206.63, 202.53, 197.83, 195.85, 194.98,
  1130. 194.86, 193.02, 191.03, 189.22, 188.77, 188.25, 187.49, 186.59,
  1131. 185.96, 185.08, 184.25, 184.19, 184.34, 183.99, 183.92, 183.57,
  1132. 183.26, 183.17, 182.81, 182.70, 182.59, 182.87, 183.16, 183.00,
  1133. 183.06, 183.28, 183.26, 183.44, 183.39, 183.57, 183.65, 183.85,
  1134. 184.11, 184.31, 184.15, 184.38, 184.55, 184.78, 185.25, 185.43,
  1135. 185.87, 186.04, 186.06, 186.19, 186.61, 186.95, 187.17, 187.46,
  1136. 187.63, 187.92, 187.99, 188.42, 188.61, 188.81, 189.10, 189.43,
  1137. 189.81, 189.93, 190.42, 190.40, 190.58, 191.20, 191.00, 191.89,
  1138. 191.55, 191.65, 192.39, 192.60, 192.70, 192.92, 193.21, 194.02,
  1139. 193.96, 193.80, 194.05, 194.42, 195.15, 194.74, 195.15, 195.86,
  1140. 195.69, 196.13, 196.48, 196.64, 197.03, 197.45, 197.77, 197.26,
  1141. 197.79, 197.93, 197.78, 197.97, 198.55, 198.90, 199.04, 199.16,
  1142. 199.54, 199.16, 199.66, 200.10, 200.39, 200.81, 200.45, 200.24,
  1143. 200.93, 201.21, 201.48, 201.68, 201.68, 201.61, 202.66, 202.52,
  1144. 202.05, 201.35, 201.89, 202.37, 203.61, 203.45, 203.75, 203.88,
  1145. 203.58, 203.09, 203.31, 204.57, 204.46, 204.89, 205.14, 205.81,
  1146. 205.79, 204.50, 206.12, 206.07, 205.86, 206.40, 206.25, 206.56,
  1147. 206.88, 206.72, 206.06, 206.89, 208.78, 206.17, 208.23, 206.99,
  1148. 208.42, 208.33, 207.79, 207.94, 207.77, 206.79, 208.69, 209.45,
  1149. 209.40, 211.47, 209.18, 209.05, 211.09, 208.65, 210.92, 211.11,
  1150. 210.44, 211.16, 210.46, 211.39, 210.92, 210.99, 210.44, 211.89,
  1151. 212.97, 211.83, 210.25, 211.96, 212.69, 213.38, 212.21, 215.57,
  1152. 213.26, 214.98, 215.59, 214.87, 216.26, 216.37, 215.11, 216.17,
  1153. 216.37, 213.79, 216.14, 214.42, 215.22, 214.40, 215.69, 215.58,
  1154. 217.22, 216.90, 219.18, 217.05, 217.72, 220.35, 220.81, 219.64,
  1155. 221.73, 220.45, 219.60, 220.21, 220.55, 216.99, 226.87, 220.02,
  1156. 223.34, 218.94, 221.12, 225.28, 226.55, 224.69, 226.13, 224.44,
  1157. 221.95, 226.09, 226.29};
  1158. std::vector<Float_t> kaonSigma{
  1159. 3134.25, 2584.22, 2206.67, 1875.84, 1560.24, 1461.16, 1237.64, 1089.92,
  1160. 1079.11, 966.25, 825.77, 823.58, 805.53, 699.73, 638.04, 639.51,
  1161. 631.65, 575.04, 518.67, 510.84, 506.91, 497.52, 463.83, 427.75,
  1162. 419.60, 417.30, 412.28, 402.13, 383.97, 359.91, 353.49, 350.35,
  1163. 347.75, 342.51, 334.93, 324.09, 309.56, 304.00, 304.51, 302.02,
  1164. 299.22, 296.38, 293.62, 286.25, 278.27, 271.87, 268.80, 268.85,
  1165. 266.46, 264.79, 263.43, 261.02, 259.30, 254.83, 251.58, 246.21,
  1166. 243.71, 242.35, 241.68, 240.91, 240.07, 238.66, 237.05, 236.49,
  1167. 234.23, 232.20, 229.64, 227.52, 224.01, 223.42, 222.02, 222.10,
  1168. 222.17, 220.40, 220.11, 219.90, 219.86, 217.23, 216.55, 215.35,
  1169. 214.89, 213.22, 212.78, 210.81, 209.93, 209.11, 208.96, 207.31,
  1170. 206.47, 205.72, 207.58, 207.22, 207.81, 208.01, 206.87, 206.50,
  1171. 206.42, 206.03, 205.53, 205.32, 204.97, 204.17, 203.76, 204.34,
  1172. 202.50, 203.19, 202.43, 201.42, 201.10, 201.23, 201.14, 201.00,
  1173. 199.74, 200.99, 197.84, 199.02, 198.92, 198.01, 198.98, 198.91,
  1174. 198.03, 198.26, 199.48, 198.51, 198.77, 197.56, 199.30, 197.04,
  1175. 196.83, 197.60, 197.61, 196.84, 196.54, 197.20, 196.71, 196.54,
  1176. 195.94, 196.31, 194.34, 195.07, 196.95, 195.58, 194.33, 195.75,
  1177. 194.91, 196.49, 196.18, 196.57, 195.60, 195.58, 195.65, 197.12,
  1178. 196.27, 195.39, 196.52, 194.10, 195.39, 195.90, 194.37, 197.73,
  1179. 195.24, 198.60, 194.33, 193.05, 193.43, 195.81, 193.50, 196.00,
  1180. 197.59, 194.61, 196.57, 196.15, 197.28, 192.05, 192.17, 195.60,
  1181. 195.71, 194.63, 198.89, 195.48, 193.99, 194.54, 196.01, 197.04,
  1182. 198.56, 193.06, 200.18, 193.05, 196.19, 195.26, 194.21, 196.31,
  1183. 194.59, 194.95, 194.26, 193.06, 194.97, 197.72, 196.07, 195.60,
  1184. 199.63, 195.30, 191.97, 199.69, 198.34, 196.06, 197.30, 196.43,
  1185. 196.07, 195.15, 196.99, 202.76, 203.14, 195.08, 199.69, 202.80,
  1186. 200.81, 198.33, 203.34, 204.34, 197.69, 201.61, 197.72, 205.31,
  1187. 200.97, 200.56, 201.52, 200.08, 207.65, 195.66, 205.47, 202.64,
  1188. 200.01, 200.90, 200.07};
  1189. std::vector<Float_t> protonSigma{
  1190. 5834.30, 5463.65, 5135.04, 4974.02, 4738.06, 4659.66, 4344.44, 4001.34,
  1191. 3321.36, 2866.88, 2546.12, 2227.61, 2018.41, 1900.60, 1822.92, 1706.99,
  1192. 1562.08, 1441.58, 1409.41, 1369.35, 1304.94, 1192.38, 1082.66, 1039.39,
  1193. 1030.27, 1011.92, 967.96, 889.66, 805.44, 779.74, 776.36, 774.98,
  1194. 767.76, 738.21, 681.22, 624.61, 601.27, 601.09, 604.94, 605.25,
  1195. 604.08, 590.53, 559.30, 520.43, 499.05, 495.32, 495.82, 495.89,
  1196. 497.52, 494.00, 487.37, 474.17, 451.02, 429.55, 419.49, 415.85,
  1197. 414.80, 415.23, 412.60, 412.08, 410.63, 405.80, 396.70, 384.24,
  1198. 370.44, 360.31, 354.42, 354.02, 353.99, 352.19, 350.52, 348.99,
  1199. 346.34, 344.66, 340.97, 335.49, 329.48, 321.13, 313.79, 308.31,
  1200. 305.99, 305.07, 305.49, 304.03, 303.40, 301.85, 301.59, 299.92,
  1201. 298.26, 297.41, 293.96, 291.60, 287.60, 283.08, 278.05, 274.27,
  1202. 271.53, 270.70, 268.35, 268.06, 267.85, 267.14, 266.01, 265.72,
  1203. 264.38, 263.13, 262.19, 260.99, 259.38, 257.30, 256.34, 253.92,
  1204. 253.50, 250.85, 251.13, 249.21, 248.20, 245.96, 245.79, 244.22,
  1205. 242.20, 240.92, 238.82, 237.26, 236.12, 235.67, 234.11, 232.93,
  1206. 233.06, 231.74, 231.12, 231.64, 230.17, 229.69, 229.22, 227.97,
  1207. 227.48, 227.11, 226.48, 226.15, 225.59, 224.21, 224.02, 222.63,
  1208. 222.56, 221.45, 220.12, 220.40, 218.58, 218.86, 217.48, 216.39,
  1209. 214.97, 215.78, 214.56, 214.06, 213.02, 213.22, 213.61, 212.78,
  1210. 212.05, 211.11, 211.98, 209.97, 210.49, 210.29, 210.48, 208.62,
  1211. 208.31, 207.24, 207.73, 207.45, 206.01, 205.58, 205.71, 206.14,
  1212. 205.70, 204.38, 204.30, 204.22, 204.89, 203.99, 203.72, 203.63,
  1213. 203.22, 204.82, 202.75, 203.20, 204.08, 203.77, 202.71, 203.79,
  1214. 201.84, 201.91, 203.54, 200.67, 202.17, 201.98, 201.37, 202.27,
  1215. 199.40, 199.81, 200.03, 199.79, 200.81, 199.39, 200.92, 201.15,
  1216. 199.94, 199.82, 198.49, 200.69, 198.03, 199.58, 200.80, 197.86,
  1217. 198.21, 198.48, 197.28, 198.41, 196.78, 196.24, 198.17, 196.89,
  1218. 198.92, 195.95, 196.78, 197.02, 197.05, 195.93, 195.33};
  1219. std::vector< Double_t > retVect; // 0-electron, 1-pion, 2-kaon, 3-proton
  1220. Double_t step = 0.01;
  1221. Double_t pRange[2] = {0.16, 1.2};
  1222. Int_t iBin = 0;
  1223. Double_t dedxMean = 0.;
  1224. Double_t dedxSigma = 1.;
  1225. // Electron
  1226. if (p >= 0.16 && p <= 1.2) {
  1227. iBin = TMath::FloorNint((p - pRange[0]) / step);
  1228. //std::cout << "e iBin: " << iBin << " size: " << electronMean.size() << std::endl;
  1229. dedxMean = electronMean.at(iBin);
  1230. dedxSigma = electronSigma.at(iBin);
  1231. } else if (p > 1.2) {
  1232. dedxMean = 3850.8 + p * 336.062 + p * p * (-115.97);
  1233. dedxSigma = 242.377 + p * 92.0334 + p * p * (-33.1556);
  1234. } else {
  1235. dedxMean = 3495.38 + p * 2068.7 + p * p * (-2379.39);
  1236. dedxSigma = 233.202 + p * 191.953 + p * p * (-244.544);
  1237. }
  1238. retVect.push_back((dEdx - dedxMean) / dedxSigma);
  1239. pRange[2] = 2.2;
  1240. // Pion
  1241. if (p >= 0.16 && p <= 2.2) {
  1242. iBin = TMath::FloorNint((p - pRange[0]) / step);
  1243. //std::cout << "pi iBin: " << iBin << " size: " << pionMean.size() << std::endl;
  1244. dedxMean = pionMean.at(iBin);
  1245. dedxSigma = pionSigma.at(iBin);
  1246. } else if (p > 2.2) {
  1247. dedxMean = 2631.15 + p * 302.052 + p * p * (-36.1923);
  1248. dedxSigma = 174.48 + p * 20.5397 + p * p * (-0.689454);
  1249. } else {
  1250. dedxMean = 10307.1 + p * (-56934.1) + p * p * 114991;
  1251. dedxSigma = 916.58 + p * (-5205.87) + p * p * 9524.84;
  1252. }
  1253. retVect.push_back((dEdx - dedxMean) / dedxSigma);
  1254. // Kaon
  1255. if (p >= 0.16 && p <= 2.2) {
  1256. iBin = TMath::FloorNint((p - pRange[0]) / step);
  1257. //std::cout << "K iBin: " << iBin << " size: " << kaonMean.size() << std::endl;
  1258. dedxMean = kaonMean.at(iBin);
  1259. dedxSigma = kaonSigma.at(iBin);
  1260. } else if (p > 2.2) {
  1261. dedxMean = 3455.91 + p * (-565.856) + p * p * 139.353;
  1262. dedxSigma = 273.621 + p * (-88.8921) + p * p * 25.0721;
  1263. }
  1264. else {
  1265. dedxMean = 68436.9 + p * (-382891) + p * p * 627999;
  1266. dedxSigma = 11429.8 + p * (-77924.7) + p * p * 142956;
  1267. }
  1268. retVect.push_back((dEdx - dedxMean) / dedxSigma);
  1269. // Proton
  1270. if (p >= 0.2 && p <= 2.2) {
  1271. iBin = TMath::FloorNint((p - pRange[0]) / step);
  1272. //std::cout << "p iBin: " << iBin << " size: " << protonMean.size() << std::endl;
  1273. dedxMean = protonMean.at(iBin);
  1274. dedxSigma = protonSigma.at(iBin);
  1275. } else if (p > 2.2) {
  1276. dedxMean = 4238.02 + p * (-788.916) + p * p * 119.585;
  1277. dedxSigma = 313.054 + p * (-83.9304) + p * p * 14.9712;
  1278. } else {
  1279. dedxMean = 78147.3 + p * (-258875) + p * p * 253721;
  1280. dedxSigma = 24925.3 + p * (-123173) + p * p * 161350;
  1281. }
  1282. retVect.push_back((dEdx - dedxMean) / dedxSigma);
  1283. return retVect;
  1284. }