#include "hparticletool.h" #include "htool.h" #include "hphysicsconstants.h" #include "hcategorymanager.h" #include "hgeomvector.h" #include "hmdcdedx2.h" #include "henergylosscorrpar.h" #include "hlinearcategory.h" //--------category definitions--------- #include "hparticledef.h" #include "haddef.h" #include "hmdcdef.h" #include "hmdctrackddef.h" #include "hmdctrackgdef.h" #include "richdef.h" #include "rpcdef.h" #include "showerdef.h" #include "tofdef.h" #include "walldef.h" #include "hstartdef.h" //------------------------------------- //-------objects----------------------- #include "hmdcgeomobj.h" #include "hgeantkine.h" #include "hgeanttof.h" #include "hgeantshower.h" #include "hgeantrpc.h" #include "hparticlecand.h" #include "hparticlecandsim.h" #include "hparticlepairmaker.h" #include "hmetamatch2.h" #include "hrichhit.h" #include "hrichhitsim.h" #include "htofhit.h" #include "htofhitsim.h" #include "htofcluster.h" #include "htofclustersim.h" #include "hrpccluster.h" #include "hrpcclustersim.h" #include "hrpchit.h" #include "hshowerhit.h" #include "hshowerhitsim.h" #include "hstart2hit.h" #include "hstart2cal.h" #include "hmdctrkcand.h" #include "hmdcseg.h" #include "hmdcsegsim.h" #include "hmdchit.h" #include "hmdcclusinf.h" #include "hmdcclusfit.h" #include "hmdcwirefit.h" #include "hmdccal1.h" #include "hmdcclus.h" #include "hmdclayer.h" //------------------------------------- #include "TMath.h" #include "TAxis.h" #include "TVector3.h" #include "TLatex.h" #include "TSystem.h" //_HADES_CLASS_DESCRIPTION /////////////////////////////////////////////////////////// // HParticleTool // // library of static functions // look here for auxiliary functions for: // physics analysis, simulation analysis, ... /////////////////////////////////////////////////////////// ClassImp(HParticleTool) Float_t HParticleTool::rpcCellHalfWidth[192]={ 11,11,11,11,11,11,11,11,11,11,11,11,11,13.5,13.5,13.5,13.5,13.5,16,16,16,16,16,18.5,18.5,18.5,18.5,21,21,21,21,0, 11,11,11,11,11,11,11,11,11,11,11,11,11,13.5,13.5,13.5,13.5,13.5,16,16,16,16,16,18.5,18.5,18.5,18.5,21,21,21,21,0, 11,11,11,11,11,11,11,11,11,11,11,11,11,13.5,13.5,13.5,13.5,13.5,16,16,16,16,16,18.5,18.5,18.5,18.5,21,21,21,21,0, 11,11,11,11,11,11,11,11,11,11,11,11,11,13.5,13.5,13.5,13.5,13.5,16,16,16,16,16,18.5,18.5,18.5,18.5,21,21,21,21,21, 11,11,11,11,11,11,11,11,11,11,11,11,11,13.5,13.5,13.5,13.5,13.5,16,16,16,16,16,18.5,18.5,18.5,18.5,21,21,21,21,0, 11,11,11,11,11,11,11,11,11,11,11,11,11,13.5,13.5,13.5,13.5,13.5,16,16,16,16,16,18.5,18.5,18.5,18.5,21,21,21,21,0 }; Float_t HParticleTool::tofCellHalfWidth[64]={ 15,15,15,15,15,15,15,15, 15,15,15,15,15,15,15,15, 15,15,15,15,15,15,15,15, 15,15,15,15,15,15,15,15, 10,10,10,10,10,10,10,10, 10,10,10,10,10,10,10,10, 10,10,10,10,10,10,10,10, 10,10,10,10,10,10,10,10 }; Double_t HParticleTool::parsSX[6][8][8][3] = { { { { 11.261, 21.546, -1.48157} , { 9.12987, 20.0815, -0.90077} , { 4.9937, 25.1566, -0.676999} , { -16.9104, 45.4421, -0.266398} , { 9.32513, 20.5801, -0.898808} , { -4.59022, 34.1841, -0.37284} , { -2.35385, 27.8868, -0.403455} , { -11.9264, 42.6951, -0.303151} , },{ { 4.60186, 19.8076, -0.587582} , { 9.14158, 14.3592, -0.894415} , { 9.90876, 16.0818, -0.991311} , { 8.46396, 16.5918, -0.843654} , { 6.32363, 17.3507, -0.664158} , { 0.0934424, 23.3906, -0.404736} , { 1.643, 21.44, -0.460577} , { 5.25953, 20.5083, -0.603152} , },{ { 2.85999, 24.5743, -0.508187} , { -10.151, 37.5395, -0.303762} , { 10.8058, 28.3838, -0.825003} , { 1.08081, 27.8462, -0.494962} , { -34.9677, 67.8661, -0.184961} , { 1.79947, 22.2413, -0.462287} , { -5.1119, 30.06, -0.374308} , { 3.37927, 22.1526, -0.557025} , },{ { 10.0322, 17.8041, -0.999337} , { 9.86126, 21.3776, -0.942591} , { 11.0032, 19.0596, -1.0397} , { 6.19078, 20.1642, -0.711014} , { 8.83717, 23.2468, -0.812281} , { -1.56601, 32.505, -0.44904} , { 9.53384, 23.7915, -0.879889} , { 8.59978, 28.9993, -0.687787} , },{ { 11.994, 14.8531, -1.20943} , { 10.7492, 13.6447, -1.1457} , { 10.95, 13.0478, -1.25738} , { 11.0514, 12.4287, -1.20815} , { 10.8234, 11.8793, -1.34116} , { 5.68873, 16.8206, -0.587501} , { 6.5672, 15.4772, -0.658939} , { 2.446, 23.0305, -0.485263} , },{ { 9.88012, 13.4981, -1.07661} , { 7.19276, 15.9658, -0.780948} , { 6.96557, 12.8213, -0.734747} , { 9.55223, 12.1896, -1.12011} , { 10.1804, 14.1742, -1.15649} , { 10.6219, 12.8026, -1.21674} , { 5.85022, 15.4349, -0.706955} , { 7.40654, 15.201, -0.767572} , },{ { 4.81399, 14.6362, -0.565544} , { 5.74664, 14.6309, -0.690463} , { 7.36668, 11.378, -0.821489} , { 5.95364, 12.972, -0.672093} , { 4.95611, 14.3107, -0.595363} , { 6.2104, 11.389, -0.724624} , { 6.66659, 11.0339, -0.727737} , { 7.63816, 13.211, -0.848946} , },{ { 8.1091, 13.0805, -0.774877} , { -25.7058, 47.7744, -0.159777} , { 5.62307, 13.9438, -0.661224} , { 12.7875, 10.2908, -1.666} , { 24.7024, 0.941701, -1.34615} , { 7.12599, 11.7903, -0.649831} , { 26.0491, 2.63838, -1.35853} , { 7.11186, 16.2064, -0.733903} , },},{ { { -69.7927, 95.6621, -0.102507} , { 11.5704, 14.7575, -0.956814} , { -3.81425, 30.7898, -0.590318} , { 3.11458, 28.9265, -0.729786} , { 1.47892, 24.3838, -0.515407} , { 6.23517, 21.564, -0.736383} , { -116.654, 145.583, -0.0751373} , { 9.90404, 15.9537, -0.938642} , },{ { -4.6888, 28.832, -0.357614} , { 2.77237, 21.3925, -0.473694} , { -59.4308, 86.4436, -0.128484} , { 2.34314, 21.0521, -0.488452} , { -2.55563, 25.2159, -0.412345} , { -7.02359, 28.9398, -0.286315} , { -4.03711, 28.8341, -0.332087} , { 1.05908, 21.9342, -0.453614} , },{ { 5.84976, 16.2704, -0.653181} , { 3.4281, 20.2133, -0.480097} , { 10.8675, 11.5562, -1.29586} , { 3.65647, 20.3825, -0.545005} , { -1.37892, 27.5623, -0.397713} , { -0.469378, 23.1732, -0.379255} , { -606.651, 633.04, -0.00706821} , { 3.69831, 19.8863, -0.55673} , },{ { 0.485102, 32.7648, -0.461632} , { 3.86282, 23.4211, -0.565574} , { 5.88116, 18.9384, -0.615225} , { 1.52918, 28.3805, -0.480573} , { -4.33597, 35.8219, -0.399291} , { -5.21571, 33.7346, -0.356766} , { -45.0971, 76.8997, -0.141939} , { -15.6641, 49.4951, -0.269886} , },{ { 3.00747, 19.6862, -0.508539} , { 8.16869, 13.8697, -0.946434} , { 4.23195, 17.045, -0.498318} , { 3.67526, 16.7018, -0.581303} , { 7.3296, 13.6046, -0.743481} , { 1.49252, 16.9971, -0.383766} , { 1.84125, 17.1814, -0.440954} , { 1.36042, 18.9722, -0.489417} , },{ { 3.65921, 17.4074, -0.485612} , { -36.2335, 61.5436, -0.133392} , { 6.7341, 12.3367, -0.740652} , { 5.87101, 13.6511, -0.624764} , { 8.08495, 17.5246, -0.94347} , { -3.47221, 21.5962, -0.305992} , { 4.49662, 15.0783, -0.515978} , { 5.63727, 14.7944, -0.675401} , },{ { -7.49479, 29.1283, -0.279881} , { 6.41769, 12.9197, -0.712292} , { 5.62732, 12.4936, -0.641615} , { 8.86575, 12.378, -1.0114} , { 8.14724, 14.0978, -0.957321} , { 4.60103, 13.3806, -0.563378} , { -0.0265961, 21.8219, -0.383856} , { 6.2296, 13.0373, -0.657273} , },{ { 5.38258, 12.5581, -0.562009} , { 5.87885, 11.5268, -0.621051} , { 3.30387, 14.705, -0.493612} , { 8.07418, 11.9762, -0.705153} , { 5.85592, 14.1069, -0.629922} , { 0.723188, 18.8489, -0.395128} , { 4.77045, 12.5237, -0.536745} , { 5.11269, 12.8383, -0.576576} , },},{ { { 7.06766, 23.2525, -1.3624} , { -7.29368, 30.838, -0.271138} , { -15772.8, 15794.3, -0.000257984} , { -1561.64, 1586.64, -0.00550739} , { -15.3113, 37.2987, -0.230762} , { -7.91799, 29.923, -0.284048} , { -10.5158, 35.6801, -0.292342} , { -10.177, 33.9929, -0.266036} , },{ { 4.65883, 19.136, -0.519569} , { 4.49183, 15.3622, -0.535768} , { 1.86781, 20.7151, -0.412712} , { -1.88466, 23.1346, -0.338045} , { 1.01231, 18.2131, -0.404546} , { -6.47853, 26.9981, -0.25901} , { -6.04414, 27.2149, -0.266319} , { 0.794279, 21.7028, -0.421218} , },{ { 2.91705, 16.7499, -0.434655} , { 2.08053, 19.5917, -0.470595} , { -0.313242, 20.801, -0.364572} , { -3.93203, 24.5866, -0.318383} , { 1.40376, 18.326, -0.397523} , { 4.10617, 14.9937, -0.49698} , { -4.71358, 25.8243, -0.270482} , { 2.95265, 17.8282, -0.469449} , },{ { 2.87899, 21.4919, -0.41543} , { 7.10942, 13.0378, -0.644465} , { 2.63346, 17.655, -0.422695} , { -6.40921, 25.2705, -0.243451} , { -0.299593, 20.9094, -0.357914} , { -4.76231, 25.0262, -0.251614} , { -2.91885, 22.3792, -0.275383} , { -802.56, 822.202, -0.00557972} , },{ { -4.72067, 26.1644, -0.308245} , { 4.63827, 15.6288, -0.530875} , { 5.60232, 13.44, -0.662811} , { 3.65147, 16.5079, -0.540377} , { 3.05993, 16.5204, -0.490319} , { 0.487647, 19.3321, -0.377115} , { 3.55346, 14.9829, -0.480615} , { 0.832636, 18.7244, -0.408394} , },{ { 4.41073, 19.1973, -0.54619} , { 5.21161, 12.3759, -0.579481} , { 2.98157, 17.1684, -0.468315} , { 5.10093, 14.2564, -0.646358} , { 0.531921, 19.6236, -0.387677} , { -1.36747, 19.3473, -0.36012} , { -13.6473, 32.9218, -0.201639} , { 5.62505, 14.0285, -0.632473} , },{ { 3.42427, 16.3945, -0.476595} , { 1.82528, 17.9768, -0.41833} , { 2.40902, 18.6981, -0.482663} , { 2.06014, 14.2975, -0.394341} , { 1.16583, 17.3155, -0.364563} , { -5.31619, 22.6872, -0.262787} , { -0.783183, 19.091, -0.349918} , { 0.493681, 19.4715, -0.394497} , },{ { 5.62189, 12.7088, -0.533233} , { 7.92305, 8.88519, -0.82643} , { 3.45005, 14.1677, -0.450129} , { 5.75466, 10.4253, -0.528091} , { 3.95846, 12.1817, -0.464105} , { 1.5955, 15.6129, -0.392609} , { 0.166719, 17.4955, -0.351714} , { -1.92506, 21.6318, -0.265856} , },},{ { { 1.43956, 26.7983, -0.492115} , { -2339.9, 2368.57, -0.00366577} , { 12.1958, 15.2955, -1.05854} , { -8046.9, 8071.73, -0.000838187} , { -929.556, 956.805, -0.00907381} , { -38.1838, 66.5739, -0.144491} , { 3.44004, 26.5457, -0.530688} , { 4.71503, 30.6348, -0.615336} , },{ { -3.8713, 29.9232, -0.302782} , { -7.35686, 34.1268, -0.27733} , { 0.446663, 25.004, -0.405858} , { -9.75136, 36.1941, -0.27054} , { -28.8047, 55.0941, -0.156446} , { -38.287, 61.6422, -0.119819} , { -26.1488, 52.2433, -0.169803} , { -5.7388, 32.843, -0.325859} , },{ { 9.19104, 23.7934, -0.737198} , { -3.12504, 32.7559, -0.341762} , { 2.5841, 25.1687, -0.464504} , { -7.36488, 31.8188, -0.277568} , { -0.410703, 30.0252, -0.390252} , { -211.625, 241.44, -0.0417897} , { -5.65592, 28.4593, -0.284852} , { -4.17497, 36.532, -0.352869} , },{ { 5.50154, 33.0517, -0.585692} , { -2589.69, 2626.15, -0.00384649} , { -81.1282, 113.464, -0.095054} , { -0.112543, 24.0797, -0.404193} , { -48.1776, 79.2123, -0.124522} , { -75.1676, 112.116, -0.107471} , { 0.663637, 28.9076, -0.431043} , { 13.5495, 16.1299, -1.18324} , },{ { 9.55522, 14.2548, -1.03044} , { 3.12998, 20.8461, -0.47705} , { 2.4138, 19.8442, -0.432118} , { 5.76255, 15.6726, -0.573903} , { -1.80161, 27.3056, -0.452487} , { 2.83756, 19.8876, -0.44414} , { -38.5196, 64.3885, -0.136232} , { 0.453318, 23.9972, -0.445565} , },{ { 3.56884, 16.7016, -0.496302} , { 1.19605, 19.3966, -0.394533} , { 2.42802, 17.5546, -0.461305} , { 0.933455, 19.8235, -0.422565} , { -5.68463, 26.1569, -0.27655} , { -7.81786, 27.0238, -0.228276} , { 0.594207, 17.9316, -0.379158} , { 0.576988, 23.8147, -0.428638} , },{ { 3.3247, 19.0785, -0.505158} , { 5.35341, 14.9904, -0.555266} , { 0.23718, 20.9901, -0.382014} , { -2.78917, 24.1714, -0.318399} , { -23.812, 45.2813, -0.158132} , { -15.5804, 35.4139, -0.181462} , { -6.73986, 29.8659, -0.283512} , { 3.76645, 18.1942, -0.510962} , },{ { -4.17168, 26.4299, -0.326553} , { 1.9555, 18.0714, -0.447893} , { -28.6188, 53.3752, -0.172926} , { 4.644, 13.2118, -0.481685} , { 3.10493, 15.9649, -0.400782} , { -10.9767, 29.9608, -0.200803} , { 6.14032, 12.3813, -0.678371} , { 1.56021, 24.1464, -0.440946} , },},{ { { -53.6498, 82.6205, -0.151528} , { -7280.12, 7308.77, -0.00109971} , { -79.589, 103.137, -0.0657872} , { -3615.35, 3639.8, -0.00171557} , { 2.04346, 24.3134, -0.513194} , { -29.0943, 52.6841, -0.147897} , { -48.391, 71.2173, -0.104596} , { 0.975112, 22.1828, -0.464883} , },{ { 5.11937, 18.5229, -0.587031} , { 0.873546, 22.131, -0.415451} , { 5.05525, 19.354, -0.568744} , { 0.576507, 21.7937, -0.432836} , { -0.965728, 26.9736, -0.403267} , { -1.87335, 26.9608, -0.366372} , { -30.1053, 53.8164, -0.146013} , { 1.51544, 21.5872, -0.426627} , },{ { 5.84296, 19.0005, -0.630305} , { 1.88012, 22.4686, -0.460473} , { 1.44526, 22.5326, -0.443196} , { 0.846063, 22.4518, -0.401213} , { -0.225704, 22.425, -0.376266} , { -12.8648, 35.6489, -0.212786} , { -13.6876, 37.7449, -0.22641} , { -1.07116, 26.6534, -0.407683} , },{ { 11.7363, 18.71, -1.03604} , { 12.6628, 17.9114, -1.0479} , { -26.412, 56.2696, -0.19631} , { -14.0258, 42.7386, -0.259523} , { -0.149337, 22.5381, -0.382891} , { 2.88327, 22.6719, -0.481628} , { -542.257, 568.36, -0.0125465} , { 2.34131, 25.2539, -0.469079} , },{ { 5.55453, 15.637, -0.643536} , { 9.83068, 13.4938, -0.978839} , { 5.46824, 14.2818, -0.584398} , { 4.74097, 15.3457, -0.530127} , { 2.16351, 17.9984, -0.421261} , { 2.33149, 20.3425, -0.427444} , { -14.1156, 33.5424, -0.18668} , { 5.09611, 14.8257, -0.578273} , },{ { 8.50185, 13.9784, -0.87673} , { 4.5894, 15.2899, -0.618536} , { 6.53476, 12.8178, -0.757674} , { 5.71812, 14.9796, -0.65994} , { 3.89662, 14.7227, -0.519097} , { -7.14286, 27.7149, -0.270248} , { -1.11071, 19.0059, -0.341487} , { 6.88412, 13.9039, -0.702333} , },{ { 4.67483, 19.968, -0.558962} , { 6.1507, 12.6341, -0.671295} , { 5.50661, 13.9816, -0.606162} , { 0.0819361, 19.2904, -0.363735} , { -15.24, 35.2309, -0.181496} , { -3.12056, 21.0821, -0.28179} , { 7.15185, 14.3027, -0.741551} , { 5.13274, 17.5356, -0.586885} , },{ { 7.19215, 19.2935, -0.653361} , { 6.0808, 12.9599, -0.603578} , { 4.79917, 14.1376, -0.538371} , { 0.423177, 19.1295, -0.336841} , { 0.067527, 18.973, -0.380874} , { 3.62671, 13.5247, -0.470833} , { 0.159459, 18.2868, -0.415175} , { 3.23734, 22.4018, -0.516815} , },},{ { { -2.64086, 30.3344, -0.388795} , { -70.2968, 96.8009, -0.0876698} , { -3.2801, 31.1341, -0.394103} , { 9.59078, 18.1993, -1.11694} , { -126.997, 156.889, -0.0700956} , { 8.36986, 17.0035, -0.83462} , { -86.5983, 111.436, -0.0733982} , { -0.0524577, 22.6739, -0.391376} , },{ { 3.40223, 21.093, -0.474667} , { -2.03144, 24.62, -0.341968} , { 0.862482, 22.8745, -0.418965} , { 0.791488, 21.5494, -0.4124} , { 0.280344, 23.738, -0.427894} , { 3.66928, 17.663, -0.475018} , { -117.203, 142.196, -0.0536026} , { 3.66428, 22.5867, -0.519825} , },{ { 4.24259, 18.5831, -0.579052} , { -37.8992, 66.7055, -0.16185} , { 5.50227, 18.3669, -0.625435} , { 1.13764, 24.0236, -0.450364} , { 2.28567, 23.682, -0.468582} , { -4.34919, 29.2843, -0.326804} , { 0.635503, 25.9999, -0.425968} , { 4.16994, 24.5464, -0.519652} , },{ { -7.18749, 44.0448, -0.326292} , { -2087.5, 2116.64, -0.00373747} , { 5.28862, 19.9115, -0.535791} , { 9.34527, 19.0638, -0.807227} , { -37.1736, 79.3995, -0.176507} , { -3004.48, 3039.29, -0.00308339} , { 1.80121, 22.637, -0.434125} , { -5.15794, 36.9021, -0.332757} , },{ { 6.67746, 14.1903, -0.624485} , { 1.06665, 25.411, -0.44613} , { 10.0152, 12.8435, -0.923435} , { 4.70145, 19.2999, -0.478278} , { 6.73914, 14.0263, -0.596819} , { -3.2595, 27.2935, -0.306938} , { 6.55423, 14.2493, -0.621017} , { -4.50553, 30.4827, -0.326369} , },{ { -0.09777, 22.8166, -0.382884} , { 4.04781, 17.7808, -0.518584} , { 3.39938, 17.0308, -0.482232} , { 4.70489, 16.0417, -0.509594} , { 4.51307, 15.1738, -0.527934} , { 4.42407, 19.2094, -0.543449} , { 4.13092, 15.6092, -0.483003} , { 8.44791, 13.39, -0.812214} , },{ { 4.34057, 15.0788, -0.484728} , { 7.99248, 12.2559, -0.769648} , { 5.69102, 14.4769, -0.54235} , { 5.68315, 12.1356, -0.589213} , { 5.64111, 14.2835, -0.56883} , { 3.99089, 14.8788, -0.473002} , { 5.82762, 14.3199, -0.576086} , { 6.38442, 12.7365, -0.607505} , },{ { 7.1873, 10.5268, -0.560181} , { 6.40901, 13.2542, -0.632216} , { 4.62558, 13.2388, -0.482767} , { 7.41846, 12.0526, -0.649268} , { 7.73648, 8.9868, -0.725813} , { 5.81736, 11.039, -0.538824} , { 5.29344, 14.8121, -0.576218} , { 0, 21.3775, -1} } } }; Double_t HParticleTool::parsDX[6][8][8][5] = { { { { 10920.9, 7171.75, -0.0313439, -18088.8, -0.0120764} , { 2.71171, 235.054, -1.77338, -220.224, -1.55159} , { -9.94265, 32.8533, -0.757599, -4624.23, -10.4416} , { 9.68661, 215.009, -0.740932, -219.353, -0.648516} , { -4.86145, 3837.65, -2.4015, -4002.42, -2.49823} , { 1.79808, 227.414, -1.94192, -235.719, -1.69912} , { -1.71727, 219.095, -1.24102, -232.967, -1.415} , { 8.58188, 695.584, -0.823929, -687.008, -0.768741} , },{ { -0.986019, 216.25, -1.41077, -235.953, -1.62675} , { 1.32734, 245.964, -1.90624, -211.215, -1.59206} , { -0.311783, 215.94, -1.76782, -237.264, -2.10211} , { 2.77371, 69.5248, -4.29675, -14.9962, -0.505348} , { -0.57256, 945.669, -1.89666, -998.262, -2.00341} , { 1.66184, 228.319, -1.11781, -219.528, -1.00249} , { -0.685493, 936.95, -1.34718, -960.371, -1.39163} , { 1.28074, 229.156, -1.52934, -222.772, -1.36138} , },{ { -1.02279, 219.049, -1.36623, -232.778, -1.55078} , { 3.08136, 225.49, -1.15489, -222.76, -1.02584} , { -1.0673, 215.692, -1.29038, -237.907, -1.52339} , { 2.00514, 225.484, -1.41274, -226.264, -1.26091} , { -2.18165, 220.251, -0.942219, -231.956, -1.07898} , { 1.67957, 230.437, -1.60522, -222.362, -1.41469} , { -1.28441, 224.655, -1.36395, -226.346, -1.5141} , { 1.68699, 224.33, -1.52287, -228.399, -1.39004} , },{ { -1.46029, 224.81, -1.5051, -225.629, -1.68358} , { 3.14972, 74.6959, -3.85783, -23.1424, -0.628695} , { -1.37198, 217.295, -1.43446, -235.908, -1.69286} , { 3.06381, 228.795, -1.18894, -219.068, -1.0392} , { -2.16669, 218.338, -1.03162, -234.838, -1.18886} , { 2.74879, 227.513, -1.25878, -221.606, -1.09226} , { -1.47411, 810.673, -1.35115, -858.131, -1.43631} , { 3.23649, 108.004, -4.32615, -16.6157, -0.333292} , },{ { -2.18632, 31.8301, -0.830277, -77.8629, -3.16749} , { 2.80897, 61.1916, -3.14552, -22.5214, -0.652499} , { -3.5874, 20.4244, -0.579296, -53.6017, -3.08114} , { 2.2418, 117.976, -6.46199, -22.3293, -0.812699} , { -2.92012, 228.168, -0.965634, -221.532, -1.03461} , { 3.5483, 220.469, -0.944079, -227.371, -0.878963} , { -2.30826, 225.146, -0.981039, -226.005, -1.06817} , { 0.55031, 238.486, -1.80503, -214.753, -1.5673} , },{ { -2.88945, 228.839, -0.633475, -221.168, -0.658604} , { 0.284656, 217.545, -0.713614, -229.584, -0.713618} , { -1.06009, 222.895, -1.48747, -228.253, -1.6655} , { 2.10625, 225.522, -1.31835, -224.402, -1.18313} , { -2.53569, 223.334, -1.23904, -227.876, -1.4077} , { 3.59289, 221.396, -1.19902, -227.218, -1.0755} , { -1.5093, 219.213, -1.46938, -232.768, -1.70781} , { 1.31633, 232.4, -1.59572, -220.125, -1.39481} , },{ { -4.75799, 229.485, -0.73053, -220.909, -0.775424} , { 0.851585, 233.927, -1.95458, -221.449, -1.72433} , { -1.82094, 216.286, -1.46719, -235.762, -1.68484} , { 1.69405, 233.323, -1.50538, -217.69, -1.3235} , { -3.59226, 226.806, -0.903884, -224.022, -0.980582} , { 0.533236, 232.435, -1.85137, -221.457, -1.66725} , { -3.84553, 226.842, -0.866162, -224.053, -0.935188} , { -0.0263842, 242.365, -2.8499, -220.695, -2.43401} , },{ { -3.65854, 231.156, -0.473427, -218.71, -0.473425} , { -1.9784, 200.818, -0.678202, -233.381, -0.75485} , { -1.85405, 223.92, -1.45522, -226.941, -1.59028} , { 1.16225, 2.86458e+06, -31.6864, -71.8657, -1.97501} , { -1.54566, 222.344, -0.574528, -229.534, -0.620804} , { -0.932238, 228.97, -1.43971, -221.729, -1.33378} , { 1.64126, 215.226, -0.860271, -239.002, -1.05921} , { -6.26831, 225.314, -0.999892, -225.144, -1.02107} , },},{ { { -10.4473, 227.013, -0.532973, -224.03, -0.620964} , { 4.67294, 614.782, -2.60042, -586.688, -2.42656} , { -12.4126, 234.484, -0.415001, -213.329, -0.423724} , { 9.15935, 208.143, -0.889162, -217.364, -0.79372} , { -2.13442, 225.391, -1.52115, -224.547, -1.66195} , { 1.35066, 210.301, -0.979296, -238.081, -0.979296} , { -2.94905, 226.009, -0.954106, -224.96, -1.07309} , { 3.19836, 214.178, -0.405672, -226.744, -0.40569} , },{ { -2.01271, 222.754, -1.21713, -228.789, -1.37872} , { 0.599154, 232.054, -2.26889, -228.043, -1.94019} , { -1.05033, 233.296, -1.14872, -216.318, -1.22295} , { 0.333085, 234.126, -2.10192, -223.668, -1.81978} , { -2.28501, 225.336, -1.08421, -225.777, -1.19676} , { 1.13348, 228.366, -1.42978, -222.506, -1.29695} , { -1.77459, 219.108, -1.23285, -232.884, -1.38633} , { 1.2362, 223.722, -1.36396, -227.373, -1.26952} , },{ { -1.61868, 222.681, -1.53835, -228.37, -1.71795} , { 0.611119, 233.602, -1.5256, -217.253, -1.36006} , { -0.949581, 239.526, -1.48911, -206.231, -1.4891} , { 1.34715, 231.475, -1.4224, -218.944, -1.25371} , { -2.23595, 222.866, -1.24788, -228.583, -1.40806} , { 1.03861, 229.623, -1.49628, -221.66, -1.35017} , { -52.8245, 354.19, -0.0597355, -298.221, -0.0701339} , { 0.724048, 216.49, -0.998519, -234.382, -0.998519} , },{ { -0.43587, 232.25, -1.41266, -217.382, -1.51609} , { 1.2369, 229.373, -1.80808, -226.521, -1.58574} , { -1.64354, 220.6, -1.46512, -230.951, -1.66619} , { 1.95451, 227.407, -1.23184, -221.474, -1.09191} , { -6.19229, 223.152, -0.656164, -228.888, -0.761577} , { 1.01131, 225.896, -1.45941, -226.009, -1.31353} , { -2.28266, 227.674, -1.01339, -222.97, -1.11923} , { 1.83593, 230.506, -1.00472, -215.04, -0.865223} , },{ { -5.1215, 232.464, -0.638398, -217.432, -0.673103} , { 2.51892, 220.587, -1.17021, -227.986, -1.08289} , { -2.21838, 226.778, -1.18125, -223.877, -1.27081} , { 2.22758, 233.719, -1.41466, -216.074, -1.24209} , { -1.62378, 223.635, -1.23193, -227.54, -1.35362} , { 1.50002, 226.48, -1.19521, -222.886, -1.11066} , { -1.47489, 222.337, -1.1975, -229.085, -1.31802} , { 1.25824, 228.564, -1.88099, -227.411, -1.66533} , },{ { -2.09495, 226.021, -0.939606, -224.858, -0.999927} , { 1.45802, 207.794, -0.737266, -235.595, -0.737246} , { -1.85307, 220.844, -1.22236, -230.86, -1.37705} , { 2.31434, 230.209, -1.30215, -218.935, -1.13823} , { -5.18055, 223.277, -0.782532, -228.253, -0.899783} , { 1.44913, 227.946, -1.4466, -223.509, -1.30223} , { -2.74655, 222.665, -1.02146, -228.868, -1.13569} , { 2.38244, 216.923, -0.530133, -230.352, -0.530128} , },{ { -0.720597, 225.748, -1.73087, -209.107, -1.86377} , { 2.69673, 31.9755, -5.7496, -16.0353, -0.643836} , { -1.85448, 224.647, -1.11725, -226.394, -1.2018} , { 2.66561, 225.911, -1.32037, -223.997, -1.18346} , { -2.47568, 222.438, -1.16156, -228.949, -1.31596} , { 2.45744, 223.025, -1.06233, -225.945, -0.996281} , { -2.10636, 223.816, -1.22747, -227.462, -1.38016} , { 1.25296, 228.646, -1.58755, -223.518, -1.44642} , },{ { -2.52517, 225.324, -0.766246, -225.641, -0.811065} , { 1.88428, 224.48, -1.21219, -225.335, -1.12897} , { -9.11946, 232.979, -0.196329, -217.519, -0.196325} , { 1.12817, 226.246, -1.68022, -226.849, -1.54266} , { -1.88053, 224.684, -1.32809, -226.24, -1.46415} , { -1.29301, 226.564, -1.62115, -226.147, -1.49219} , { -0.700307, 221.558, -1.04414, -230.238, -1.18143} , { -2.01262, 232.922, -1.89831, -221.227, -1.70886} , },},{ { { 36.0114, 215.343, -0.414381, -224.78, -0.25242} , { 1.40997, 221.245, -2.14814, -218.551, -2.01684} , { -3.18016, 222.195, -1.07673, -229.398, -1.19871} , { 3.16919, 226.919, -1.89574, -227.336, -1.71048} , { -2.0229, 224.427, -1.17125, -226.655, -1.2839} , { 1.13062, 232.596, -2.16113, -223.326, -1.91358} , { -4.68687, 226.861, -0.741107, -224.246, -0.82604} , { 0.370807, 1406.52, -1.93987, -1402.99, -1.89538} , },{ { -8.06559, 229.297, -0.404653, -221.627, -0.442481} , { 0.467499, 229.862, -1.59129, -221.659, -1.46439} , { -0.887441, 214.619, -1.53342, -237.737, -1.79901} , { 1.18556, 230.697, -1.52264, -220.778, -1.36317} , { -1.23127, 226.99, -1.07728, -223.585, -1.13402} , { 0.570031, 227.716, -1.57323, -224.296, -1.44624} , { -1.49978, 221.51, -1.20526, -229.963, -1.30823} , { 1.39225, 224.127, -1.16791, -225.443, -1.09216} , },{ { -1.75068, 227.97, -0.802415, -222.3, -0.822587} , { 1.86597, 222.78, -1.25983, -227.698, -1.17544} , { -1.22247, 225.547, -1.18896, -225.232, -1.2697} , { 0.905307, 231.798, -1.84648, -222.723, -1.64062} , { -1.61408, 223.087, -1.23557, -228.167, -1.35528} , { 1.81932, 230.876, -1.28493, -218.107, -1.1465} , { -1.2994, 221.118, -1.12898, -230.591, -1.24158} , { 0.834589, 229.104, -1.79668, -224.277, -1.64008} , },{ { -1.41685, 218.538, -1.17526, -233.494, -1.29993} , { 4.1474, 228.146, -1.03187, -218.015, -0.901792} , { -0.571046, 222.008, -1.58215, -228.876, -1.6868} , { 0.440803, 231.64, -1.9311, -223.047, -1.7469} , { -1.47342, 221.95, -1.2046, -229.492, -1.31578} , { 0.776127, 228.154, -1.28967, -221.585, -1.1923} , { -1.48807, 221.107, -0.930545, -230.687, -1.01058} , { 4.91245, 116.914, -0.5439, -113.288, -0.443251} , },{ { -2.33914, 225.006, -0.950831, -226.332, -1.07074} , { 2.04328, 229.168, -1.43416, -221.963, -1.26802} , { -1.10915, 218.579, -1.37418, -233.345, -1.56215} , { 4.41073, 215.668, -0.491442, -230.354, -0.485419} , { -1.97165, 225.543, -1.28877, -225.229, -1.43124} , { 4.58842, 221.214, -0.522259, -225.686, -0.491307} , { -1.93997, 226.822, -1.13034, -223.766, -1.19798} , { 1.06229, 228.909, -1.79156, -225.277, -1.60886} , },{ { -0.517388, 224.521, -1.06864, -226.548, -1.14829} , { 3.80313, 221.204, -0.747589, -225.525, -0.697573} , { -1.87407, 225.686, -1.14818, -225.227, -1.25697} , { 4.27869, 222.188, -1.00387, -224.917, -0.90519} , { -0.996095, 228.195, -1.49673, -222.063, -1.63367} , { 1.76876, 225.422, -1.38877, -225.784, -1.26448} , { -9.98939, 234.889, -0.201789, -214.941, -0.201785} , { 1.13816, 224.789, -1.47993, -227.066, -1.37058} , },{ { -3.51978, 227.015, -0.735545, -223.883, -0.78877} , { 4.72942, 216.201, -0.774231, -230.062, -0.73111} , { -1.52378, 227.659, -1.4186, -222.814, -1.57858} , { 2.35644, 222.933, -0.896181, -225.289, -0.848523} , { -2.05752, 220.605, -1.13326, -231.214, -1.26643} , { 2.10405, 227.23, -1.1852, -221.644, -1.08166} , { -2.32703, 221.448, -1.1386, -230.226, -1.27063} , { 2.7436, 221.386, -0.906814, -226.765, -0.856607} , },{ { -0.718245, 226.186, -1.62612, -224.102, -1.73488} , { 1.12092, 228.83, -1.60574, -223.34, -1.4713} , { -1.63837, 224.758, -1.51789, -226.059, -1.68171} , { 1.83744, 223.307, -0.847296, -224.865, -0.805671} , { -1.61927, 223.579, -1.31641, -227.435, -1.42734} , { -1.067, 229.256, -1.27755, -220.08, -1.16674} , { 0.537486, 222.114, -1.22994, -229.436, -1.38141} , { -4.53293, 225.898, -1.21759, -224.177, -1.18394} , },},{ { { -3.00772, 214.28, -1.12287, -239.244, -1.29978} , { 5.03654, 220.856, -0.904892, -227.372, -0.852899} , { -3.22487, 213.465, -1.75347, -239.431, -2.04766} , { 0.357594, 229.604, -1.74044, -222.539, -1.61584} , { -2.57441, 226.16, -0.828441, -224.822, -0.891781} , { 0.803157, 227.245, -1.37067, -223.422, -1.23329} , { -2.21335, 216.415, -0.97004, -237.74, -1.12985} , { 0.634854, 245.146, -1.64219, -209.002, -1.32961} , },{ { -0.148042, 218.098, -1.52822, -233.698, -1.74907} , { 1.01958, 228.447, -1.32885, -221.454, -1.20501} , { -1.6397, 225.586, -0.934409, -225.475, -1.00694} , { 2.16687, 227.346, -1.0831, -220.075, -0.966315} , { -2.10292, 222.329, -0.880475, -229.528, -0.981649} , { 1.84788, 229.12, -1.19207, -218.907, -1.05185} , { -1.71528, 221.896, -1.04156, -229.876, -1.15934} , { 2.45984, 221.969, -0.893445, -225.821, -0.8366} , },{ { -0.877447, 216.227, -1.47452, -236.799, -1.72739} , { 1.98016, 233.233, -1.35281, -215.872, -1.16591} , { -6.41604, 20.7859, -0.336873, -30.4858, -3.91354} , { 1.54771, 232.957, -1.50604, -217.851, -1.3395} , { -1.51259, 222.065, -1.14936, -229.574, -1.28781} , { 3.40055, 224.18, -1.03741, -222.456, -0.913622} , { -1.63225, 225.215, -0.999986, -225.851, -1.07509} , { 2.48454, 43.2347, -3.4756, -18.3881, -0.586504} , },{ { -1.2843, 224.327, -1.07725, -226.805, -1.20716} , { 1.47831, 221.97, -1.09833, -226.866, -1.00742} , { -0.283995, 230.324, -1.35231, -219.601, -1.44109} , { 1.33297, 227.696, -1.87024, -228.073, -1.66535} , { -1.79685, 226.028, -0.912218, -224.972, -1.00956} , { 1.01842, 220.494, -0.881847, -227.339, -0.831313} , { -0.944663, 222.72, -1.1096, -228.778, -1.22537} , { 1.91399, 211.912, -0.983501, -238.675, -0.983499} , },{ { -0.870527, 234.163, -1.03176, -213.368, -1.03179} , { 3.03388, 219.245, -0.778025, -227.8, -0.733838} , { -4.00202, 227.276, -0.903705, -223.503, -0.996956} , { 2.46834, 222.22, -1.15519, -227.336, -1.06796} , { -4.50313, 33.9416, -0.723327, -42.9443, -7.3419} , { 2.63743, 227.237, -1.18535, -221.098, -1.05015} , { -1.54816, 222.745, -1.19183, -228.807, -1.34699} , { 2.01112, 220.532, -1.29186, -230.918, -1.19223} , },{ { -1.1443, 233.343, -0.86568, -215.883, -0.865688} , { 2.49983, 219.954, -0.98023, -228.964, -0.930655} , { -1.50991, 224.53, -1.22854, -226.543, -1.35897} , { 3.8495, 220.498, -0.74719, -225.807, -0.695145} , { -1.56131, 224.288, -1.28148, -226.852, -1.4324} , { 3.95265, 224.609, -0.963532, -222.152, -0.867026} , { -1.08217, 222.321, -1.24567, -228.983, -1.35381} , { 2.88146, 223.974, -1.04178, -223.887, -0.949107} , },{ { -3.65294, 228.075, -0.760829, -222.675, -0.814312} , { 1.76451, 231.991, -1.74816, -222.527, -1.53091} , { -1.61955, 223.499, -1.2995, -227.657, -1.44353} , { 2.48414, 226.186, -1.3322, -224.275, -1.19757} , { -1.19691, 226.765, -0.93687, -224.006, -0.997821} , { 4.3695, 221.036, -0.870861, -225.079, -0.792701} , { -4.66407, 227.312, -0.92848, -223.216, -1.04667} , { 5.47802, 214.74, -0.345547, -230.004, -0.345546} , },{ { -3.52474, 226.816, -0.96626, -224.003, -1.05543} , { 2.6898, 215.686, -0.658691, -232.478, -0.65869} , { -2.61111, 228.671, -1.38904, -221.746, -1.53506} , { 0.795911, 228.856, -1.975, -225.925, -1.80148} , { -0.222993, 219.339, -1.33694, -232.208, -1.44088} , { -0.093059, 226.935, -0.96251, -219.306, -0.855086} , { -1.66524, 220.611, -0.822399, -231.671, -0.942479} , { -2.78885, 235.605, -1.7154, -219.443, -1.43742} , },},{ { { -6.38707, 1765.66, -1.41963, -1844.26, -1.48181} , { 3.86115, 226.102, -1.19586, -222.923, -1.08131} , { -2.62574, 222.278, -1.30912, -228.914, -1.44272} , { 0.51772, 225.42, -1.28017, -224.911, -1.22406} , { -1.66358, 17.6336, -0.755425, -90.6772, -4.03838} , { 0.97403, 222.8, -0.97259, -225.976, -0.924671} , { -2.61598, 222.897, -0.938596, -228.823, -1.04103} , { 2.3216, 226.557, -1.28542, -223.212, -1.15403} , },{ { -0.911372, 216.444, -1.27195, -235.862, -1.45053} , { 1.71743, 224.443, -1.30787, -226.019, -1.19722} , { -0.76583, 942.44, -1.58182, -970.88, -1.64645} , { 1.31991, 227.685, -1.286, -221.983, -1.16763} , { -1.24747, 221.074, -1.18634, -230.757, -1.34534} , { 1.75494, 231.54, -1.20668, -216.303, -1.05761} , { -1.41044, 220.812, -1.02474, -231.069, -1.16368} , { 1.14971, 233.018, -1.44677, -217.41, -1.26758} , },{ { -1.75279, 225.448, -1.18655, -225.457, -1.288} , { 2.01222, 226.703, -1.27171, -222.979, -1.15314} , { -1.21371, 222.837, -1.35124, -228.366, -1.4912} , { 1.56716, 233.265, -1.48757, -217.389, -1.31098} , { -2.54646, 225.489, -0.867887, -225.64, -0.932951} , { 1.56691, 228.123, -1.21412, -220.794, -1.09921} , { -1.10247, 219.798, -1.57581, -231.688, -1.79849} , { 2.54603, 221.048, -0.989239, -227.641, -0.931508} , },{ { -1.10405, 210.139, -1.46343, -244.742, -1.77612} , { 2.43748, 72.3865, -3.41681, -18.2317, -0.605751} , { -1.51491, 221.023, -0.911438, -230.844, -0.998382} , { 2.89559, 223.893, -0.968145, -223.439, -0.883261} , { -1.64823, 221.214, -1.04409, -230.602, -1.15149} , { 2.93871, 225.82, -0.99061, -221.258, -0.893091} , { -0.504367, 232.028, -1.20221, -217.671, -1.2757} , { 2.90495, 230.003, -0.836523, -215.954, -0.731784} , },{ { -2.31114, 227.365, -1.03861, -223.445, -1.15036} , { 2.41984, 95.1153, -5.33189, -17.4685, -0.651637} , { -1.15241, 223.574, -1.20649, -227.582, -1.31184} , { 2.6873, 226.03, -1.20268, -223.114, -1.09255} , { -1.32145, 225.181, -1.05767, -225.785, -1.13715} , { 2.31298, 229.352, -1.32908, -220.372, -1.19266} , { -2.29236, 222.202, -1.00063, -229.509, -1.10926} , { 2.36424, 224.747, -1.0769, -223.613, -0.988918} , },{ { -3.5282, 224.829, -0.90383, -226.311, -1.00002} , { 1.51094, 228.767, -1.54412, -223.383, -1.38829} , { -2.97202, 228.868, -0.771996, -221.577, -0.812049} , { 2.7179, 223.67, -1.23998, -226.367, -1.12342} , { -1.54891, 226.316, -1.37658, -224.212, -1.46898} , { 3.72693, 225.031, -0.995584, -221.612, -0.889084} , { -2.30539, 224.619, -1.03874, -226.503, -1.13078} , { 1.35313, 226.725, -1.7173, -227.109, -1.55898} , },{ { -4.14499, 227.603, -0.786357, -223.21, -0.858873} , { 2.20772, 18.5647, -2.60532, -15.4385, -0.617832} , { -3.46711, 227.511, -0.75865, -223.32, -0.811921} , { 3.07143, 224.618, -0.939503, -222.645, -0.861271} , { -3.49119, 222.955, -0.946961, -228.711, -1.06722} , { 3.01431, 223.415, -0.742329, -223.581, -0.690361} , { -2.95303, 218.398, -0.982271, -234.933, -1.12773} , { 1.35309, 228.375, -1.83753, -227.243, -1.6317} , },{ { -1.93029, 222.904, -1.03861, -228.336, -1.17836} , { 1.12709, 228.064, -1.60167, -223.989, -1.47961} , { -1.13442, 227.451, -1.00077, -222.918, -1.03516} , { 1.75555, 224.143, -0.720701, -224.291, -0.692396} , { -2.46369, 225.85, -1.0022, -225.075, -1.08274} , { -0.345161, 226.132, -1.02001, -222.404, -0.956868} , { 0.715912, 219.442, -1.14306, -232.537, -1.29185} , { -2.78633, 225.085, -1.42037, -227.496, -1.26018} , },},{ { { -2.60883, 224.337, -1.40492, -226.22, -1.55519} , { 3.79, 230.121, -1.195, -215.123, -1.05065} , { -3.08331, 213.405, -1.46472, -240.512, -1.77481} , { 5.2067, 804.357, -1.93539, -747.094, -1.78387} , { -8.1654, 229.798, -0.616982, -220.543, -0.686652} , { 2.81701, 123.265, -7.26523, -16.2166, -0.600496} , { -1.25038, 221.961, -1.23158, -229.61, -1.37728} , { 2.01475, 225.412, -1.17729, -223.646, -1.07379} , },{ { -2.2573, 229.225, -0.104209, -221.146, -0.09786} , { 1.95414, 224.589, -1.03249, -223.963, -0.963318} , { -1.55058, 224.821, -1.21195, -226.135, -1.304} , { 0.545925, 228.086, -1.59981, -224.03, -1.46542} , { -1.86003, 222.378, -1.11678, -229.21, -1.23853} , { 0.86565, 226.971, -1.35563, -223.471, -1.25287} , { -1.6233, 218.868, -1.33382, -232.989, -1.50545} , { 1.55372, 225.489, -1.15507, -223.726, -1.07383} , },{ { -1.47817, 222.635, -1.42412, -228.55, -1.58439} , { 3.4981, 215.491, -0.7928, -232.282, -0.765305} , { -1.31446, 221.612, -1.56164, -229.651, -1.77362} , { 1.98797, 224.311, -1.12318, -224.505, -1.03147} , { -1.63099, 224.736, -1.29825, -226.307, -1.4419} , { 1.91809, 227.038, -1.14538, -221.404, -1.0375} , { -1.6853, 217.266, -1.17385, -235.196, -1.34853} , { 0.458547, 228.541, -1.38702, -221.874, -1.25286} , },{ { -1.60793, 225.065, -0.916617, -226.171, -1.00749} , { 3.67834, 215.68, -0.41411, -229.937, -0.414109} , { -1.61574, 220.023, -1.36279, -231.829, -1.55359} , { 2.74604, 218.749, -0.817824, -224.709, -0.772087} , { 0.627665, 221.874, -1.01069, -229.893, -1.13153} , { 0.574993, 222.28, -1.14939, -227.349, -1.07312} , { -1.82642, 220.726, -1.17751, -231.045, -1.31374} , { 1.1534, 227.215, -1.18954, -221.316, -1.06212} , },{ { -1.73329, 225.856, -1.29913, -225.076, -1.45511} , { 4.7482, 220.367, -0.791581, -224.681, -0.713808} , { -1.86862, 224.315, -1.10141, -226.545, -1.21285} , { 1.79336, 227.373, -1.28845, -222.296, -1.15971} , { -1.5334, 223.473, -1.21047, -227.726, -1.32373} , { 2.54803, 223.067, -1.09194, -225.815, -1.01134} , { -2.05748, 222.812, -1.0907, -228.658, -1.20159} , { 3.39996, 221.31, -0.871281, -225.865, -0.80859} , },{ { -3.27812, 225.19, -0.823215, -226.047, -0.905632} , { 2.33084, 224.514, -1.52432, -228.412, -1.38237} , { -0.764462, 224.552, -1.46616, -226.189, -1.58402} , { 2.31699, 224.681, -1.14854, -224.431, -1.05986} , { -2.39156, 226.856, -0.951984, -223.911, -1.0227} , { 2.74236, 224.936, -1.27121, -225.143, -1.12766} , { -1.60253, 224.285, -1.20597, -226.764, -1.31316} , { 2.07611, 224.481, -1.27114, -225.707, -1.16843} , },{ { -1.98271, 224.732, -0.920454, -226.41, -0.986612} , { 6.63079, 35.7484, -3.2716, -16.0726, -0.286612} , { -1.54135, 220.974, -1.35441, -230.556, -1.5239} , { 2.1892, 225.265, -1.12724, -223.738, -1.04704} , { -1.80035, 224.594, -1.17088, -226.43, -1.27457} , { 1.80205, 228.707, -1.5126, -222.966, -1.37015} , { -1.25572, 220.873, -1.30921, -230.63, -1.45686} , { 5.52119, 217.076, -0.510687, -228.267, -0.487367} , },{ { -0.597878, 223.911, -1.00885, -227.014, -1.05396} , { 4.27705, 219.588, -0.942211, -228.063, -0.874841} , { -1.55622, 223.107, -1.15155, -228.132, -1.24101} , { 1.86067, 225.516, -1.16795, -223.917, -1.09481} , { -0.968742, 223.498, -1.48328, -227.336, -1.59333} , { 1.06456, 228.122, -1.57557, -224.001, -1.43923} , { 2.2248, 221.768, -1.42353, -229.433, -1.56903} , { 0, 222.746, -1.21728, -228.784, -1.37889} } } }; TF1* HParticleTool::gf1 = 0; TF1* HParticleTool::gf2 = 0; TF1* HParticleTool::gfsum = 0; TF1* HParticleTool::gScaledy=0; TF1* HParticleTool::fdxoffset=0; TF1* HParticleTool::fdxsigma =0; Double_t HParticleTool::scaleDyPars[4] = {3.5,35./20.,9.,2.}; // a + slope + part2 + gain HParticleTool::HParticleTool() { } HParticleTool::~HParticleTool() { } Float_t HParticleTool::phiSecToLabDeg(Int_t sec, Float_t phiRad) { // Convert phi (rad) angle from sec to Lab phi in deg phiRad *= TMath::RadToDeg(); if(sec < 5) return phiRad += 60.0f * sec; else return phiRad -= 60.0f; } Float_t HParticleTool::thetaToLabDeg(Float_t thetaRad) { // Convert theta angle (rad) to the coordinate lab system return TMath::RadToDeg() * thetaRad; } Float_t HParticleTool::phiLabToPhiSecDeg(Float_t phiLabDeg) { // Convert phi angle (deg,lab: 0-360) to the coordinate sector system (60,120) return fmod(phiLabDeg,60.F) + 60; } Int_t HParticleTool::phiLabToSec(Float_t phiLabDeg) { // Convert phi angle (deg,lab: 0-360) to the sector number Int_t sec; if(phiLabDeg >= 60) sec = ((Int_t) phiLabDeg/60) - 1; else sec = 5; return sec; } Float_t HParticleTool::calcRichQA(HMdcSeg* seg, HRichHit* richHit) { // return -1 if fails. if(!seg || ! richHit) return -1; Float_t rphi = fmod(richHit->getPhi(),60.F) + 60; Float_t dTheta = richHit->getTheta() - seg->getTheta()*TMath::RadToDeg(); Float_t dPhi = ( rphi - seg->getPhi()*TMath::RadToDeg() ) * TMath::Sin(seg->getTheta()); return sqrt(dPhi*dPhi + dTheta*dTheta); } Float_t HParticleTool::calcRichQA(HMdcSeg* seg,Float_t richTheta,Float_t richPhi) { // return -1 if fails. richTheta and richPhi in deg in lab system // if(!seg ) return -1; Float_t rphi = fmod(richPhi,60.F) + 60; Float_t dTheta = richTheta - seg->getTheta()*TMath::RadToDeg(); Float_t dPhi = ( rphi - seg->getPhi()*TMath::RadToDeg() ) * TMath::Sin(seg->getTheta()); return sqrt(dPhi*dPhi + dTheta*dTheta); } Float_t HParticleTool::getOpeningAngle(Float_t phi1,Float_t theta1,Float_t phi2,Float_t theta2) { // phi and theta angles of particle 1 and 2 in lab coordinates [deg] // returns opening angle [deg] //Temporary vector objects for scalar product computation TVector3 vec1; TVector3 vec2; // convert angles to radians phi1 *= TMath::DegToRad() ; theta1 *= TMath::DegToRad() ; phi2 *= TMath::DegToRad() ; theta2 *= TMath::DegToRad() ; //Above all angles are in degree! We need them in radians! vec1.SetMagThetaPhi(1.0,theta1,phi1); vec2.SetMagThetaPhi(1.0,theta2,phi2); return TMath::RadToDeg() * vec1.Angle(vec2); } Float_t HParticleTool::getOpeningAngle(TLorentzVector& vec1,TLorentzVector& vec2) { // returns opening angle [deg] return vec1.Angle(vec2.Vect())*TMath::RadToDeg(); } Float_t HParticleTool::getOpeningAngle(HParticleCand* cand1,HParticleCand* cand2) { // returns opening angle [deg] return HParticleTool::getOpeningAngle(cand1->getPhi(),cand1->getTheta(),cand2->getPhi(),cand2->getTheta()); } Float_t HParticleTool::getOpeningAngle(HGeantKine* kine1,HGeantKine* kine2) { // returns opening angle [deg] Float_t theta1,theta2,phi1,phi2; kine1->getPhiThetaDeg(theta1,phi1); kine2->getPhiThetaDeg(theta2,phi2); return HParticleTool::getOpeningAngle(phi1,theta1,phi2,theta2); } Bool_t HParticleTool::setCloseCandidates(Float_t oACut, Bool_t sameSector, Bool_t skipSameSeg) { // loops over HParticleCand category and finds the closest partners of each candidate. set // setAngleToNearbyFittedInner() and setAngleToNearbyUnfittedInner() of the candidates // convention : oAngle is positive if the closest partner had a ring match otherwise negative. // Only opening Angles smaller oACut [Deg] are taken into account. Default opening Angles // are set to -99. if sameSector = kTRUE (default) closest partners are searched in the // same sector only. if skipSameSeg = kTRUE (default) candidates which share the same // inner segment are ignored. HCategory* catCand = HCategoryManager::getCategory(catParticleCand); if(catCand){ UInt_t size = catCand->getEntries(); if(size < 2) return kTRUE; Float_t oAngle = -1; Float_t oAngleSel = -1; Float_t oASmallFitted = -99; Float_t oASmallUnfitted = -99; HParticleCand* cand = 0; HParticleCand* cand2 = 0; //---------------------------------------------------------- for(UInt_t i = 0; i < size ; i ++){ // outer loop oASmallFitted = -99; oASmallUnfitted = -99; cand = HCategoryManager::getObject(cand,catCand,i); cand->setAngleToNearbyFittedInner (oASmallFitted); cand->setAngleToNearbyUnfittedInner(oASmallUnfitted); //---------------------------------------------------------- for(UInt_t j = 0 ; j < size; j ++){ // inner loop if(i==j) continue; cand2 = HCategoryManager::getObject(cand2,catCand,j); if(sameSector && cand->getSector() != cand2->getSector()) continue; // only in same sector if(skipSameSeg && cand->getInnerSegInd() == cand2->getInnerSegInd()) continue; // don't take into account tracks which share same inner seg oAngle = HParticleTool::getOpeningAngle(cand,cand2); // opening angle beteween candiates without sign if(oAngle < oACut){ // only for small angle if(cand2->getInnerSegmentChi2()<=0){ // no fitted inner segment oAngleSel = oAngle < fabs(oASmallUnfitted) ? oAngle: fabs(oASmallUnfitted); if(cand2->getRingCorr()!=0) { // ring match seg or rk oASmallUnfitted = oAngleSel; } else { oASmallUnfitted = -oAngleSel; } } else { // fitted inner segment oAngleSel = oAngle < fabs(oASmallFitted) ? oAngle: fabs(oASmallFitted); if(cand2->getRingCorr()!=0) { // ring match seg or rk oASmallFitted = oAngleSel; } else { oASmallFitted = -oAngleSel; } } } // oACut } // end inner loop cand->setAngleToNearbyFittedInner (oASmallFitted); cand->setAngleToNearbyUnfittedInner(oASmallUnfitted); } // end outer loop return kTRUE; } return kFALSE; } Int_t HParticleTool::getCloseCandidates(HParticleCand* cand, vector& vcand, vector& vopeninAngle, Float_t oACut, Bool_t sameSector, Bool_t skipSameSeg) { // loops over HParticleCand category and finds the closest partners of cand // and fill them to vcand. the opening angles are filled in same order to vopeningAngles. // Only opening Angles smaller oACut [Deg] are taken into account. // if sameSector = kTRUE (default) closest partners are searched in the // same sector only. if skipSameSeg = kTRUE (default) candidates which share the same // inner segment are ignored. The candidates are sorted by inreasing opening angle. vcand .clear(); vopeninAngle.clear(); if(!cand) return 0; HCategory* catCand = HCategoryManager::getCategory(catParticleCand); if(catCand){ UInt_t size = catCand->getEntries(); if(size < 2) return 0; vector vindex; vector vtmpOA; vector vtmpC; Float_t oAngle = -1; HParticleCand* cand2 = 0; //---------------------------------------------------------- for(UInt_t j = 0 ; j < size; j ++){ cand2 = HCategoryManager::getObject(cand2,catCand,j); if(cand == cand2) continue; // skip same cand if(sameSector && cand->getSector() != cand2->getSector()) continue; // only in same sector if(skipSameSeg && cand->getInnerSegInd() == cand2->getInnerSegInd()) continue; // don't take into account tracks which share same inner seg oAngle = HParticleTool::getOpeningAngle(cand,cand2); // opening angle beteween candiates without sign if(oAngle < oACut){ // only for small angle vtmpC .push_back(cand2); vopeninAngle.push_back(oAngle); } // oACut } // end loop //---------------------------------------------------------- vindex.assign(vtmpC.size(),0); HTool::sort(vopeninAngle.size(),vopeninAngle.data(),vindex.data(),kFALSE,kTRUE); //sort increasing by OA , overwrite the vopeningAngle vcand.assign(vtmpC.size(),0); for(UInt_t i=0;i& vSeg,Float_t oACut, Bool_t sameSector, Bool_t skipSameSeg) { // loops over HParticleCand category and finds the closest partners of cand // and fill the inner seg Ind to vSeg. // Only opening Angles smaller oACut [Deg] are taken into account. // if sameSector = kTRUE (default) closest partners are searched in the // same sector only. if skipSameSeg = kTRUE (default) candidates which share the same // inner segment are ignored. // The seg ind are sorted by inreasing opening angle. vSeg.clear(); if(!cand) return 0; vector vopeninAngle; vectorvcand; HParticleTool::getCloseCandidates(cand,vcand,vopeninAngle,oACut,sameSector,skipSameSeg); for(UInt_t i=0;igetInnerSegInd())==vSeg.end()) vSeg.push_back(vcand[i]->getInnerSegInd()); } return vSeg.size(); } void HParticleTool::getTLorentzVector(HGeantKine* kine, TLorentzVector& vec,Int_t pid) { // fills TLorentzVector vec from HGeantKine (if no id is provide id from // kine will be used). // CAUTION : vec.Phi() will be -pi to pi (-180 to +180 deg) // To be compatibel with the HADES lab system phi (0-360 deg) // one has to use // HParticleTool::getLabPhiDeg(TLorentzVector& vec) if(kine == NULL) { vec.SetXYZM(-1000,-1000,-1000,-1); return; } Int_t id = ( pid != -1 ) ? pid : kine->getID(); Float_t mass = HPhysicsConstants::mass(id); Float_t xmom,ymom,zmom; kine->getMomentum(xmom,ymom,zmom); vec.SetPxPyPzE(xmom,ymom,zmom,TMath::Sqrt(mass*mass + xmom*xmom + ymom*ymom + zmom*zmom)); } void HParticleTool::fillTLorentzVector(TLorentzVector& v,HParticleCand* cand,Float_t mass) { // fill v from cand. cand will be unmodified. // user has to provide the mass. if the candidate has // no reconstructed momentum the vector will be filled // by angles only (enough for opening angle calculation, but invMass // for such vectors will be wrong). FOR PARTICLES WITH CHARGE != 1 // USE fillTLorentzVector(TLorentzVector& v,HParticleCand* cand,Int_t pid,Bool_t correctMom=kTRUE) if(cand->getMomentum() != -1){ v.SetXYZM( cand->getMomentum() * TMath::Sin( TMath::DegToRad() * cand->getTheta()) * TMath::Cos( TMath::DegToRad() * cand->getPhi()), cand->getMomentum() * TMath::Sin( TMath::DegToRad() * cand->getTheta() ) * TMath::Sin( TMath::DegToRad() * cand->getPhi() ), cand->getMomentum() * TMath::Cos( TMath::DegToRad() * cand->getTheta() ), mass ); } else { TVector3 v1; v1.SetMagThetaPhi(1.0,cand->getTheta() * TMath::DegToRad(),cand->getPhi() * TMath::DegToRad()); v.SetVect(v1); } } void HParticleTool::fillTLorentzVector(TLorentzVector& v,HParticleCand* cand,Int_t pid,Bool_t correctMom) { // fill v from cand. cand will be unmodified. // user has to provide the PID. if the candidate has // no reconstructed momentum the vector will be filled // by angles only (enough for opening angle calculation, but invMass // for such vectors will be wrong). correctMom will // switch to correction of momentum (HEnergylossCorrPar needed) // The function takes care about particle charge != 1 when // fill the momenta. Float_t mass = HPhysicsConstants::mass(pid); Float_t mom = cand->getMomentumPID(pid); if(correctMom) mom = cand->getCorrectedMomentumPID(pid); if(mom != -1){ v.SetXYZM( mom * TMath::Sin( TMath::DegToRad() * cand->getTheta()) * TMath::Cos( TMath::DegToRad() * cand->getPhi()), mom * TMath::Sin( TMath::DegToRad() * cand->getTheta() ) * TMath::Sin( TMath::DegToRad() * cand->getPhi() ), mom * TMath::Cos( TMath::DegToRad() * cand->getTheta() ), mass ); } else { TVector3 v1; v1.SetMagThetaPhi(1.0,cand->getTheta() * TMath::DegToRad(),cand->getPhi() * TMath::DegToRad()); v.SetVect(v1); } } Float_t HParticleTool::getLabPhiDeg(TLorentzVector& vec) { // CAUTION : vec.Phi() will be -pi to pi (-180 to +180 deg) // To be compatibel with the HADES lab system phi (0-360 deg) // one has to use shift negative values by 360deg return vec.Phi() < 0 ? vec.Phi()*TMath::RadToDeg()+360 : vec.Phi()*TMath::RadToDeg(); } HGeomVector HParticleTool::getGlobalVertex(Int_t v,Bool_t warn) { // return the global vertex position // v can be Particle::kVertexCluster, Particle::kVertexSegment, // Particle::kVertexParticle. In case of Particle::kVertexSegment // and Particle::kVertexParticle the chi2 of the vertex fit is // checked. If the the fit Particle::kVertexSegment failed the position of // Particle::kVertexCluster is returned instead. If Particle::kVertexParticle // fails Particle::kVertexSegment is used if possible otherwise fallback // is Particle::kVertexCluster. if(!gHades) { HGeomVector vec; cout<<"HParticleTool::getGlobalVertex() : No Hades object found !"<getCurrentEvent()) { HGeomVector vec; cout<<"HParticleTool::getGlobalVertex() : No Event structure found !"<getCurrentEvent()->getHeader()->getVertexCluster().getPos(); } else if (v == Particle::kVertexSegment){ if(gHades->getCurrentEvent()->getHeader()->getVertex().getChi2() < 0){ if(warn) cout<<"HParticleTool::getGlobalVertex() : Vertex option v = "<getCurrentEvent()->getHeader()->getVertexCluster().getPos(); } return gHades->getCurrentEvent()->getHeader()->getVertex().getPos(); } else if (v == Particle::kVertexParticle ){ if(gHades->getCurrentEvent()->getHeader()->getVertexReco().getChi2() < 0){ if(gHades->getCurrentEvent()->getHeader()->getVertex().getChi2() < 0){ if(warn) cout<<"HParticleTool::getGlobalVertex() : Vertex option v = "<getCurrentEvent()->getHeader()->getVertexCluster().getPos(); } if(warn) cout<<"HParticleTool::getGlobalVertex() : Vertex option v = "<getCurrentEvent()->getHeader()->getVertex().getPos(); } return gHades->getCurrentEvent()->getHeader()->getVertexReco().getPos(); } else { cout<<"HParticleTool::getGlobalVertex() : Vertex option v = "<getCurrentEvent()->getHeader()->getVertex().getPos(); } }; Double_t HParticleTool:: getMinimumDistToVertex(HParticleCand* cand,HGeomVector& vertex) { // returns distance of closest approach to vertex if(!cand) return -1; HGeomVector base, dir; HParticleTool::calcSegVector(cand->getZ(),cand->getR(),cand->getPhi(),cand->getTheta(),base,dir); return HParticleTool::calculateMinimumDistanceStraightToPoint(base,dir,vertex); } HGeomVector HParticleTool:: getPointOfClosestApproachToVertex(HParticleCand* cand,HGeomVector& vertex) { // returns distance of closest approach to vertex HGeomVector base, dir; HParticleTool::calcSegVector(cand->getZ(),cand->getR(),cand->getPhi(),cand->getTheta(),base,dir); return HParticleTool::calculatePointOfClosestApproachStraightToPoint(base,dir,vertex); } Double_t HParticleTool::scaledy (Double_t* x, Double_t* par) { // scaling function for dy METAMATCH cut with 1/p [GeV/c] // par[0] a : dy at 1/p=0; // par[1] slope per 1/p unit // par[2] part2 : open window for very slow particles // par[3] gain : slope = slope*gain for 1/p > part2 Double_t one_over_p = x[0]; Double_t a = par[0]; Double_t slope = par[1]; Double_t part2 = par[2]; Double_t gain = par[3]; Double_t dy; if(one_over_p part2 // TF1 is created with initial parameters HParticleTool::scaleDyPars[4] at // first call of the function, dyCut [mm] will be used instead of scaleDyPars[0] // if dyCut is positive. Double_t dy = (dyCut>0) ? dyCut : scaleDyPars[0]; if(gScaledy == NULL){ gScaledy = new TF1("scaleDY",scaledy,0,100,4); gScaledy->SetParameter(0,dy); gScaledy->SetParameter(1,scaleDyPars[1]); gScaledy->SetParameter(2,scaleDyPars[2]); gScaledy->SetParameter(3,scaleDyPars[3]); } else { gScaledy->SetParameter(0,dy); } if(!c) return dy; if(c->getMomentum()<0) return dy; Double_t mom = c->getMomentum() * 0.001; // GeV/c if(mom < 0.01) mom = 0.01; if(mom > 1000) mom = 1000; return gScaledy->Eval(1./mom); } Float_t HParticleTool::getRpcCellHalfWidth(Int_t col,Int_t cell) { // returns the half cell width for rpc cell in column col // return -1 if out of bounds (col 0-5,cell 0-31) if(col<0 || cell < 0) return -1; if(col<6 && cell < 32){ return rpcCellHalfWidth[col*32+cell]; } else { return -1; } } Float_t HParticleTool::getTofCellHalfWidth(Int_t mod,Int_t cell) { // returns the half cell width for tof cell in module mod (addressing like data, not geom!) // return -1 if out of bounds (mod 0-7, cell 0-7) if(mod<0 || cell < 0) return -1; if(mod<8 && cell < 8){ return tofCellHalfWidth[mod*8+cell]; } else { return -1; } } Bool_t HParticleTool::isGoodMetaCell(HParticleCand* c,Double_t bound,Bool_t doScaling) { // returns kTRUE is the candidate matches in dy to the meta cell. // returns kFALSE if the candidate has no meta. If the candidate has shower // only it returns alway kTRUE. If doScaling = kTRUE a scaling with 1/p // of bound is applied to account for multiple scattering of slow particles if(c==0) return kFALSE; if(bound<0) bound = 0; if(doScaling) bound = getScaledDy(c,bound); Int_t sys = c->getSystemUsed(); if(sys < 0) return kFALSE; // no meta hit if(c->isTofHitUsed()||c->isTofClstUsed()){ // TOF :if cluster has to different cell sizes take the bigger one Float_t w1 = getTofCellHalfWidth(c->getMetaModule(0),c->getMetaCell(0)); Float_t w2 = getTofCellHalfWidth(c->getMetaModule(1),c->getMetaCell(1)); Float_t dy = fabs(c->getRkMetaDy()); if(w2isRpcClstUsed()) { // RPC : if cluster has to different cell sizes take the bigger one Float_t w1 = getRpcCellHalfWidth(c->getMetaModule(0),c->getMetaCell(0)); Float_t w2 = getRpcCellHalfWidth(c->getMetaModule(1),c->getMetaCell(1)); Float_t dy = fabs(c->getRkMetaDy()); if(w2isShowerUsed()) { // allways good return kTRUE; } return kFALSE; } Bool_t HParticleTool::normDX(HParticleCand* c) { // This function changes // c->setMetaMatchQuality(metaDxnorm) (quality for dx only, dx unchanged); // RK dx values will be normalized by correct width to alow // for general treatment of cuts (RPC (clustersize 1 and 2, TOF and SIM)). if(!fdxoffset) fdxoffset = new TF1("fdxoffset","[0]+[1]*exp(x*[2])+[3]*exp(x*[4])",0,25); // create only once if(!fdxsigma) fdxsigma = new TF1("fdxsigma" ,"[0]+[1]*x**[2]",0,25); // create only once if(!c) return kFALSE; if(c->getBeta()<0) return kTRUE; // no tof or rpc Float_t metaDxnorm = c->getRkMetaDx(); if(metaDxnorm<-999) return kTRUE; // no RK Int_t sector = c->getSector(); Int_t module = c->getMetaModule(0); Int_t cell = c->getMetaCell(0); if(cell<0||module<0) { cout<<"module or cell ==-1, should not happen"< (c); if(c->isTofHitUsed()==1||c->isTofClstUsed()==1) { Bool_t isGeant= kFALSE; if(csim){ // we are in simulation // now find out if it was embedding real hit or pure sim for(Int_t i=0;i<4;i++){ Int_t tr = csim->getGeantTrackMeta(i); if(tr==-1) continue; if(tr > 0) isGeant = kTRUE; } } if(isGeant){ metaDxnorm /= TOFSIGSIM; // fixed width in sim } else { Float_t dedx =c->getTofdEdx(); if(dedx>25) dedx = 25; if(dedx<0) { cout<<"dedx out of bounce "< SetParameters( parsDX[sector][module][cell] ); fdxsigma -> SetParameters( parsSX[sector][module][cell] ); metaDxnorm -= fdxoffset -> Eval(c->getTofdEdx()) ; metaDxnorm /= fdxsigma -> Eval(c->getTofdEdx()) ; } } else if(c->isRpcClstUsed()==1) { Float_t sigma = RPCSIG1; if(c->getMetaModule(1)!=-1) sigma = RPCSIG2 ; // clustersize 2 : sigma/sqrt(2) if(csim){ Bool_t isGeant= kFALSE; for(Int_t i=0;i<4;i++){ Int_t tr = csim->getGeantTrackMeta(i); if(tr==-1) continue; if(tr > 0) isGeant = kTRUE; } if(isGeant){ if(c->getMetaModule(1)!=-1) sigma = RPCSIG2SIM; else sigma = RPCSIG1SIM; } } metaDxnorm /= sigma; } //c->setRkMetaDx(metaDxnorm); c->setMetaMatchQuality(fabs(metaDxnorm)); return kTRUE; } Float_t HParticleTool::getNormDX(HParticleCand* c) { // This function does not change the candidate. // it returns normalized RK dx value (including sign) to alow // for general treatment of cuts (RPC (clustersize 1 and 2, TOF and SIM)). if(!fdxoffset) fdxoffset = new TF1("fdxoffset","[0]+[1]*exp(x*[2])+[3]*exp(x*[4])",0,25); // create only once if(!fdxsigma) fdxsigma = new TF1("fdxsigma" ,"[0]+[1]*x**[2]",0,25); // create only once if(!c) return kFALSE; if(c->getBeta()<0) return kTRUE; // no tof or rpc Float_t metaDxnorm = c->getRkMetaDx(); if(metaDxnorm<-999) return metaDxnorm; // no RK Int_t sector = c->getSector(); Int_t module = c->getMetaModule(0); Int_t cell = c->getMetaCell(0); if(cell<0||module<0) { cout<<"module or cell ==-1, should not happen"< (c); if(c->isTofHitUsed()==1||c->isTofClstUsed()) { Bool_t isGeant= kFALSE; if(csim){ // we are in simulation // now find out if it was embedding real hit or pure sim for(Int_t i=0;i<4;i++){ Int_t tr = csim->getGeantTrackMeta(i); if(tr==-1) continue; if(tr > 0) isGeant = kTRUE; } } if(isGeant){ metaDxnorm /= TOFSIGSIM; // fixed width in sim } else { Float_t dedx =c->getTofdEdx(); if(dedx>25) dedx = 25; if(dedx<0) { cout<<"dedx out of bounce "< SetParameters( parsDX[sector][module][cell] ); fdxsigma -> SetParameters( parsSX[sector][module][cell] ); metaDxnorm -= fdxoffset -> Eval(c->getTofdEdx()) ; metaDxnorm /= fdxsigma -> Eval(c->getTofdEdx()) ; } } else if(c->isRpcClstUsed()==1) { Float_t sigma = RPCSIG1; if(c->getMetaModule(1)!=-1) sigma = RPCSIG2 ; // clustersize 2 : sigma/sqrt(2) if(csim){ Bool_t isGeant= kFALSE; for(Int_t i=0;i<4;i++){ Int_t tr = csim->getGeantTrackMeta(i); if(tr==-1) continue; if(tr > 0) isGeant = kTRUE; } if(isGeant){ if(c->getMetaModule(1)!=-1) sigma = RPCSIG2SIM; else sigma = RPCSIG1SIM; } } metaDxnorm /= sigma; } return metaDxnorm; } Float_t HParticleTool::getSigmaDX(HParticleCand* c) { // This function does not change the candidate. // it returns sigma RK dx value to alow // for general treatment of cuts (RPC (clustersize 1 and 2, TOF and SIM)). if(!fdxoffset) fdxoffset = new TF1("fdxoffset","[0]+[1]*exp(x*[2])+[3]*exp(x*[4])",0,25); // create only once if(!fdxsigma) fdxsigma = new TF1("fdxsigma" ,"[0]+[1]*x**[2]",0,25); // create only once if(!c) return kFALSE; if(c->getBeta()<0) return kTRUE; // no tof or rpc Float_t sigma = -1; Float_t metaDxnorm = c->getRkMetaDx(); if(metaDxnorm<-999) return sigma; // no RK Int_t sector = c->getSector(); Int_t module = c->getMetaModule(0); Int_t cell = c->getMetaCell(0); if(cell<0||module<0) { cout<<"module or cell ==-1, should not happen"< (c); if(c->isTofHitUsed()==1||c->isTofClstUsed()) { Bool_t isGeant= kFALSE; if(csim){ // we are in simulation // now find out if it was embedding real hit or pure sim for(Int_t i=0;i<4;i++){ Int_t tr = csim->getGeantTrackMeta(i); if(tr==-1) continue; if(tr > 0) isGeant = kTRUE; } } if(isGeant){ return TOFSIGSIM; // fixed width in sim } else { Float_t dedx =c->getTofdEdx(); if(dedx>25) dedx = 25; if(dedx<0) { cout<<"dedx out of bounce "< SetParameters( parsDX[sector][module][cell] ); fdxsigma -> SetParameters( parsSX[sector][module][cell] ); sigma = fdxsigma -> Eval(c->getTofdEdx()) ; } } else if(c->isRpcClstUsed()==1) { sigma = RPCSIG1; if(c->getMetaModule(1)!=-1) sigma = RPCSIG2 ; // clustersize 2 : ~ sigma/sqrt(2) if(csim){ Bool_t isGeant= kFALSE; for(Int_t i=0;i<4;i++){ Int_t tr = csim->getGeantTrackMeta(i); if(tr==-1) continue; if(tr > 0) isGeant = kTRUE; } if(isGeant){ if(c->getMetaModule(1)!=-1) sigma = RPCSIG2SIM; else sigma = RPCSIG1SIM; } } } return sigma; } Float_t HParticleTool::getCorrectedMomentum(HParticleCand* c) { // returns the corrected momenta (not changing the candidate) // corrections: 1. systematic deviation for high momenta (obtained // from simulation, applied for SIM+REAL) 2. correction for field map // (obtained from apr12 data applied for REAL only). // REMARK: NO ENERGYLOSS CORRECTION HERE! if(!c) return 0; Float_t dmom = c->getMomentum(); if(dmom < 0) return dmom; Float_t theta = c->getTheta(); //-------------------------------------------------- // do correction for systematics in mom at high mom // (SIM + REAL DATA) static TF1* ftheta_mom_sys = 0; static TF1* fpot_mom_sys = 0; if(!ftheta_mom_sys){ ftheta_mom_sys = new TF1("ftheta_mom_sys","[0]/([1]-exp([2]*x+[3]))",0,90); ftheta_mom_sys -> SetParameters(3.33755e-02,-5.96582e+00,-7.78283e-02,5.42673e+00); } if(!fpot_mom_sys) { fpot_mom_sys = new TF1("fpot_mom_sys","[1]*x**[0]+[2]",0,100); fpot_mom_sys -> SetParameters(-3.,-0.003,0.001); } Float_t par = ftheta_mom_sys -> Eval(theta); fpot_mom_sys -> SetParameter(1,par); Float_t dev = fpot_mom_sys -> Eval(1./dmom*1000.); if(fabs(dev)>0.05) return dmom; dmom *= (1.+dev); //-------------------------------------------------- HParticleCandSim* csim = dynamic_cast (c); if(csim){ // we are in simulation Int_t tr =-1; Bool_t found=kFALSE; // for embedding we check who creates the track for(Int_t i=0;i<2;i++){ tr = csim->getGeantTrackInnerMdc(i); if(tr>0) {found = kTRUE; break;} tr = csim->getGeantTrackOuterMdc(i); if(tr>0) {found = kTRUE; break;} } if(found) return dmom; // track is from GEANT } //-------------------------------------------------- // field scaling (REAL DATA ONLY) Float_t moms = 0.954297 + theta* 0.00584318 + theta*theta * -0.000218986 + theta*theta*theta * 3.29399e-06 + theta*theta*theta*theta * -1.69156e-08; dmom *= moms; //-------------------------------------------------- return dmom; } Float_t HParticleTool::setCorrectedMomentum(HParticleCand* c) { // returns the corrected momenta (candidate is changed) // corrections: 1. systematic deviation for high momenta (obtained // from simulation, applied for SIM+REAL) 2. correction for field map // (obtained from apr12 data applied for REAL only). // REMARK: NO ENERGYLOSS CORRECTION HERE! if(!c) return 0; Float_t dmom = HParticleTool::getCorrectedMomentum(c); c->setMomentum(dmom); return dmom; } Bool_t HParticleTool::isParticledEdx(Int_t PID, HParticleCand* pCand, Float_t& deloss, Float_t& dsigma) { // return the delta eloss and sigma of particle assuming PID from // the calibrated Mdc dedx vs beta curve. The functions takes // into account the dependency of beta from energyloss (as function // of theta .charge =1 and charge =2 are treated correct. // The function makes use of HEnergyLossPar. You have to init // the container before. Params are taken from apr12. if(!pCand) return kFALSE; if(pCand->getBeta() < 0 ) return kFALSE; if(pCand->getMdcdEdx() <= 0 ) return kFALSE; if(HPhysicsConstants::charge(PID) == 0 ) return kFALSE; if((HPhysicsConstants::charge(PID)/pCand->getCharge())<0) return kFALSE; static TF1* fdedx_vs_beta = 0; if(!fdedx_vs_beta) fdedx_vs_beta = new TF1("f_dedx_vs_beta",HParticleTool::dedxfunc,0.15,1.5,7); Float_t momQscale = 1.; if( HPhysicsConstants::charge(PID) == 2 ) momQscale = 2.; Float_t theta = pCand->getTheta(); Float_t distMeta = pCand->getDistanceToMetaHit(); Float_t beta = pCand->getBeta(); Float_t dedx = pCand->getMdcdEdx(); fdedx_vs_beta -> SetParameter(0,PID); fdedx_vs_beta -> SetParameter(1,1.); fdedx_vs_beta -> SetParameter(2,0.); fdedx_vs_beta -> SetParameter(3,1.); fdedx_vs_beta -> SetParameter(4,1); fdedx_vs_beta -> SetParameter(5,theta); fdedx_vs_beta -> SetParameter(5,(0.0783057+theta*0.00145768)*2100./distMeta); Float_t meandEdx = fdedx_vs_beta -> Eval(beta); fdedx_vs_beta -> SetParameter(1,0.4); Float_t mindEdx = fdedx_vs_beta -> Eval(beta); fdedx_vs_beta -> SetParameter(1,2.05); Float_t maxdEdx = fdedx_vs_beta -> Eval(beta); if(dedx>mindEdx && dedxgetCharge()== 0) return kFALSE; if(pCand->getBeta() < 0) return kFALSE; Float_t chrg=HPhysicsConstants::charge(PID); if(chrg== 0) return kFALSE; Int_t sys = pCand->getSystemUsed(); if(sys < 0) return kFALSE; if((chrg/pCand->getCharge())<0) return kFALSE; Float_t momQscale = 1.; if( chrg == 2 ) momQscale = 2.; if(pCand->getCharge()>0 && (PID == 9) ) { return kFALSE; } if(pCand->getCharge()<0 && (PID == 8 || PID == 14) ) { return kFALSE; } Float_t theta = pCand->getTheta(); Float_t mom = pCand->getMomentum(); Float_t b = pCand->getBeta(); Float_t maxbeta=-1.0; Float_t minbeta=-1.0; if(mommomMax) return kFALSE; static TF1* f_beta_mom = 0; if(!f_beta_mom) f_beta_mom = new TF1("f_beta_mom", HParticleTool::betaandgammafunc,0,10000,7); Float_t m = HPhysicsConstants::mass(PID); // mom/be*sqrt(1.-be*be); Float_t distMeta = pCand->getDistanceToMetaHit(); f_beta_mom->SetParameter(0,PID); f_beta_mom->SetParameter(1,1.); f_beta_mom->SetParameter(2,0.); f_beta_mom->SetParameter(3,1.00*momQscale); f_beta_mom->SetParameter(4,0); f_beta_mom->SetParameter(5,theta); f_beta_mom->SetParameter(6,(0.0783057+theta*0.00145768)*2100./distMeta); // Will be set to a constant overall! 0.1 Float_t be = f_beta_mom->Eval(mom); Float_t sT = 0.15; if(sys == 0){ sT = 1./sqrt(2.)*(0.056+be*be*be*0.353-be*be*be*be*0.287); if(pCand->getMetaModule(0)>=0 && pCand->getMetaCell(1)>=0 ) sT /= sqrt(2.); // cluster SYS0 } if(sys==1) { sT = HParticleTool::getSigmaDX(pCand)/144.; } Float_t sMom = 0.01; // sigma Mom 1% sMom = 1.25978e-02 * pow(1./mom*1000.,-9.41923e-01) + 1.45266e-02; // sigma Mom 1% Float_t mT = pCand->getDistanceToMetaHit()/b/299.792458; // time of flight used for beta uncertanty calculation dtime = distMeta/b/(TMath::C()*1e-6) - distMeta/be/(TMath::C()*1e-6); Float_t dsigmab = sqrt( pow((sT/((mT/mom*sqrt(mom*mom+m*m) )-sT)),2) + pow(( (mom*(1.+sMom))/sqrt((mom*(1.+sMom))*(mom*(1.+sMom))+m*m) -mom/sqrt(mom*mom+m*m) ),2))*1.; Float_t sigbeta2 = sqrt( pow( m*m/mom/mom/mom*pow(m*m/mom/mom+1.,-3/2.) * sMom*mom,2) + pow(distMeta/mT/mT/(TMath::C()*1e-6)*sT,2) ); dsigma = distMeta/be/be/(TMath::C()*1e-6)*dsigmab; Float_t dsigma2 = distMeta/be/be/(TMath::C()*1e-6)*sigbeta2; dsigma = dsigma2; if(PID==2 || PID==3) { dsigma = sT; // for electrons no mom dependence } if(PID!=2 && PID!=3) { maxbeta = be + nsigma*sigbeta2; minbeta = be - nsigma*sigbeta2; } else { if(sys==0) { minbeta = 0.97; maxbeta = 1.04; } if(sys==1) { minbeta = 0.95; maxbeta = 1.1; } } if(b>minbeta && b< maxbeta) return kTRUE; else return kFALSE; } Bool_t HParticleTool::correctPathLength(HParticleCand* pCand, HGeomVector& vertex, const HMdcPlane* planes, const HGeomVector& targetMidPoint, Double_t beamEnergy ) { if(!pCand) return kFALSE; if(pCand->getBeta()>0) return kFALSE; if(pCand->getChi2()<0) return kFALSE; if(vertex.X()==-1000. && vertex.Y()==-1000. && vertex.Z()==-1000. ) return kFALSE; const Double_t c = 299.792458 ; Double_t beamVelocity = HParticleTool::beta(14,HParticleTool::kinEToMom(14,beamEnergy)); Double_t xVm = targetMidPoint.getX(); Double_t yVm = targetMidPoint.getY(); Double_t zVm = targetMidPoint.getZ(); Float_t mom = pCand->getMomentum(); Float_t newPath = pCand->getDistanceToMetaHit(); Float_t newBeta = pCand->getBeta(); Float_t tof = newPath/newBeta/c; //in (ns) Double_t xP, yP, zP; planes[pCand->getSector()].calcSegIntersec(pCand->getZ(), pCand->getR(), pCand->getTheta()*TMath::DegToRad(), pCand->getPhi()*TMath::DegToRad(), xP, yP, zP ); Float_t rktl = sqrt((xP-xVm)*(xP-xVm)+(yP-yVm)*(yP-yVm)+(zP-zVm)*(zP-zVm) ); // And the real was: Float_t rktl_real = sqrt( (xP-vertex.getX())*(xP-vertex.getX()) + (yP-vertex.getY())*(yP-vertex.getY()) + (zP-vertex.getZ())*(zP-vertex.getZ()) ); // And the new PL is....: newPath -= rktl ; newPath += rktl_real; tof += (zVm-vertex.getZ()) /beamVelocity/c; newBeta = newPath/tof/c; Float_t newMass2 = mom*mom*(1.-newBeta*newBeta)/newBeta/newBeta; HParticleCandSim* csim = dynamic_cast (pCand); if(csim){ // we are in simulation Int_t tr =-1; Bool_t found=kFALSE; // for embedding we check who creates the track for(Int_t i=0;i<2;i++){ tr = csim->getGeantTrackInnerMdc(i); if(tr>0) {found = kTRUE; break;} tr = csim->getGeantTrackOuterMdc(i); if(tr>0) {found = kTRUE; break;} } if(!found) { // sim , but real track from embedding //Tof correction (ToF of the incident Ion to the center of the target. pCand->setMass2(newMass2); pCand->setBeta(newBeta); } else { // simulation : here we do not correct the beam ion tof = pCand->getTof(); newBeta = newPath/tof/c; newMass2 = mom*mom*(1.-newBeta*newBeta)/newBeta/newBeta; pCand->setMass2(newMass2); pCand->setBeta(newBeta); } } else { // for real data we correct tof+beta+mass //Tof correction (ToF of the incident Ion to the center of the target. pCand->setMass2(newMass2); pCand->setBeta(newBeta); } pCand->setDistanceToMetaHit(newPath); return kTRUE; } Bool_t HParticleTool::checkCropedLayer(HGeantKine* kine,HMdcLayer* mdcLayer, Bool_t* croped) { // this function checks for particle kine (and his decayed daughter it it exist) // if the track is in the croped area of the layer. if the track has a hit // but not in the layers to check (0 or 5, depending on first or last module in segment) // it will be counted as croped. The function returns true if anay of the mdcs // was croped. if an pointer to Bool_t croped[4] is provided the crop results // per module will be returned to the array. if(!kine) return kFALSE; if(!mdcLayer) return kFALSE; //------------------------------------------------------- // has the particle a decayed daughter ? pion->muon vector daughters; HGeantKine::getDaughters(kine,daughters); HGeantKine* d = 0; for(UInt_t i=0;igetMechanism()==5 && d1->getFirstMdcHit()!=-1) { d = d1; break;} } //------------------------------------------------------- //------------------------------------------------------- // check layer hits of particle (and + decay daugher when it exist) vector v; kine->getMdcHits(v); Int_t mLay [4] ={0,5,0,5}; // first /last lay of module i Float_t ax, ay, atof, ptof; Int_t s,m,l; Bool_t isCroped [4] = {kFALSE,kFALSE,kFALSE,kFALSE}; // is croped Bool_t isChecked [4] = {kFALSE,kFALSE,kFALSE,kFALSE}; // is in module Bool_t isCheckedL[4] = {kFALSE,kFALSE,kFALSE,kFALSE}; // is in layer to check TVector2 p; for(UInt_t i=0;igetSector(); m = mdc->getModule(); l = mdc->getLayer(); isChecked[m] = kTRUE; if(l==mLay[m]){ isCheckedL[m] = kTRUE; mdc->getHit(ax,ay,atof,ptof); p.Set(ax,ay); if(mdcLayer->isCroped(p,s,m,l)) { isCroped[m] = kTRUE; } else { isCroped[m] = kFALSE; } } } if(d){ // on daugher if exists d->getMdcHits(v); for(UInt_t i=0;igetSector(); m = mdc->getModule(); l = mdc->getLayer(); isChecked[m] = kTRUE; if(l==mLay[m]){ isCheckedL[m] = kTRUE; mdc->getHit(ax,ay,atof,ptof); p.Set(ax,ay); if(mdcLayer->isCroped(p,s,m,l)) { isCroped[m] = kTRUE; } else { isCroped[m] = kFALSE; } } } } //------------------------------------------------------- //------------------------------------------------------- // evaluate the track Bool_t isCropedAll=kFALSE; for(Int_t i=0;i<4;i++ ){ // of MDC if(isChecked[i] && (!isCheckedL[i] || isCroped[i]) ) { isCropedAll =kTRUE; } } if(croped) { // per MDC for(Int_t i=0;i<4;i++ ){ if(isChecked[i]) { croped[i] = (isCroped[i] || !isCheckedL[i]) ? kTRUE : kFALSE; } else { croped[i] = kFALSE; } } } //------------------------------------------------------- return isCropedAll; } //-------------------------------------------------- // kinematic helper functions Double_t HParticleTool::beta(Int_t id,Double_t p) { // Id : see HPysicsConstants // p [MeV/c] // return beta Double_t mass = HPhysicsConstants::mass(id); if(mass <= 0) return -1; return sqrt(1. / (((mass*mass)/(p*p)) + 1.)); } Double_t HParticleTool::betaToP(Int_t id,Double_t beta) { // Id : see HPysicsConstants // return p [MeV/c] Double_t mass = HPhysicsConstants::mass(id); if(mass <= 0) return -1; Double_t p = sqrt( mass*mass / (1./(beta*beta)-1)); return p; } Double_t HParticleTool::gamma(Int_t id,Double_t p) { // Id : see HPysicsConstants // p [MeV/c] // return gamma Double_t mass = HPhysicsConstants::mass(id); if(mass <= 0) return -1; Double_t beta2 = 1. / (((mass*mass)/(p*p)) + 1.); return sqrt(1./ (1 - beta2)); } Double_t HParticleTool::gammaToBeta(Int_t id,Double_t gamma) { // Id : see HPysicsConstants // return beta Double_t mass = HPhysicsConstants::mass(id); if(mass <= 0) return -1; if(gamma*gamma<1) return 0; Double_t beta = sqrt( (-1./(gamma*gamma))+1) ; return beta; } Double_t HParticleTool::gammaToP(Int_t id,Double_t gamma) { // Id : see HPysicsConstants // return p [MeV/c] Double_t mass = HPhysicsConstants::mass(id); if(mass <= 0) return -1; if(gamma*gamma<1) return 0; Double_t beta2 = (1./(gamma*gamma))-1 ; Double_t p = sqrt( mass*mass / (1./(beta2)-1)); return p; } Double_t HParticleTool::betagamma(Int_t id,Double_t p) { // Id : see HPysicsConstants // p [MeV/c] // return beta*gamma Double_t mass = HPhysicsConstants::mass(id); if(mass <= 0) return -1; return p/mass; } Double_t HParticleTool::betagammaToP(Int_t id,Double_t betagamma) { // Id : see HPysicsConstants // return p [MeV/c] Double_t mass = HPhysicsConstants::mass(id); if(mass <= 0) return -1; Double_t p = betagamma*mass; return p; } Double_t HParticleTool::kinEToMom(Int_t id, Double_t Ekin){ // Id : see HPysicsConstants input kinetic energy // return p [MeV/c] Double_t m = HPhysicsConstants::mass(id); return sqrt( pow(Ekin+m,2) - pow(m,2)); } Double_t HParticleTool::momToKinE(Int_t id, Double_t p){ // Id : see HPysicsConstants // return kinetic energy [MeV] Double_t m = HPhysicsConstants::mass(id); return sqrt( pow(p,2)+pow(m,2) ) - m; } Double_t HParticleTool::dedxfunc(Double_t *x, Double_t *par) { // function of MdcDedx used in TF1 // par[0] Id : see HPysicsConstants // par[1] scale return val // par[2] shift x (after scaling) // par[3] scale x // par[4] option: 0:p,1:beta,2:gamma,3:betagamma // par[5] theta [deg]: do dedx corr for fixed theta .if <=0 don't do anything // par[6]; frac of dedx corr // return MDC dEdx Double_t p = fabs(x [0]); Int_t id = par[0]; Double_t scale = par[1]; Double_t offset= par[2]; Double_t scaleX= par[3]; Int_t opt = par[4]; Double_t theta = par[5]; // theta in deg : if <=0 don't do anything Double_t frac = par[6]; // frac of dedx corr p*=scaleX; p+=offset; if (opt==1){ // beta Double_t mass = HPhysicsConstants::mass(id); Double_t beta=p; if(beta>=1) beta=0.999999; p = sqrt( mass*mass / (1./(beta*beta)-1)); } else if(opt==2){ // gamma Double_t mass = HPhysicsConstants::mass(id); Double_t gamma=p; if(gamma*gamma<1) return 0; Double_t beta2 = (-1./(gamma*gamma))+1 ; p = sqrt( mass*mass / (1./(beta2)-1)); } else if(opt==3){ // betagamma Double_t mass = HPhysicsConstants::mass(id); Double_t betagamma=p; p = betagamma*mass; } HEnergyLossCorrPar* corrpar = HEnergyLossCorrPar::getObject(); Double_t dp=0; if(theta>0 && corrpar) { dp = corrpar->getDeltaMom(id,p,theta); } p+=frac*dp; Double_t loss = HMdcDeDx2::energyLoss(id,p,0.6); return scale*loss; } Double_t HParticleTool::betaandgammafunc(Double_t *x, Double_t *par) { // function of beta / gamma used in TF1 // par[0] Id : see HPysicsConstants // par[1] scale return val // par[2] shift x (after scaling) // par[3] scale x // par[4] option: 0:beta,1:gamma,2:betagamma // par[5] theta [deg]: do dedx corr for fixed theta .if <=0 don't do anything // par[6]; frac of dedx corr // return beta or gamma or betagamma Double_t p = fabs(x [0]); Int_t id = par[0]; Double_t scale = par[1]; Double_t offset= par[2]; Double_t scaleX= par[3]; Int_t opt = par[4]; Double_t theta = par[5]; // theta in deg : if <=0 don't do anything Double_t frac = par[6]; // frac of dedx corr Double_t mass = HPhysicsConstants::mass(id); p*=scaleX; HEnergyLossCorrPar* corrpar = HEnergyLossCorrPar::getObject(); Double_t dp=0; if(theta>0 && corrpar) { dp = corrpar->getDeltaMom(id,p,theta); } p+=frac*dp; p+=offset; if (opt==0){ // beta return scale*sqrt(1. / (((mass*mass)/(p*p)) + 1.)); } else if(opt==1){ // gamma return scale*gamma(id,p); } else if(opt==2){ // betagamma return scale*betagamma(id,p); } return 0; } TCutG* HParticleTool::makeCut(TF1* lower,TF1* upper,TString name, Double_t xlow,Double_t xup,Double_t ymin,Double_t ymax,Int_t npoint, Int_t linecolor,Int_t linestyle) { // creates TCutG from upper and lower TF1 in range xlow to xup // the x interval is splitted into npoints. User has to delete // the TCutG objects if(!lower || !upper) return 0; TCutG* cut = new TCutG(); cut->SetName(name.Data()); if(npoint<2) npoint=50; Double_t x, y; Int_t ct = 0 ; for(Int_t i=0;iEval(x); if(y>ymax || ySetPoint(i,x,y); ct++; } for(Int_t i=0;iEval(x); if(y>ymax || ySetPoint(i+npoint,x,y); ct++; } cut->GetPoint(0,x,y); cut->SetPoint(ct,x,y); cut->SetLineStyle(linestyle); cut->SetLineColor(linecolor); return cut; } Double_t HParticleTool::ptyfunc (Double_t* x, Double_t* par) { // return pt as function of y // for fix theta or momenta of a particle Int_t id = par[0]; Double_t M =HPhysicsConstants::mass(id); if(M<=0) return 100000; Double_t midRap= par[2]; Int_t opt = par[3]; Double_t y = x[0]-midRap; if(opt==0){// angles Double_t ThetaLab = par[1]*TMath::DegToRad(); Double_t Part_A = 0.0; Double_t pt = 0.0; Part_A = pow(1/(tan(ThetaLab)*sinh(y)),2.0)-1.0; if(Part_A >= 0.0) { pt = M*pow(Part_A,-0.5); } else return 1000000.0; if(pt > 0.0 && pt < 10000.0) { return pt; } else { return 1000000.0; } return pt; } else if(opt==1){// momenta Double_t PLab = par[1]; Double_t pt = sqrt((pow(PLab,2)-pow(M*sinh(y),2))/(pow(sinh(y),2)+1.0)); if(pt > 0.0 && pt < 10000.0){ return pt; } else { return -1000000.0; } return pt; } else return 100000; } Double_t HParticleTool::fcross(Double_t* xin, Double_t* par){ if(!gf1 || !gf2) return -1000; Double_t x = xin[0]; return gf1->Eval(x) - gf2->Eval(x); } Bool_t HParticleTool::getIntersectionPoint(TF1* f1,TF1* f2,Double_t &xout,Double_t& yout,Double_t xlow,Double_t xup,Int_t n) { // calculates the intersection point of f1 and f2 in interval xlow,xup // to vars xout, yout // return kFALSE if no intersection was found (xout=yout=-10000) // higher number of points n give better precision at the higher // calculation time if(!f1 || !f2) return kFALSE; xout=-10000; yout=-10000; if(gf1){ delete gf1;} gf1 = new TF1(); gf1 ->SetName("cross_f1"); f1->Copy(*gf1); if(gf2){ delete gf2 ;} gf2 = new TF1(); gf2 ->SetName("cross_f2"); f2->Copy(*gf2); if(!gfsum) { gfsum = new TF1("cross_f12",fcross); } gfsum->SetNpx(n); Double_t x = gfsum->GetX(0,xlow,xup); // Double_t epsilon = 1.E-10, Int_t maxiter = 100, Bool_t logx = false) if(x!=xup){ xout=x; yout=gf1->Eval(x); return kTRUE; } ::Warning("HParticleTool::getIntersectionPoint()","No intersection of %s and %s in interval [%f,%f] found!",f1->GetName(),f2->GetName(),xlow,xup); return kFALSE; } TF1* HParticleTool::energyLossTF1(Int_t id, TString name,TString opt , Double_t scaleY,Double_t xoffset,Double_t scaleX,Double_t theta,Double_t frac, Double_t xmin,Double_t xmax, Int_t linecolor,Int_t linestyle,Int_t npoints) { // MdcDedx TF1 // name : name of TF1 (default dedx, the id of the paricle will be appended) // opt : p beta gamma betagamma // // par[0] Id : see HPysicsConstants // par[1] scale return val // par[2] shift x (after scaling) // par[3] scale x // par[4] option: 0:p,1:beta,2:gamma,3:betagamma // par[5] theta [deg]: do dedx corr for fixed theta .if <=0 don't do anything // par[6] fraction of dedx corr // return TF1 // The user has to delete the TF1 object if(HPhysicsConstants::mass(id)<=0) return 0; if(name.CompareTo("")==0) name="dedx"; TF1* f = new TF1(Form("%s_%s",name.Data(),HPhysicsConstants::pid(id)),dedxfunc,xmin,xmax,7); f->SetNpx(npoints); f->SetParameter(0,id); f->SetParameter(1,scaleY); f->SetParameter(2,xoffset); f->SetParameter(3,scaleX); if (opt.CompareTo("p")==0) f->SetParameter(4,0); else if(opt.CompareTo("beta")==0) f->SetParameter(4,1); else if(opt.CompareTo("gamma")==0) f->SetParameter(4,2); else if(opt.CompareTo("betagamma")==0) f->SetParameter(4,3); else { f->SetParameter(4,0);} f->SetParameter(5,theta); f->SetParameter(6,frac); f->SetLineColor(linecolor); f->SetLineStyle(linestyle); return f; } TF1* HParticleTool::betaAndGammaTF1(Int_t id, TString name,TString opt, Double_t scaleY,Double_t xoffset,Double_t scaleX,Double_t theta,Double_t frac, Double_t xmin,Double_t xmax, Int_t linecolor,Int_t linestyle,Int_t npoints) { // beta / gamma TF1 // name : name of TF1 (default equals opt, the id of the paricle will be appended) // opt : beta gamma betagamma // // par[0] Id : see HPysicsConstants // par[1] scale return val // par[2] shift x (after scaling) // par[3] scale x // par[4] option: 0:p,1:beta,2:gamma,3:betagamma // par[5] theta [deg]: do dedx corr for fixed theta .if <=0 don't do anything // par[6] fraction of dedx corr // return TF1 // The user has to delete the TF1 object if(HPhysicsConstants::mass(id)<=0) return 0; if(name.CompareTo("")==0) name=opt; TF1* f = new TF1(Form("%s_%s",name.Data(),HPhysicsConstants::pid(id)),betaandgammafunc,xmin,xmax,7); f->SetNpx(npoints); f->SetParameter(0,id); f->SetParameter(1,scaleY); f->SetParameter(2,xoffset); f->SetParameter(3,scaleX); if (opt.CompareTo("beta")==0) f->SetParameter(4,0); else if(opt.CompareTo("gamma")==0) f->SetParameter(4,1); else if(opt.CompareTo("betagamma")==0) f->SetParameter(4,2); else { f->SetParameter(4,0);} f->SetParameter(5,theta); f->SetParameter(6,frac); f->SetLineColor(linecolor); f->SetLineStyle(linestyle); return f; } TF1* HParticleTool::ptyTF1(Int_t id,Double_t val, TString name, TString opt, Double_t xmin,Double_t xmax,Double_t midRap, Int_t linecolor,Int_t linestyle) { // ptY TF1 // name : name of TF1 (default opt, the d of the paricle will be appended) // opt : theta momentum // // par[0] Id : see HPysicsConstants // par[1] val // par[2] midRap (default 0) // par[3] option: 0:theta,1:momentum // return TF1 // The user has to delete the TF1 object opt.ToLower(); if(HPhysicsConstants::mass(id)<=0) return 0; if(name.CompareTo("")==0) name=opt; TF1* f = new TF1(Form("prty_%s_%s_%5.2f",name.Data(),HPhysicsConstants::pid(id),val),ptyfunc,xmin,xmax,4); f->SetTitle(Form("%5.1f",val)); f->SetNpx(500); f->SetParameter(0,id); f->SetParameter(1,val); f->SetParameter(2,midRap); if (opt.CompareTo("theta")==0) f->SetParameter(3,0); else if(opt.CompareTo("momentum")==0) f->SetParameter(3,1); else { f->SetParameter(3,0);} f->SetLineColor(linecolor); f->SetLineStyle(linestyle); return f; } vector HParticleTool::ptyGrid(Int_t id,vector& vtheta,vector& vmom, TString name, TString opt, Double_t xmin,Double_t xmax,Double_t midRap, Int_t linecolorTheta,Int_t linestyleTheta, Int_t linecolorMom,Int_t linestyleMom) { // draws a pt -y grid for a give particle at given list of theta and momenta // opt : draw drawcopy // range xmin to xmax // The created TF1 objects are return in a vector. The use has to delet them. opt.ToLower(); vector grid; for(UInt_t i=0;iDraw("Lsame"); if(opt.CompareTo("drawcopy")==0)f->DrawCopy("Lsame"); } for(UInt_t i=1;iDraw("Lsame"); if(opt.CompareTo("drawcopy")==0)f->DrawCopy("Lsame"); } return grid; } vector HParticleTool::ptyGrid(Int_t id,TString setup, TString name, TString opt, Double_t xmin,Double_t xmax,Double_t midRap, Int_t linecolorTheta,Int_t linestyleTheta, Int_t linecolorMom,Int_t linestyleMom, TString labels) { // Draw pt - y grid for particle id: lines for constant #theta and constant momentum //------------------------------------ // setup: // theta : format string "@theta:8,0,10" : draw 8 lines starting from 0 deg with stepsize 10 deg // "@theta:-1,10,20,45,60,75" : draw lines at list values (first entry has to be -1) // momentum : format string "@momentum:8,100,200" : draw 8 lines starting from 100 MeV/c with stepsize 200 MeV/c // "@momentum:-1,500,1000,1500" : draw lines at list values (first entry has to be -1) //------------------------------------ // name : optional string for the name construction of TF1 // opt : "" no draw // "draw" draw in current canvas // "drawcopy" draw a copy in current canvas // xmin,xmax : rapidity range (default 0 - 2) // midRap : midRapidity (default 0) // linecolorTheta,lienstyleTheta : line attributes of theta (default black,dashed) // linecolorMom ,lienstyleMom : line attributes of momentum (default black,dashed) // //------------------------------------ // labels="@theta:draw=yes,format=%5.1f#circ,textsize=0.021,angle=0,align=-1@momentum:draw=yes,format=%5.1f MeV/c,textsize=0.023,angle=0,align=-1,yoff=0,02,xoff=0.01 // draw : yes (default) // format: print format // textsize: 0.021 // angle : 0 // align : -1(default: build in) sum of TTextAtributes: kHAlignLeft = 10, kHAlignCenter = 20, kHAlignRight = 30 ,kVAlignBottom = 1, kVAlignCenter = 2, kVAlignTop = 3 // yoff : 0.02 offset from yaxis // xoff : 0.01 offset from xaxis vector grid; setup.ToLower(); opt.ToLower(); TLatex latex; Int_t nTh=0; Double_t strtTh=0; Double_t stepTh=0; Int_t nMom=0; Double_t strtMom=0; Double_t stepMom=0; Bool_t lab_theta = kTRUE; TString format_theta = "%5.1f#circ"; Double_t textsize_theta = 0.023; Double_t angle_theta = 0; Int_t align_theta =-1; Double_t yoffset_theta =-0.02; Double_t xoffset_theta =-0.01; Bool_t lab_momentum = kTRUE; TString format_momentum = "%5.1f MeV/c"; Double_t textsize_momentum= 0.021; Double_t angle_momentum = 0; Int_t align_momentum =-1; Double_t yoffset_momentum =0.01; Double_t xoffset_momentum =0.01; vector vtheta; vector vmom; Int_t nst = 500; // nsteps bounddary search TObjArray* ar = labels.Tokenize("@"); if(ar) { //----------------------------------------- // extract labels first for(Int_t j=0;jGetEntries();j++){ TString token=((TObjString*)(ar->At(j)))->GetString(); if(token.Contains("theta:")){ token.ReplaceAll("theta:",""); TObjArray* ar1 = token.Tokenize(","); for(Int_t k=0;kGetEntries();k++){ TString token1=((TObjString*)(ar1->At(k)))->GetString(); if(token1.Contains("draw=")){ token1.ReplaceAll("draw=",""); if(token1.CompareTo("yes")==0) lab_theta=kTRUE; else lab_theta=kFALSE; } else if(token1.Contains("format=")){ token1.ReplaceAll("format=",""); format_theta=token1; } else if(token1.Contains("textsize=")){ token1.ReplaceAll("textsize=",""); textsize_theta=token1.Atof(); } else if(token1.Contains("angle=")){ token1.ReplaceAll("angle=",""); angle_theta=token1.Atof(); } else if(token1.Contains("align=")){ token1.ReplaceAll("align=",""); align_theta=token1.Atoi(); } else if(token1.Contains("yoff=")){ token1.ReplaceAll("yoff=",""); yoffset_theta=token1.Atof(); } else if(token1.Contains("xoff=")){ token1.ReplaceAll("xoff=",""); xoffset_theta=token1.Atof(); } else { cout<<"ptyGrid() : unknown token for theta label "<GetEntries();k++){ TString token1=((TObjString*)(ar1->At(k)))->GetString(); if(token1.Contains("draw=")){ token1.ReplaceAll("draw=",""); if(token1.CompareTo("yes")==0) lab_momentum=kTRUE; else lab_momentum=kFALSE; } else if(token1.Contains("format=")){ token1.ReplaceAll("format=",""); format_momentum=token1; } else if(token1.Contains("textsize=")){ token1.ReplaceAll("textsize=",""); textsize_momentum=token1.Atof(); } else if(token1.Contains("angle=")){ token1.ReplaceAll("angle=",""); angle_momentum=token1.Atof(); } else if(token1.Contains("align=")){ token1.ReplaceAll("align=",""); align_momentum=token1.Atoi(); } else if(token1.Contains("yoff=")){ token1.ReplaceAll("yoff=",""); yoffset_momentum=token1.Atof(); } else if(token1.Contains("xoff=")){ token1.ReplaceAll("xoff=",""); xoffset_momentum=token1.Atof(); } else { cout<<"ptyGrid() : unknown token for momentum label "<Delete(); delete ar; ar=NULL; } } ar = setup.Tokenize("@"); if(ar) { for(Int_t j=0;jGetEntries();j++){ TString token=((TObjString*)(ar->At(j)))->GetString(); //----------------------------------------- if(token.Contains("theta:")){ token.ReplaceAll("theta:",""); Bool_t doForLoop = kTRUE; TObjArray* ar1 = token.Tokenize(","); for(Int_t k=0;kGetEntries();k++){ TString token1=((TObjString*)(ar1->At(k)))->GetString(); Int_t n = -1; if(k==0){ n = token1.Atoi(); // positve: forloop, -1: list of values if(n==-1){ doForLoop = kFALSE; continue; } } if(doForLoop) { if(k==0) nTh =token1.Atoi(); if(k==1) strtTh=token1.Atof(); if(k==2) stepTh=token1.Atof(); } else { vtheta.push_back(token1.Atof()); } } // end token loop if(doForLoop){ for(Int_t i=0;i0?vtheta[0]:-1)<< " end "<0?vtheta[n-1]:-1)<Delete(); delete ar1; } } // end theta //----------------------------------------- else if(token.Contains("momentum:")){ token.ReplaceAll("momentum:",""); Bool_t doForLoop = kTRUE; TObjArray* ar1 = token.Tokenize(","); for(Int_t k=0;kGetEntries();k++){ TString token1=((TObjString*)(ar1->At(k)))->GetString(); Int_t n = -1; if(k==0){ n = token1.Atoi(); // positve: forloop, -1: list of values if(n==-1){ doForLoop = kFALSE; continue; } } if(doForLoop){ if(k==0) nMom =token1.Atoi(); if(k==1) strtMom=token1.Atof(); if(k==2) stepMom=token1.Atof(); } else { vmom.push_back(token1.Atof()); } } if(doForLoop){ for(Int_t i=0;i0?vmom[0]:-1)<< " end "<0?vmom[n-1]:-1)<Delete(); delete ar1; } } // end theta else { cout<<"ptyGrid() : unknow token "<Delete(); delete ar; if((opt.CompareTo("draw")==0||opt.CompareTo("drawcopy")==0) && lab_theta) { gPad->Modified(); gPad->Update(); gSystem->ProcessEvents(); } //----------------------------------------- // Now draw for(UInt_t i=0;iDraw("Lsame"); if(opt.CompareTo("drawcopy")==0)f->DrawCopy("Lsame"); if((opt.CompareTo("draw")==0||opt.CompareTo("drawcopy")==0) && lab_theta) { Double_t x1 = gPad->GetUxmin(); Double_t x2 = gPad->GetUxmax(); Double_t y1 = gPad->GetUymin(); Double_t y2 = gPad->GetUymax(); Double_t xmax = -1000; Double_t ymax = -1000; Bool_t Yout = kFALSE; Double_t yoff = (y2-y1)*yoffset_theta; Double_t xoff = (x2-x1)*xoffset_theta; for(Int_t j=0;jEval(x); if(valy1){ xmax=x; ymax=val;} if(x>midRap&&val>y2) Yout=kTRUE; } if(Yout) ymax = y2 +yoff; else xmax = xmax+xoff; if(align_theta<0){ if(!Yout) latex.SetTextAlign(30+1);// right 30, left 10 bottom 1 top 3 hcent 20 vcent 2 else latex.SetTextAlign(10+3);// right 30, left 10 bottom 1 top 3 hcent 20 vcent 2 } else { latex.SetTextAlign(align_theta); } latex.SetTextSize(textsize_theta); TString var = f->GetTitle(); Double_t mom = var.Atof(); latex.DrawLatex(xmax,ymax,Form(format_theta,mom)); } } //----------------------------------------- for(UInt_t i=0;iDraw("Lsame"); if(opt.CompareTo("drawcopy")==0)f->DrawCopy("Lsame"); if((opt.CompareTo("draw")==0||opt.CompareTo("drawcopy")==0) && lab_momentum) { Double_t x1 = gPad->GetUxmin(); Double_t x2 = gPad->GetUxmax(); Double_t y1 = gPad->GetUymin(); Double_t y2 = gPad->GetUymax(); Double_t xmax = -1000; Double_t ymax = -1000; Bool_t Yout = kTRUE; Double_t yoff = (y2-y1)*yoffset_momentum; Double_t xoff = (x2-x1)*xoffset_momentum; if(midRap!=0) // draw left side vertical { for(Int_t j=0;jEval(x); if(valy1&&xmax==-1000){ xmax=x; ymax=val; break;} if(valEval(xmax); if(ymax>0.95*y2) continue; } latex.SetTextSize(textsize_momentum); if(align_momentum<0){ if(midRap==0){ latex.SetTextAlign(10+1);// // right 30, left 10 bottom 1 top 3 hcent 20 vcent 2 latex.SetTextAngle(angle_momentum); } else { latex.SetTextAlign(20+1);// // right 30, left 10 bottom 1 top 3 hcent 20 vcent 2 //latex.SetTextAngle(90+angle_momentum); latex.SetTextAngle(angle_momentum); } } else { if(midRap==0){ latex.SetTextAlign(align_momentum);// // right 30, left 10 bottom 1 top 3 hcent 20 vcent 2 latex.SetTextAngle(angle_momentum); } else { latex.SetTextAlign(align_momentum);// // right 30, left 10 bottom 1 top 3 hcent 20 vcent 2 //latex.SetTextAngle(90+angle_momentum); latex.SetTextAngle(angle_momentum); } } TString var = f->GetTitle(); Double_t mom = var.Atof(); latex.DrawText(xmax,ymax,Form(format_momentum,mom)); } } //----------------------------------------- } return grid; } void HParticleTool::drawPtyGrid(vector& grid,TString opt) { // draws an already existing pt-y grid (see HParticleTool::ptyGrif()) // opt : draw drawcopy opt.ToLower(); for(UInt_t i=0;iDraw("Lsame"); if(opt.CompareTo("drawcopy")==0)grid[i]->DrawCopy("Lsame"); } } //-------------------------------------------------- //------------------------------------------------------------------------------- // Functions from GeomFunct.h (A.Schmah) void HParticleTool::calcSegVector(Double_t z, Double_t rho, Double_t phi, Double_t theta, HGeomVector &base, HGeomVector &dir) { base.setXYZ(rho*cos(phi+TMath::PiOver2()), rho*sin(phi+TMath::PiOver2()), z); dir.setXYZ(sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta)); } Double_t HParticleTool::calcRMS(const Double_t* valArr, Double_t Mean,Int_t valNum) { // Calculates RMS of valNum numbers in valArr using the Mean of the values Double_t RMS = 0.0; Double_t sum = 0.0; for(Int_t i = 0; i < valNum; i++) { sum += pow(valArr[i]-Mean,2.0); } RMS = sqrt(sum/valNum); return RMS; } Double_t HParticleTool::calcDeterminant(HGeomVector& v1, HGeomVector& v2, HGeomVector& v3) { // calculating the Determinant of a 3 x 3 Matrix // with the column vectors [v1, v2, v3] // using the RULE of SARRUS // // | v1(0) v2(0) v3(0) | | v1(0) v2(0) v3(0)|v1(0) v2(0) . // | | | \\ \\ X | / / . // | | | \\ \\ / \\ | / / . // | | | \\ \\ / \\ |/ / . // | | | \\ X \\/ / . // | | | \\ / \\ /\\ / . // | | | \\ / \\ / |\\ / . // | v1(1) v2(1) v3(1) | = | v1(1) v2(1) v3(1)|v1(1) v2(1) . // | | | / \\ /\\ | /\\ . // | | | / \\ / \\ |/ \\ . // | | | / \\/ \\/ \\ . // | | | / /\\ /\\ \\ . // | | | / / \\ / |\\ \\ . // | | | / / \\/ | \\ \\ . // | v1(2) v2(2) v3(2) | | v1(2) v2(2) v3(2)| v1(2) v2(2) . // / / / \\ \\ \\ . // // - - - + + + . return ( v1(0) * v2(1) * v3(2) + v2(0) * v3(1) * v1(2) + v3(0) * v1(1) * v2(2) - v3(0) * v2(1) * v1(2) - v1(0) * v3(1) * v2(2) - v2(0) * v1(1) * v3(2)); } Double_t HParticleTool::calculateMinimumDistanceStraightToPoint(HGeomVector &base, HGeomVector &dir,HGeomVector &point) { //calculates the minimum distance of a point to a straight //given as parametric straight x = base + n * dir if (!(dir.length()>0)) { return -1000000.; } HGeomVector diff = base-point; HGeomVector cross = dir.vectorProduct(diff); return cross.length()/dir.length(); } HGeomVector HParticleTool::calculatePointOfClosestApproachStraightToPoint(HGeomVector &base, HGeomVector &dir,HGeomVector &point) { //calculates the point of minimum distance of a point to a straight //given as parametric straight x = base + n * dir HGeomVector ndir = dir; ndir /= dir.length(); HGeomVector diff = point-base; Double_t l = ndir.scalarProduct(diff); HGeomVector dcap = ndir; dcap *=l; dcap +=base; return dcap; } Double_t HParticleTool::calculateMinimumDistance(HGeomVector &base1, HGeomVector &dir1, HGeomVector &base2, HGeomVector &dir2) { // calculates the minimum distance of two tracks given // as parametric straights x = base + n * dir HGeomVector cross = dir1.vectorProduct(dir2); HGeomVector ab = base1 - base2; if ( !( fabs(cross.length())>0.)) // dir1 || dir2 { return HParticleTool::calculateMinimumDistanceStraightToPoint(base1, dir1, base2); } return fabs(ab.scalarProduct(cross)/cross.length()); } HGeomVector HParticleTool::calculatePointOfClosestApproach(HGeomVector &base1, HGeomVector &dir1, HGeomVector &base2, HGeomVector &dir2) { // calculating point of closest approach // // from the equations of the straight lines of g and h // g: x1 = base1 + l * dir1 // h: x2 = base2 + m * dir2 // // you can construct the following planes: // // E1: e1 = base1 + a * dir1 + b * (dir1 x dir2) // E2: e2 = base2 + s * dir2 + t * (dir1 x dir2) // // now the intersection point of E1 with g2 = {P1} // and the intersection point of E2 with g1 = {P2} // // form the base points of the perpendicular to both straight lines. // // The point of closest approach is the middle point between P1 and P2: // // vertex = (p2 - p1)/2 // // E1 ^ g2: // // e1 = x2 // --> base1 + a * dir1 + b * (dir1 x dir2) = base2 + m * dir2 // --> base1 - base2 = m * dir2 - a * dir1 - b * (dir1 x dir2) // (m) // --> [ dir2, -dir1, -(dir1 x dir2)] (a) = base1 - base2 // (b) // // using CRAMER's RULE you can find the solution for m (a,b, not used) // // using the rules for converting determinants: // // D12 = det [dir2, -dir1, -(dir1 x dir2)] // = det [dir2, dir1, (dir1 x dir2)] // // Dm = det [base1 - base2, -dir1, -(dir1 x dir2)] // = det [base1 - base2, dir1, (dir1 x dir2)] // // m = Dm/D12 // // P1: p1 = x2(m) // = base2 + Dm/D12 * dir2 // // E2 ^ g1: // // e2 = x1 // --> base2 + s * dir2 + t * (dir1 x dir2) = base1 + l * dir1 // --> base2 - base1 = l * dir1 - s * dir2 - t * (dir1 x dir2) // (l) // --> [ dir1, -dir2, -(dir1 x dir2)] (s) = base2 - base1 // (t) // // again using CRAMER's RULE you can find the solution for m (a,b, not used) // // using the rules for converting determinants: // // D21 = det [dir1, -dir2, -(dir1 x dir2)] // = det [dir1, dir2, (dir1 x dir2)] // = -det [dir2, dir1, (dir1 x dir2)] // = -D12 // // Dl = det [base2 - base1, -dir2, -(dir1 x dir2)] // = det [base2 - base1, dir1, (dir1 x dir2)] // = -det [base1 - base2, dir1, (dir1 x dir2)] // // l = Dl/D21 // = - Dl/D12 // // P2: p2 = x1(m) // = base1 - Dl/D12 * dir1 // // // vertex = p1 + 1/2 * (p2 - p1) // = 1/2 * (p2 + p1) // = 1/2 *( (base1 + base2) + 1/D12 * ( Dm * dir2 - Dl * dir1) ) // HGeomVector cross = dir1.vectorProduct(dir2); // cross product: dir1 x dir2 // straight lines are either skew or have a cross point HGeomVector diff = base1; diff-=base2; // Difference of two base vectors base1 - base2 Double_t D; D = calcDeterminant(dir2, dir1 ,cross); if (!(fabs(D) > 0.)) { ::Warning(":calculatePointOfClosestApproach","Dirs and cross-product are lin. dependent: returning default Vertex (-20000,-20000,-20000)"); return HGeomVector(-20000.,-20000.,-20000.); } Double_t Dm = calcDeterminant(diff , dir1, cross); Double_t Dl = -calcDeterminant(diff , dir2, cross); HGeomVector vertex; HGeomVector dm; HGeomVector dl; dm = dir2; dm *= Dm; dl = dir1; dl *= Dl; vertex = dm - dl; vertex *= ((1.)/D); vertex+=base1; vertex+=base2; vertex*=0.5; return HGeomVector(vertex); } HGeomVector HParticleTool::calculateCrossPoint(HGeomVector &base1, HGeomVector &dir1, HGeomVector &base2, HGeomVector &dir2) { // calculating cross point // taking all three equations into account solving the overdetermined set of lin. equations // of // base1 + l * dir2 = base1 + m * dir2 // // set of lin. equations: // // base1(0) + l * dir1(0) = base2(0) + m * dir2(0) // base1(1) + l * dir1(1) = base2(1) + m * dir2(1) // base1(2) + l * dir1(2) = base2(2) + m * dir2(2) this line is ignored // // written in matrix form // // l \\ . // M * | | = base2 - base1 // \\ m / // // M is a 3x2 matrix // // to solve multiply the equation by the transposed Matrix of M from the left: M // // T / l \\ . // M * M * | | = M * (base2 - base1) // \\ -m / // MIND THE '-' of m // // / dir1(0) dir2(0) \\ . // | | T / dir1(0) dir1(1) dir1(2) \\ . // M = | dir1(1) dir2(1) |, M = | | // | | \\ dir2(0) dir2(1) dir2(2) / . // \\ dir1(2) dir2(2) / // // T / (dir1(0)*dir1(0) + dir1(1)*dir1(1) + dir1(2)*dir1(2)) (dir1(0)*dir2(0) + dir1(1)*dir2(1) + dir1(2)*dir2(2)) \\ . // // M * M = | | // // \\ (dir1(0)*dir2(0) + dir1(1)*dir2(1) + dir1(2)*dir2(2)) (dir2(0)*dir2(0) + dir2(1)*dir2(1) + dir2(2)*dir2(2)) / // // // T / d1d1 d1d2 \\ . // M * M = | | // \\ d1d2 d2d2 / // // diff = base2 - base1 // // T / (dir1(0)*diff(0) + dir1(1)*diff(1) + dir1(2)*diff(2)) \\ . // M * diff = | | // \\ (dir2(0)*diff(0) + dir2(1)*diff(1) + dir2(2)*diff(2)) / // // T / d1diff \\ . // M * diff = | | // \\ d2diff / // // now the new Matrix set is to be solved by CRAMER'S Rule: // // / d1d1 d1d2 \\ / l \\ / d1diff \\ . // | | * | | = | | // \\ d1d2 d2d2 / \\ -m / \\ d2diff / // // | d1d1 d1d2 | // D = | | = d1d1*d2d2 - d1d2*d1d2; // | d1d2 d2d2 | // // | d1diff d1d2 | // Dl= | | = d1diff*d2d2 - d1d2*d2diff; // | d2diff d2d2 | // // l = Dl/D = l_cross // // vertex = base1 + l_cross * dir1 // Double_t d1d1 = dir1(0)*dir1(0) + dir1(1)*dir1(1) + dir1(2)*dir1(2); Double_t d2d2 = dir2(0)*dir2(0) + dir2(1)*dir2(1) + dir2(2)*dir2(2); Double_t d1d2 = dir1(0)*dir2(0) + dir1(1)*dir2(1) + dir1(2)*dir2(2); Double_t D = d1d1*d2d2 - (d1d2*d1d2); if (!(fabs(D) > 0.)) { ::Warning("calculateCrossPoint","Error while calculating cross point ... eqns are lin. dependent:returning default Vertex (-20000,-20000,-20000)"); return HGeomVector(-20000.,-20000.,-20000.); } Double_t d1diff = dir1(0)*(base2(0)-base1(0))+dir1(1)*(base2(1)-base1(1))+dir1(2)*(base2(2)-base1(2)); Double_t d2diff = dir2(0)*(base2(0)-base1(0))+dir2(1)*(base2(1)-base1(1))+dir2(2)*(base2(2)-base1(2)); Double_t Dlambda = d1diff*d2d2-d1d2*d2diff; Double_t lambda = Dlambda/D; HGeomVector vertex; vertex += dir1; vertex *= lambda; vertex += base1; return HGeomVector(vertex); } HGeomVector HParticleTool::calcVertexAnalytical(HGeomVector &base1, HGeomVector &dir1, HGeomVector &base2, HGeomVector &dir2) { // Calculates the Vertex of two straight lines defined by the vectors base and dir // // g: x1 = base1 + l * dir1 // h: x2 = base2 + m * dir2 , where l,m are real numbers // h // 1. are g and h // parallel / identical, i.e. are dir1 and dir2 linear dependent? // // /- // | // | = 0 linear dependent, no unique solution, returning dummy // => cross product : dir1 x dir2 = -| // | != 0 linear independent // | // \\- // // 2. are g and h // skew or do they have a crossing point, i.e are dir1, dir2 and (base1 - base2) linear dependent ? // // /- // | // | = 0 linear dependent // | g and h are intersecting // | calculating vertex as point of intersection // | // => determinant: det[ dir1, dir2, base1-base2]= -| // | != 0 linear independent // | g and h are skew // | calulating vertex as point of closest approach // | // \\- // // 3. // (a) calculating intersection point // (b) calculating point of closest approach // 1. exists a unique solution ? if ((dir1.vectorProduct(dir2)).length()> 0.) // dir1 and dir2 linear independent { // straight lines are either skew or have a cross point HGeomVector diff = base1; diff-=base2; // Difference of two base vectors base1 - base2 // 2. skew or intersecting ? if (fabs(HParticleTool::calcDeterminant(dir2, dir1 ,diff))>0.) { // 3. (b) skew return HGeomVector(HParticleTool::calculatePointOfClosestApproach(base1, dir1, base2, dir2)); } else { // 3. (a) intersection return HGeomVector(HParticleTool::calculateCrossPoint(base1 ,dir1, base2 ,dir2)); } } else { // dir1 and dir2 linear dependent -> g1 and g2 identical or parallel return HGeomVector(-10000000.,-10000000.,-10000000.); } return HGeomVector(-10000000.,-10000000.,-10000000.); } // end functions from GeomFunct.h //------------------------------------------------------------------------------- Int_t HParticleTool::findFirstHitInTof(Int_t trackID,Int_t modeTrack) { //------------------------------------------------- // find the first track ID entering the TOF // Used to suppress the secondaries created in the // TOF itself. // 0 = realistic (secondaries included) // 1 primary particle is stored // 2 (default) the first track number entering the tof in SAME SECTOR is stored // 3 as 2 but condition on SAME SECTOR && MOD // 4 as 2 but SAME SECTOR && MOD && ROD Int_t numTrack = trackID; if(numTrack <= 0) return numTrack; // nothing to do for negative track numbers //------------------------------------------ // getting categories HLinearCategory* fGeantKineCat = (HLinearCategory*)HCategoryManager::getCategory(catGeantKine, kFALSE); if(!fGeantKineCat){ return trackID; } HLinearCategory* fGeantCat = (HLinearCategory*) HCategoryManager::getCategory(catTofGeantRaw,kFALSE); if(!fGeantCat){ return trackID; } //------------------------------------------ HGeantTof *poldTof; Int_t first = 0; HGeantKine* kine = 0; kine = HCategoryManager::getObject(kine,fGeantKineCat,numTrack - 1); if(kine){ first=kine->getFirstTofHit(); if(first != -1){ poldTof = (HGeantTof*)fGeantCat->getObject(first); } else { ::Error("findFirstHitInTof()","No first tof hit!"); return numTrack; } } else { ::Error("findFirstHitInTof()","Received Zero pointer for kine index = %i !",numTrack); return numTrack; } if(numTrack != poldTof->getTrack()){ ::Error("findFirstHitInTof()","First tof hit not same trackID!"); return numTrack; } //-------------------------------------------------------- // return the track number for // the selected option storeFirstTrack //-------------------------------------- // case 0 if(modeTrack == 0) return numTrack; //-------------------------------------- // case 1 if(modeTrack == 1) { // return track number of primary particle // of the given track kine = (HGeantKine*)fGeantKineCat->getObject(numTrack-1); kine = HGeantKine::getPrimary(numTrack,fGeantKineCat); return kine->getTrack(); } //-------------------------------------- // case 2 and 3 kine = (HGeantKine*)fGeantKineCat->getObject(numTrack - 1); if(kine->getParentTrack() == 0){return numTrack;} // nothing to do Int_t s,m,c; s = poldTof->getSector(); // GEANT 15-22 (TOF), 23-26 (TOFINO) // numbering reverse ! m = 21-poldTof->getModule(); c = 7 -poldTof->getCell(); first = 0; Int_t tempTrack = numTrack; while((kine = kine->getParent(tempTrack,fGeantKineCat)) != 0) { first = kine->getFirstTofHit(); if(first != -1) { // track is still in TOF // now we have to check if it is in TOF or TOFINO HGeantTof* gtof = (HGeantTof*)fGeantCat->getObject(first); Int_t s1,m1,c1; s1 = m1 = c1 = 0; m1 = 21 - gtof->getModule(); if(m1 >= 0) { // inside TOF s1 = gtof->getSector(); c1 = 7-gtof->getCell(); if(modeTrack == 2 && s == s1) { // case 2 :: check only sector tempTrack = kine->getTrack(); } if(modeTrack == 3 && s == s1 && m == m1) { // case 3 :: check only sector,module tempTrack = kine->getTrack(); } else if(modeTrack == 4 && s == s1 && m == m1 && c == c1) { // case 4 :: check for same rod tempTrack = kine->getTrack(); } else { // track has no TOF hit any longer // which fulfills the condition break; } } else { // track is in TOFINO break; } } else { // track has no TOF hit any longer, // so the previous object was the one we // searched for break; } } return tempTrack; } Int_t HParticleTool::findFirstHitShowerInTofino(Int_t trackID, Int_t modeTrack) { //------------------------------------------------- // Used to suppress the secondaries created in the // SHOWER itself. // 0 = realistic (secondaries included) // 1 primary particle is stored // 2 the first track number entering the SHOWER in SAME SECTOR is stored // 3 the first track number entering the TOFINO in SAME SECTOR is stored // or the primary track if no TOFINO was found // 4 (default) the first track number entering the TOFINO in SAME SECTOR is stored // or the first track in SHOWER if no TOFINO was found Int_t numTrack = trackID; if(numTrack <= 0) return numTrack; // nothing to do for negative track numbers //-------------------------------------- // case 0 : in = out if(modeTrack == 0) { return numTrack; } HGeantShower *poldShower; Int_t first = 0; Int_t parent= 0; //------------------------------------------ // getting categories HLinearCategory* fGeantKineCat = (HLinearCategory*)HCategoryManager::getCategory(catGeantKine, kFALSE); if(!fGeantKineCat){ return trackID; } HLinearCategory* fGeantShowerCat = (HLinearCategory*) HCategoryManager::getCategory(catShowerGeantRaw,kFALSE); if(!fGeantShowerCat){ return trackID; } HLinearCategory* fGeantTofCat = (HLinearCategory*) HCategoryManager::getCategory(catTofGeantRaw,kFALSE); if(!fGeantTofCat){ return trackID; } //------------------------------------------ HGeantKine* kine = (HGeantKine*)fGeantKineCat->getObject(numTrack - 1); if(kine){ parent = kine->getParentTrack(); if( parent == 0) { return numTrack; } // nothing todo first = kine->getFirstShowerHit(); if(first != -1){ poldShower = (HGeantShower*)fGeantShowerCat->getObject(first); } else { ::Error("findFirstHitShowerInTofino()","No first shower hit!"); return numTrack; } } else { ::Error("findFirstHitShowerInTofino()","Received Zero pointer for kine!"); return numTrack; } if(numTrack != poldShower->getTrack()){ ::Error("findFirstHitShowerInTofino()","First shower hit not same trackID!"); return numTrack; } //-------------------------------------------------------- // return the track number for // the selected option modeTrack Int_t s = poldShower->getSector(); //-------------------------------------- // case 1 : in -> primary if(modeTrack == 1) { // return track number of primary particle // of the given track kine = HGeantKine::getPrimary(numTrack,fGeantKineCat); return kine->getTrack(); } //-------------------------------------- // case 2 and 3 : entering SHOWER/TOF Int_t tempTrack = numTrack; first = 0; //-------------------------------------- // case 2 : entering SHOWER if(modeTrack == 2) { while((kine = kine->getParent(tempTrack,fGeantKineCat)) != 0) { first = kine->getFirstShowerHit(); if(first != -1) { // track is still in SHOWER HGeantShower* gshower = (HGeantShower*)fGeantShowerCat->getObject(first); Int_t s1 = gshower->getSector(); if(s == s1) { // check only sector tempTrack = kine->getTrack(); } else { // track has no SHOWER hit any longer // which fulfills the condition break; } } else { // track has no SHOWER hit any longer, // so the previous object was the one we // searched for break; } } return tempTrack; } //-------------------------------------- //-------------------------------------- // case 3 : entering TOFINO if(modeTrack >= 3) { Bool_t foundTof = kFALSE; Int_t tempTrack2 = tempTrack; do { first = kine->getFirstTofHit(); if(first != -1) { // we are in TOF HGeantTof* gtof = (HGeantTof*)fGeantTofCat->getObject(first); Int_t s1 = gtof->getSector(); Int_t m = gtof->getModule(); if(s == s1 && m > 21 ) { // check only sector + TOFINO foundTof = kTRUE; tempTrack = tempTrack2; } } tempTrack2 = kine->getParentTrack(); } while( tempTrack2 > 0 && (kine = (HGeantKine*)fGeantKineCat->getObject(tempTrack2 - 1)) != 0); if(foundTof) { tempTrack += 10000; } else { kine = (HGeantKine*)fGeantKineCat->getObject(numTrack - 1); if( modeTrack == 3 ){ // store primaries if no TOFino was found kine = HGeantKine::getPrimary(numTrack,fGeantKineCat); tempTrack = kine->getTrack() + 100000000; } else if (modeTrack == 4){ //-------------------------------------- // recover first particle entering SHOWER while((kine = kine->getParent(tempTrack,fGeantKineCat)) != 0) { first = kine->getFirstShowerHit(); if(first != -1) { // track is still in SHOWER HGeantShower* gshower = (HGeantShower*)fGeantShowerCat->getObject(first); Int_t s1 = gshower->getSector(); if(s == s1) { // check only sector tempTrack = kine->getTrack(); } else { // track has no SHOWER hit any longer // which fulfills the condition break; } } else { // track has no SHOWER hit any longer, // so the previous object was the one we // searched for break; } } tempTrack += 100000000; //-------------------------------------- } } } //-------------------------------------- return tempTrack; } Int_t HParticleTool::findFirstHitShowerInRpc(Int_t trackID,Int_t modeTrack) { //------------------------------------------------- // Used to suppress the secondaries created in the // SHOWER itself. // 0 = realistic (secondaries included) // 1 primary particle is stored // 2 the first track number entering the SHOWER in SAME SECTOR is stored // 3 the first track number entering the RPC in SAME SECTOR is stored // or the primary track if no RPC was found // 4 (default) the first track number entering the RPC in SAME SECTOR is stored // or the first track in SHOWER if no RPC was found Int_t numTrack = trackID; if(numTrack <= 0) return numTrack; // nothing to do for negative track numbers //-------------------------------------- // case 0 : in = out if(modeTrack == 0) { return numTrack; } HGeantShower *poldShower; Int_t first = 0; Int_t parent= 0; //------------------------------------------ // getting categories HLinearCategory* fGeantKineCat = (HLinearCategory*)HCategoryManager::getCategory(catGeantKine, kFALSE); if(!fGeantKineCat){ return trackID; } HLinearCategory* fGeantShowerCat = (HLinearCategory*) HCategoryManager::getCategory(catShowerGeantRaw,kFALSE); if(!fGeantShowerCat){ return trackID; } HLinearCategory* fGeantRpcCat = (HLinearCategory*) HCategoryManager::getCategory(catRpcGeantRaw,kFALSE); if(!fGeantRpcCat){ return trackID; } //------------------------------------------ HGeantKine* kine = (HGeantKine*)fGeantKineCat->getObject(numTrack - 1); if(kine){ parent = kine->getParentTrack(); if( parent == 0) { return numTrack; } // nothing todo first = kine->getFirstShowerHit(); if(first != -1){ poldShower = (HGeantShower*)fGeantShowerCat->getObject(first); } else { ::Error("findFirstHitShowerInRpc()","No first shower hit!"); return numTrack; } } else { ::Error("findFirstHitShowerInRpc()","Received Zero pointer for kine!"); return numTrack; } if(numTrack != poldShower->getTrack()){ ::Error("findFirstHitShowerInRpc()","First shower hit not same trackID!"); return numTrack; } //-------------------------------------------------------- // return the track number for // the selected option modeTrack Int_t s = poldShower->getSector(); //-------------------------------------- // case 1 : in -> primary if(modeTrack == 1) { // return track number of primary particle // of the given track kine = HGeantKine::getPrimary(numTrack,fGeantKineCat); return kine->getTrack(); } //-------------------------------------- // case 2 and 3 : entering SHOWER/RPC Int_t tempTrack = numTrack; first = 0; //-------------------------------------- // case 2 : entering SHOWER if(modeTrack == 2) { while((kine = kine->getParent(tempTrack,fGeantKineCat)) != 0) { first = kine->getFirstShowerHit(); if(first != -1) { // track is still in SHOWER HGeantShower* gshower = (HGeantShower*)fGeantShowerCat->getObject(first); Int_t s1 = gshower->getSector(); if(s == s1) { // check only sector tempTrack = kine->getTrack(); } else { // track has no SHOWER hit any longer // which fulfills the condition break; } } else { // track has no SHOWER hit any longer, // so the previous object was the one we // searched for break; } } return tempTrack; } //-------------------------------------- //-------------------------------------- // case 3 : entering RPC if(modeTrack >= 3) { Bool_t foundRPC = kFALSE; Int_t tempTrack2 = tempTrack; do { first = kine->getFirstRpcHit(); if(first != -1) { // we are in TOF HGeantRpc* grpc = (HGeantRpc*)fGeantRpcCat->getObject(first); Int_t s1 = grpc->getSector(); if(s == s1) { // check only sector + RPC foundRPC = kTRUE; tempTrack = tempTrack2; } } tempTrack2 = kine->getParentTrack(); } while( tempTrack2 > 0 && (kine = (HGeantKine*)fGeantKineCat->getObject(tempTrack2 - 1)) != 0); if(foundRPC) { tempTrack += 10000; } else { kine = (HGeantKine*)fGeantKineCat->getObject(numTrack - 1); if( modeTrack == 3 ){ // store primaries if no RPC was found kine = HGeantKine::getPrimary(numTrack,fGeantKineCat); tempTrack = kine->getTrack() + 100000000; } else if (modeTrack == 4){ //-------------------------------------- // recover first particle entering SHOWER while((kine = kine->getParent(tempTrack,fGeantKineCat)) != 0) { first = kine->getFirstShowerHit(); if(first != -1) { // track is still in SHOWER HGeantShower* gshower = (HGeantShower*)fGeantShowerCat->getObject(first); Int_t s1 = gshower->getSector(); if(s == s1) { // check only sector tempTrack = kine->getTrack(); } else { // track has no SHOWER hit any longer // which fulfills the condition break; } } else { // track has no SHOWER hit any longer, // so the previous object was the one we // searched for break; } } tempTrack += 100000000; //-------------------------------------- } } } //-------------------------------------- return tempTrack; } Float_t HParticleTool::getInterpolatedValue(TH1* h, Float_t xVal, Bool_t warn) { // retrieve content of 1-dim Hist corresponding to x val. // The values will be linear interpolated. If warn == kTRUE // warning will be emitted if the xval is out of the range // of the hist. In case the value is out of range the value // of first or last bin will be returned respectively. if(!h) return 0; TAxis *pA = h->GetXaxis(); Int_t binX = pA->FindBin(xVal); // We need to start interpolation from the next-lower bin boundary if(binX > 1 && binX <= pA->GetNbins()) { // only here it make sense to interpolate Int_t binXLower = binX - 1; Float_t startVal = h->GetBinContent(binXLower); Float_t endVal = h->GetBinContent(binX); Float_t startx = pA->GetBinUpEdge(binXLower); Float_t slope = (endVal - startVal)/pA->GetBinWidth(binX); Float_t calc = startVal + slope * (xVal - startx); return calc; } else { // noting to interpolate.... if(binX == 1) { return h->GetBinContent(1); } else if(binX < 1) { // lower boundary violation ...return lowest bin if(warn) ::Warning("HParticleTool::getInterpolatedValue()","Non-intercepted bin-out-of-range-situation (bin too low)"); return h->GetBinContent(1); } else { // upper boundary violation ...return highest bin if(warn) ::Warning("HParticleTool::getInterpolatedValue()","Non-intercepted bin-out-of-range-situation (bin too high)"); return h->GetBinContent(pA->GetNbins()); } } } Stat_t HParticleTool::getValue(TH1* h,Float_t xVal, Float_t yVal, Float_t zVal) { // retrieve value from Histogram (1,2 and 3 dim) corresponding // to x val, y val and z val. If the values are out of range of // the histogram axis the lowest/highest bin of the axis will be // used. No interpolation performed. if(!h) return 0; TAxis *pA = h->GetXaxis(); Int_t binX, binY, binZ; binY = binZ = 0; pA = h->GetXaxis(); binX = pA->FindBin(xVal); if(binX <= 0) { binX = 1; } else { if(binX > pA->GetNbins()) binX = pA->GetNbins(); } pA = h ->GetYaxis(); binY = pA->FindBin(yVal); if(binY <= 0) { binY = 1; } else { if(binY > pA->GetNbins()) binY = pA->GetNbins(); } pA = h ->GetZaxis(); binZ = pA->FindBin(zVal); if(binZ <= 0) { binZ = 1; } else { if(binZ > pA->GetNbins()) binZ = pA->GetNbins(); } return h->GetBinContent(h->GetBin(binX, binY, binZ)); } HRichHit* HParticleTool::getRichHit(Int_t richind) { // return HRichHit pointer from the object index // in the category. In case of no success returns NULL HRichHit* richhit = 0; richhit = HCategoryManager::getObject(richhit,catRichHit,richind,kTRUE); return richhit; } HTofHit* HParticleTool::getTofHit(Int_t tofind) { // return HTofHit pointer from the object index // in the category. In case of no success returns NULL HTofHit* tofhit=0; tofhit = HCategoryManager::getObject(tofhit,catTofHit,tofind,kTRUE); return tofhit; } HTofCluster* HParticleTool::getTofCluster(Int_t tofind) { // return HTofCluster pointer from the object index // in the category. In case of no success returns NULL HTofCluster* tofclst=0; tofclst = HCategoryManager::getObject(tofclst,catTofCluster,tofind,kTRUE); return tofclst; } HRpcCluster* HParticleTool::getRpcCluster(Int_t rpcind) { // return HRpcCluster pointer from the object index // in the category. In case of no success returns NULL HRpcCluster* rpcclst=0; rpcclst = HCategoryManager::getObject(rpcclst,catRpcCluster,rpcind,kTRUE); return rpcclst; } HShowerHit* HParticleTool::getShowerHit(Int_t showerind) { // return HShowerHit pointer from the object index // in the category. In case of no success returns NULL HShowerHit* showerhit=0; showerhit = HCategoryManager::getObject(showerhit,catShowerHit,showerind,kTRUE); return showerhit; } HMetaMatch2* HParticleTool::getMetaMatch(Int_t metaind) { // return HMetaMatch2 pointer from the object index // in the category. In case of no success returns NULL HMetaMatch2* meta = 0; meta = HCategoryManager::getObject(meta,catMetaMatch,metaind,kTRUE); return meta; } HMdcTrkCand* HParticleTool::getMdcTrkCand(Int_t metaind) { // return HMdcTrkCand pointer from the object index of MetaMatch2 // in the category (indMeta->HMetaMatch2->HMdcTrkCand). In case of // no success returns NULL HMetaMatch2* meta = 0; HMdcTrkCand* trk = 0; meta = HCategoryManager::getObject(meta,catMetaMatch,metaind); if(meta){ Int_t mdctrkInd = meta->getTrkCandInd(); if(mdctrkInd >= 0){ trk = HCategoryManager::getObject(trk,catMdcTrkCand,mdctrkInd,kTRUE); } } return trk; } HMdcSeg* HParticleTool::getMdcSeg(Int_t segind) { // return HMdcSeg pointer from the object index // in the category. In case of no success returns NULL HMdcSeg* seg = 0; seg = HCategoryManager::getObject(seg,catMdcSeg,segind,kTRUE); return seg; } HMdcHit* HParticleTool::getMdcHit(Int_t segind,Int_t nhit ) { // return HMdcHit pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcHit). nhit (0 or 1) // is the number of the hit inside the segment which should // be retrived. In case of no success returns NULL HMdcSeg* seg = 0; HMdcHit* hit = 0; seg = HCategoryManager::getObject(seg,catMdcSeg,segind,kTRUE); if(seg){ Int_t hitind = seg->getHitInd(nhit); if(hitind != -1) { hit = HCategoryManager::getObject(hit,catMdcHit,hitind,kTRUE); } } return hit; } HMdcClusInf* HParticleTool::getMdcClusInf(Int_t segind,Int_t nhit) { // return HMdcClusInf pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcClusInf). nhit (0 or 1) // is the number of the hit inside the segment which should // be retrived. In case of no success returns NULL HMdcSeg* seg = 0; HMdcClusInf* clusinf = 0; seg = HCategoryManager::getObject(seg,catMdcSeg,segind,kTRUE); if(seg){ Int_t hitind = seg->getHitInd(nhit); if(hitind != -1) { clusinf = HCategoryManager::getObject(clusinf,catMdcClusInf,hitind,kTRUE); } } return clusinf; } HMdcClusFit* HParticleTool::getMdcClusFit(Int_t segind) { // return HMdcClusFit pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcClsuInf->HMdcClusFit). // In case of no success returns NULL HMdcSeg* seg = 0; HMdcClusInf* clusinf = 0; HMdcClusFit* clusfit = 0; seg = HCategoryManager::getObject(seg,catMdcSeg,segind,kTRUE); if(seg){ Int_t hitind = seg->getHitInd(0); if(hitind != -1) { clusinf = HCategoryManager::getObject(clusinf,catMdcClusInf,hitind,kTRUE); } else { hitind = seg->getHitInd(1); clusinf = HCategoryManager::getObject(clusinf,catMdcClusInf,hitind,kTRUE); } if(clusinf){ Int_t clusfitind = clusinf ->getClusFitIndex(); if(clusfitind !=-1){ clusfit = HCategoryManager::getObject(clusfit,catMdcClusFit,clusfitind,kTRUE); } } } return clusfit; } HMdcClus* HParticleTool::getMdcClus(Int_t segind) { // return HMdcClus pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcClus). // In case of no success returns NULL HMdcSeg* seg = 0; HMdcClus* clus = 0; seg = HCategoryManager::getObject(seg,catMdcSeg,segind,kTRUE); if(seg){ Int_t clusind = seg->getClusIndex(); if(clusind != -1) { clus = HCategoryManager::getObject(clus,catMdcClus,clusind,kTRUE); } } return clus; } TObjArray* HParticleTool::getMdcWireFitSeg(Int_t segind) { // return TObjArray pointer to an array of HMdcWireFit // pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcClusInf->HMdcClusFit->HMdcWireFit). // The user has to delete the TObjArray object // (not the content) by himself. TObjArray* ar = new TObjArray(); HMdcClusFit* clusfit = HParticleTool::getMdcClusFit(segind); if(clusfit){ Int_t first = clusfit->getFirstWireFitInd(); Int_t last = clusfit->getLastWireFitInd(); HCategory *wirefitcat = HCategoryManager::getCategory(catMdcWireFit,2);// silent if(wirefitcat){ for(Int_t i = first; i<= last; i++){ HMdcWireFit* wf = 0; wf = HCategoryManager::getObject(wf,wirefitcat,i,kTRUE); if(wf) ar->Add(wf); } } } return ar; } Int_t HParticleTool::getMdcWireFitSeg(Int_t segind, vector& v, Bool_t clear) { // fill vector of HMdcWireFit // pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcClusInf->HMdcClusFit->HMdcWireFit). // The user has to delete the TObjArray object // (not the content) by himself. Returns the number of // collected objects .the vector will be cleared // if clear=kTRUE if(clear) v.clear(); HMdcClusFit* clusfit = HParticleTool::getMdcClusFit(segind); if(clusfit){ Int_t first = clusfit->getFirstWireFitInd(); Int_t last = clusfit->getLastWireFitInd(); HCategory *wirefitcat = HCategoryManager::getCategory(catMdcWireFit,2);// silent if(wirefitcat){ for(Int_t i = first; i<= last; i++){ HMdcWireFit* wf = 0; wf = HCategoryManager::getObject(wf,wirefitcat,i,kTRUE); v.push_back(wf); } } } return v.size(); } TObjArray* HParticleTool::getMdcCal1Seg(Int_t segind) { // return TObjArray pointer to an array of HMdcCal1 // pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcCal1). // user has to delete the TObjArray object // (not the content) by himself. TObjArray* ar = new TObjArray(); HMdcSeg* seg = HParticleTool::getMdcSeg(segind); HCategory *cal1cat = HCategoryManager::getCategory(catMdcCal1,2);// silent if(seg && cal1cat){ HLocation loccal; loccal.set(4,0,0,0,0); // t2-t1 from segment for(Int_t l = 0; l < 12; l ++) // loop over layers in segment { Int_t nCell = seg->getNCells(l); loccal[0] = seg->getSec(); Int_t io = seg->getIOSeg(); for(Int_t ic = 0; ic < nCell; ic ++) { //--------------------------------------------- // find cal1 coordinates // find module if(io == 0){ (l < 6)? loccal[1] = 0 : loccal[1] = 1; } else if(io == 1) { (l < 6)? loccal[1] = 2 : loccal[1] = 3; } // find layer (l < 6)? loccal[2] = l : loccal[2] = l - 6; // find cell number loccal[3] = seg->getCell(l,ic); //--------------------------------------------- HMdcCal1* cal1 = (HMdcCal1*)cal1cat->getObject(loccal); if(cal1) { ar->Add(cal1); } } } } return ar; } Int_t HParticleTool::getMdcCal1Seg(Int_t segind, vector& v, Bool_t clear) { // fill vector of HMdcCal1 // pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcCal1). // user has to delete the TObjArray object // (not the content) by himself. Returns number // of collected objects. the vector will be cleared // if clear=kTRUE if(clear) v. clear(); HMdcSeg* seg = HParticleTool::getMdcSeg(segind); HCategory *cal1cat = HCategoryManager::getCategory(catMdcCal1,2);// silent if(seg && cal1cat){ HLocation loccal; loccal.set(4,0,0,0,0); // t2-t1 from segment for(Int_t l = 0; l < 12; l ++) // loop over layers in segment { Int_t nCell = seg->getNCells(l); loccal[0] = seg->getSec(); Int_t io = seg->getIOSeg(); for(Int_t ic = 0; ic < nCell; ic ++) { //--------------------------------------------- // find cal1 coordinates // find module if(io == 0){ (l < 6)? loccal[1] = 0 : loccal[1] = 1; } else if(io == 1) { (l < 6)? loccal[1] = 2 : loccal[1] = 3; } // find layer (l < 6)? loccal[2] = l : loccal[2] = l - 6; // find cell number loccal[3] = seg->getCell(l,ic); //--------------------------------------------- HMdcCal1* cal1 = (HMdcCal1*)cal1cat->getObject(loccal); if(cal1) { v.push_back(cal1); } } } } return v.size(); } TObjArray* HParticleTool::getMdcCal1Cluster(Int_t segind) { // return TObjArray pointer to an array of HMdcCal1 // pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcClus->HMdcCal1). // user has to delete the TObjArray object // (not the content) by himself. TObjArray* ar = new TObjArray(); HMdcClus* clus = HParticleTool::getMdcClus(segind); HCategory *cal1cat = HCategoryManager::getCategory(catMdcCal1,2);// silent if(clus && cal1cat){ HLocation loccal; loccal.set(4,0,0,0,0); // t2-t1 from segment for(Int_t l = 0; l < 12; l ++) // loop over layers in segment { Int_t nCell = clus->getNCells(l); loccal[0] = clus->getSec(); Int_t io = clus->getIOSeg(); for(Int_t ic = 0; ic < nCell; ic ++) { //--------------------------------------------- // find cal1 coordinates // find module if(io == 0){ (l < 6)? loccal[1] = 0 : loccal[1] = 1; } else if(io == 1) { (l < 6)? loccal[1] = 2 : loccal[1] = 3; } // find layer (l < 6)? loccal[2] = l : loccal[2] = l - 6; // find cell number loccal[3] = clus->getCell(l,ic); //--------------------------------------------- HMdcCal1* cal1 = (HMdcCal1*)cal1cat->getObject(loccal); if(cal1) { ar->Add(cal1); } } } } return ar; } Int_t HParticleTool::getMdcCal1Cluster(Int_t segind, vector& v, Bool_t clear) { // fill vector of HMdcCal1 // pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcClus->HMdcCal1). // user has to delete the TObjArray object // (not the content) by himself. Returns the number // of collected objects. The vector will be cleared // if clear=kTRUE if(clear) v. clear(); HMdcClus* clus = HParticleTool::getMdcClus(segind); HCategory *cal1cat = HCategoryManager::getCategory(catMdcCal1,2);// silent if(clus && cal1cat){ HLocation loccal; loccal.set(4,0,0,0,0); // t2-t1 from segment for(Int_t l = 0; l < 12; l ++) // loop over layers in segment { Int_t nCell = clus->getNCells(l); loccal[0] = clus->getSec(); Int_t io = clus->getIOSeg(); for(Int_t ic = 0; ic < nCell; ic ++) { //--------------------------------------------- // find cal1 coordinates // find module if(io == 0){ (l < 6)? loccal[1] = 0 : loccal[1] = 1; } else if(io == 1) { (l < 6)? loccal[1] = 2 : loccal[1] = 3; } // find layer (l < 6)? loccal[2] = l : loccal[2] = l - 6; // find cell number loccal[3] = clus->getCell(l,ic); //--------------------------------------------- HMdcCal1* cal1 = (HMdcCal1*)cal1cat->getObject(loccal); if(cal1) { v.push_back(cal1); } } } } return v.size(); } Bool_t HParticleTool::getSimTracks(HParticleCandSim* cand, vector&tracksMeta, vector&tracksShower, vector&tracksRich, vector&weightRich, vector&tracksInnerMdc, vector&weightInnerMdc, vector&tracksOuterMdc, vector&weightOuterMdc, Bool_t print ){ tracksMeta.clear(); tracksShower.clear(); tracksRich.clear(); tracksInnerMdc.clear(); tracksOuterMdc.clear(); weightInnerMdc.clear(); weightOuterMdc.clear(); Bool_t ret =kTRUE; if(print)cout<<"SIM tracks----------------------------------------------"<getInnerSegInd()); seg [1]= (HMdcSegSim*) HParticleTool::getMdcSeg(cand->getOuterSegInd()); for(Int_t io=0;io<2;io++){ if(seg[io]==0) continue; for(Int_t tr=0;trgetNTracks(); tr++){ if(io==0){ tracksInnerMdc.push_back((Int_t)seg[io]->getTrack(tr)); weightInnerMdc.push_back((Int_t)seg[io]->getNTimes(tr)); } else { tracksOuterMdc.push_back(seg[io]->getTrack(tr)); weightOuterMdc.push_back(seg[io]->getNTimes(tr)); } } } if(print){ if(seg[0]){ cout<<"seg : "<<0<<" , ntr "<getRichInd()>-1){ HRichHitSim* richHit = (HRichHitSim*) HParticleTool::getRichHit(cand->getRichInd()); if(richHit){ if(richHit->track1>0) { tracksRich.push_back(richHit->track1); weightRich.push_back(richHit->weigTrack1); } if(richHit->track2>0) { tracksRich.push_back(richHit->track2); weightRich.push_back(richHit->weigTrack2); } if(richHit->track3>0) { tracksRich.push_back(richHit->track3); weightRich.push_back(richHit->weigTrack3); } if(print){ cout<<"Rich : "<track1 <<" ntimes "<weigTrack1<track2 <<" ntimes "<weigTrack2<track3 <<" ntimes "<weigTrack3<isTofHitUsed()) { tofhit =(HTofHitSim*) HParticleTool::getTofHit(cand->getTofHitInd()); if(tofhit){ tracksMeta.push_back(tofhit->getNTrack1()); tracksMeta.push_back(tofhit->getNTrack2()); if(print){ cout<<"TOF hit : "<getNTrack1()<<" "<getNTrack2() <isTofClstUsed()){ tofclst=(HTofClusterSim*)HParticleTool::getTofCluster(cand->getTofClstInd()); if(tofclst){ if(print){ cout<<"TOF clst : "<getNTrack1()<<" "<getNTrack2() <getNTrack1()); tracksMeta.push_back(tofclst->getNTrack2()); } else { if(print)cout<<"HParticleTool::getSimTracks() : retrieved zero pointer of tofclst!"<isRpcClstUsed()){ Int_t tracks[4]; rpcclst=(HRpcClusterSim*)HParticleTool::getRpcCluster(cand->getRpcInd()); if(rpcclst){ rpcclst->getTrackList(tracks); for(Int_t i=0;i<4;i++){ if(tracks[i]>0)tracksMeta.push_back(tracks[i]); } if(print){ cout<<"RPC clst : "<getShowerInd()>-1) { shower=(HShowerHitSim*)HParticleTool::getShowerHit(cand->getShowerInd()); if(shower){ for(Int_t i=0;igetNTracks();i++){ tracksShower.push_back(shower->getTrack(i)); if(cand->isShowerUsed()) tracksMeta.push_back(shower->getTrack(i)); } if(print) { cout<<"ShowerHit : "<getTrack(i)<<" "<getInnerSegInd()); seg [1]= (HMdcSegSim*) HParticleTool::getMdcSeg(cand->getOuterSegInd()); for(Int_t io=0;io<2;io++){ if(seg[io]==0) continue; cout<<"seg : "<getNTracks()<<" : "<getNTracks(); tr++){ cout<<" "<getTrack(tr)<<" ntimes "<getNTimes(tr)<getRichInd()>-1){ HRichHitSim* richHit = (HRichHitSim*) HParticleTool::getRichHit(cand->getRichInd()); if(richHit){ cout<<"Rich : "<track1 <<" ntimes "<weigTrack1<track2 <<" ntimes "<weigTrack2<track3 <<" ntimes "<weigTrack3<isTofHitUsed()) { tofhit =(HTofHitSim*) HParticleTool::getTofHit(cand->getTofHitInd()); if(tofhit){ cout<<"TOF hit : "<getNTrack1()<<" "<getNTrack2() <isTofClstUsed()){ tofclst=(HTofClusterSim*)HParticleTool::getTofCluster(cand->getTofClstInd()); if(tofclst){ cout<<"TOF clst : "<getNTrack1()<<" "<getNTrack2() <isRpcClstUsed()){ Int_t tracks[4]; rpcclst=(HRpcClusterSim*)HParticleTool::getRpcCluster(cand->getRpcInd()); if(rpcclst){ rpcclst->getTrackList(tracks); cout<<"RPC clst : "<getShowerInd()>-1){ shower = (HShowerHitSim*)HParticleTool::getShowerHit(cand->getShowerInd()); if(shower){ cout<<"ShowerHit : "<getNTracks();i++){ cout<getTrack(i)<<" "<getRichInd() >= 0 && cand2->getRichInd() >= 0) || (cand1->getRichInd() >= 0 && cand2->getRichInd() < 0)){ if( cand1->getRichInd() == cand2->getRichInd() ) flag = kSameRICH; else flag = kNoSameRICH; } if(cand1->getInnerSegInd() >= 0 && cand2->getInnerSegInd() >= 0){ if( cand1->getInnerSegInd() == cand2->getInnerSegInd() ) flag = flag|kSameInnerMDC; else flag = flag|kNoSameInnerMDC; } if(cand1->getOuterSegInd() >= 0 && cand2->getOuterSegInd() >= 0 ){ if( cand1->getOuterSegInd() == cand2->getOuterSegInd() ) flag = flag|kSameOuterMDC; else flag = flag|kNoSameOuterMDC; } if(cand1->getSystemUsed() != -1 && cand2->getSystemUsed() != -1 ){ if( cand1->isTofHitUsed() && cand2->isTofHitUsed()){ if( cand1->getTofHitInd() == cand2->getTofHitInd() ) flag = flag|kSameMETA; else flag = flag|kNoSameMETA; } else if(cand1->isTofClstUsed() && cand2->isTofClstUsed() ){ if(cand1->getTofClstInd() == cand2->getTofClstInd() ) flag = flag|kSameMETA; else flag = flag|kNoSameMETA; } else if(cand1->isRpcClstUsed() && cand2->isRpcClstUsed() ) { if (cand1->getRpcInd() == cand2->getRpcInd()) flag = flag|kSameMETA; else flag = flag|kNoSameMETA; } else if(cand1->isShowerUsed() && cand2->isShowerUsed()){ if(cand1->getShowerInd() == cand2->getShowerInd() ) flag = flag|kSameMETA; else flag = flag|kNoSameMETA; } else { flag = flag|kNoSameMETA;} } if(cand1->getCharge() != 0 && cand2->getCharge() != 0){ if (cand1->getCharge() > 0 && cand2->getCharge() > 0) flag = flag|kSamePosPolarity|kSamePolarity; else if(cand1->getCharge() < 0 && cand2->getCharge() < 0) flag = flag|kSameNegPolarity|kSamePolarity; else if(cand1->getCharge() < 0 && cand2->getCharge() > 0) flag = flag|kNoSamePosPolarity|kNoSameNegPolarity|kNoSamePolarity; else if(cand1->getCharge() > 0 && cand2->getCharge() < 0) flag = flag|kNoSamePosPolarity|kNoSameNegPolarity|kNoSamePolarity; } } //---------------------------------------------------------------------- if(cand2->getRichInd() >= 0) flag = flag|kRICH2; else flag = flag|kNoRICH2; if(cand2->getInnerSegmentChi2() >= 0) flag = flag|kFittedInnerMDC2; else flag = flag|kNoFittedInnerMDC2; if(cand2->getOuterSegmentChi2() >= 0) flag = flag|kFittedOuterMDC2; else flag = flag|kNoFittedOuterMDC2; if(cand2->getOuterSegInd() >= 0 ) flag = flag|kOuterMDC2; else flag = flag|kNoOuterMDC2; if(cand2->getChi2() >= 0 && cand2->getChi2() < 1e6) flag = flag|kRK2; else flag = flag|kNoRK2; if(cand2->getSystemUsed() != -1 ) flag = flag|kMETA2; else flag = flag|kNoMETA2; if(cand2->isFlagBit(kIsLepton)) flag = flag|kIsLepton2; else flag = flag|kNoIsLepton2; if(cand2->isFlagBit(kIsUsed)) flag = flag|kIsUsed2; else flag = flag|kNoIsUsed2; return kTRUE; } Bool_t HParticleTool::evalPairsFlags(UInt_t flag,UInt_t fl) { // evalutates the bitwise flags defined in flag (see eClosePairSelect + ePairCase (hparticledef.h) // against fl. return kTRUE if all set bit in flag are set in fl too. /* if( (flag&kSameRICH) == kSameRICH && (fl&kSameRICH) != kSameRICH ) return kFALSE; if( (flag&kSameInnerMDC) == kSameInnerMDC && (fl&kSameInnerMDC) != kSameInnerMDC) return kFALSE; if( (flag&kSameOuterMDC) == kSameOuterMDC && (fl&kSameOuterMDC) != kSameOuterMDC) return kFALSE; if( (flag&kSameMETA) == kSameMETA && (fl&kSameMETA) != kSameMETA ) return kFALSE; if( (flag&kSamePosPolarity) == kSamePosPolarity && (fl&kSamePosPolarity) != kSamePosPolarity ) return kFALSE; if( (flag&kSameNegPolarity) == kSameNegPolarity && (fl&kSameNegPolarity) != kSameNegPolarity ) return kFALSE; if( (flag&kSamePolarity) == kSamePolarity && (fl&kSamePolarity) != kSamePolarity ) return kFALSE; if( (flag&kRICH2) == kRICH2 && (fl&kRICH2) != kRICH2 ) return kFALSE; if( (flag&kFittedInnerMDC2) == kFittedInnerMDC2 && (fl&kFittedInnerMDC2) != kFittedInnerMDC2 ) return kFALSE; if( (flag&kFittedOuterMDC2) == kFittedOuterMDC2 && (fl&kFittedOuterMDC2) != kFittedOuterMDC2 ) return kFALSE; if( (flag&kOuterMDC2) == kOuterMDC2 && (fl&kOuterMDC2) != kOuterMDC2 ) return kFALSE; if( (flag&kRK2) == kRK2 && (fl&kRK2) != kRK2 ) return kFALSE; if( (flag&kMETA2) == kMETA2 && (fl&kMETA2) != kMETA2 ) return kFALSE; if( (flag&kIsLepton2) == kIsLepton2 && (fl&kIsLepton2) != kIsLepton2 ) return kFALSE; if( (flag&kIsUsed2) == kIsUsed2 && (fl&kIsUsed2) != kIsUsed2 ) return kFALSE; if( (flag&kNoSameRICH) == kNoSameRICH && (fl&kNoSameRICH) != kNoSameRICH ) return kFALSE; if( (flag&kNoSameInnerMDC) == kNoSameInnerMDC && (fl&kNoSameInnerMDC) != kNoSameInnerMDC) return kFALSE; if( (flag&kNoSameOuterMDC) == kNoSameOuterMDC && (fl&kNoSameOuterMDC) != kNoSameOuterMDC) return kFALSE; if( (flag&kNoSameMETA) == kNoSameMETA && (fl&kNoSameMETA) != kNoSameMETA ) return kFALSE; if( (flag&kNoSamePosPolarity) == kNoSamePosPolarity && (fl&kNoSamePosPolarity) != kNoSamePosPolarity ) return kFALSE; if( (flag&kNoSameNegPolarity) == kNoSameNegPolarity && (fl&kNoSameNegPolarity) != kNoSameNegPolarity ) return kFALSE; if( (flag&kNoSamePolarity) == kNoSamePolarity && (fl&kNoSamePolarity) != kNoSamePolarity ) return kFALSE; if( (flag&kNoRICH2) == kNoRICH2 && (fl&kNoRICH2) != kNoRICH2 ) return kFALSE; if( (flag&kNoFittedInnerMDC2) == kNoFittedInnerMDC2 && (fl&kNoFittedInnerMDC2) != kNoFittedInnerMDC2 ) return kFALSE; if( (flag&kNoFittedOuterMDC2) == kNoFittedOuterMDC2 && (fl&kNoFittedOuterMDC2) != kNoFittedOuterMDC2 ) return kFALSE; if( (flag&kNoOuterMDC2) == kNoOuterMDC2 && (fl&kNoOuterMDC2) != kNoOuterMDC2 ) return kFALSE; if( (flag&kNoRK2) == kNoRK2 && (fl&kNoRK2) != kNoRK2 ) return kFALSE; if( (flag&kNoMETA2) == kNoMETA2 && (fl&kNoMETA2) != kNoMETA2 ) return kFALSE; if( (flag&kNoIsLepton2) == kNoIsLepton2 && (fl&kNoIsLepton2) != kNoIsLepton2 ) return kFALSE; if( (flag&kNoIsUsed2) == kNoIsUsed2 && (fl&kNoIsUsed2) != kNoIsUsed2 ) return kFALSE; if( (flag&kNoUseRICH) == kNoUseRICH && (fl&kNoUseRICH) != kNoUseRICH ) return kFALSE; */ if ( (flag&fl) == flag ) return kTRUE; else return kFALSE; // return kTRUE; } Bool_t HParticleTool::isPairsFlagsBit(UInt_t flag,UInt_t fl) { // evalutates the bitwise flags defined fl (see eClosePairSelect + ePairCase (hparticledef.h)) // return kTRUE if all set bits in fl are set in flags too if ( (flag&fl) == fl ) return kTRUE; else return kFALSE; } Bool_t HParticleTool::evalPairsFlags(UInt_t flag,HParticleCand* cand1,HParticleCand* cand2) { // evalutates the bitwise flags defined in eClosePairSelect + ePairCase (hparticledef.h) // return kTRUE if all set bit are identical. cand1 as reference for Pair. if flag == 0 // result will be kTRUE if (flag == 0) return kTRUE; UInt_t fl = 0; HParticleTool::setPairFlags(fl,cand2,cand1); return evalPairsFlags(flag,fl); } Bool_t HParticleTool::evalPairsFlags(UInt_t flag,HParticlePair& pair) { // evalutates the bitwise flags defined in eClosePairSelect + ePairCase (hparticledef.h) // return kTRUE if all set bit are identical . if flag == 0 // result will be kTRUE if (flag == 0) return kTRUE; UInt_t fl = pair.getPairFlags(); return evalPairsFlags(flag,fl); return kTRUE; } Bool_t HParticleTool::evalPairsFlags(vector& flags,vector& results,HParticlePair& pair) { // evalutates the bitwise flags defined in eClosePairSelect + ePairCase (hparticledef.h) // return kTRUE if one of the flags in vector flags are identical. The result for // each individual flag is store in vector results (which will be automaticly cleared before). results.clear(); Bool_t val =kFALSE; Bool_t val2=kFALSE; for(UInt_t i = 0 ; i < flags.size(); i ++){ val = HParticleTool::evalPairsFlags(flags[i],pair); if(val) val2 = kTRUE; results.push_back(val); } return val2; } Bool_t HParticleTool::evalPairsFlags(vector& flags,HParticlePair& pair) { // evalutates the bitwise flags defined in eClosePairSelect + ePairCase (hparticledef.h) // return kTRUE if one of the flags in vector flags are identical. for(UInt_t i = 0 ; i < flags.size(); i ++){ if(HParticleTool::evalPairsFlags(flags[i],pair)) return kTRUE; } return kFALSE; } Bool_t HParticleTool::isGoodClusterVertex(Float_t minZ) { // request good chi2 and a z pos > minZ (apr12 typically -65 mm, jul14/aug14 -200) // minZ is given in mm // if minZ > 0 gHades is used to retrieve the beamTimeID default values Float_t minCut = minZ; if(gHades && minZ > 0){ // try to get it from gHades if(gHades->getBeamTimeID()==Particle::kUnknownBeam) { ::Error("isGoodClusterVertex()","minZ > 0 but no beamTimeID set in gHades"); return kFALSE; } else if (gHades->getBeamTimeID()==Particle::kApr12) { minCut = -65; } else if (gHades->getBeamTimeID()==Particle::kJul14||gHades->getBeamTimeID()==Particle::kAug14) { minCut = -200; } } else if(!gHades && minZ > 0) { ::Error("isGoodClusterVertex()","minZ > 0 but gHades==NULL"); return kFALSE; } HVertex& vertexClust = gHades->getCurrentEvent()->getHeader()->getVertexCluster(); if(vertexClust.getChi2()>-1 && vertexClust.getZ()>minCut) return kTRUE; return kFALSE; } Bool_t HParticleTool::isGoodRecoVertex(Float_t minZ) { // request good chi2 and a z pos > minZ (apr12 typically -65 mm) // minZ is given in mm. For good chi2 at least 2 particle cands are needed. // if minZ > 0 gHades is used to retrieve the beamTimeID default values Float_t minCut = minZ; if(gHades && minZ > 0){ // try to get it from gHades if(gHades->getBeamTimeID()==Particle::kUnknownBeam) { ::Error("isGoodRecoVertex()","minZ > 0 but no beamTimeID set in gHades"); return kFALSE; } else if (gHades->getBeamTimeID()==Particle::kApr12) { minCut = -65; } else if (gHades->getBeamTimeID()==Particle::kJul14||gHades->getBeamTimeID()==Particle::kAug14) { minCut = -200; } } else if(!gHades && minZ > 0) { ::Error("isGoodRecoVertex()","minZ > 0 but gHades==NULL"); return kFALSE; } HVertex& vertexReco = gHades->getCurrentEvent()->getHeader()->getVertexReco(); if(vertexReco.getChi2()>-1 && vertexReco.getZ()>minCut) return kTRUE; return kFALSE; } Bool_t HParticleTool::isGoodSTART(Int_t minFlag) { // request existing start hit object and starthit->getCorrFlag()>=minFlag // REMARK : corrFlags are filled by HParticleStart2HitF (-1,0,1,2). // In this case HStart2Hit::getIteration() (0,1,2) gives 2. For Iteration // <2 corrFlags are not checked HStart2Hit* starthit = 0; starthit = HCategoryManager::getObject(starthit,catStart2Hit,0,kTRUE); if(starthit && starthit->getIteration() > 1 && starthit->getCorrFlag()>=minFlag) return kTRUE; if(starthit && starthit->getIteration() < 2) return kTRUE; return kFALSE; } Bool_t HParticleTool::isNoVETO(Float_t minStart,Float_t maxStart) { // return kTRUE if no VETO hit inside timewindow minStart-maxStart [ns] // was found. otherwise kFALSE. // requires catStart2Cal // return kFALSE if category is not available HStart2Cal* pCal = NULL; HCategory* fCatCal = HCategoryManager::getCategory(catStart2Cal,-1); if(!fCatCal) return kFALSE; Bool_t notfound = kTRUE; for (Int_t entry = 0; entry < fCatCal->getEntries(); ++entry) { if (NULL == (pCal = static_cast(fCatCal->getObject(entry)))) { return kFALSE; } if (3 != pCal->getModule()) { continue; } for (Int_t i = 0; i < pCal->getMultiplicity() && i < pCal->getMaxMultiplicity(); ++i) { Float_t t = pCal->getTime(i + 1); if(t < minStart || t > maxStart ) continue; notfound = kFALSE; } } return notfound; } Bool_t HParticleTool::isNoSTART(Float_t minStart,Float_t maxStart,Float_t window) { // return kTRUE if no START hit inside timewindow minStart-maxStart [ns] of the // START time was found. otherwise kFALSE. (startime +- window excluded) // requires catStart2Hit and catStart2Cal // return kFALSE if categories are not available HStart2Hit* starthit = 0; starthit = HCategoryManager::getObject(starthit,catStart2Hit,0,kTRUE); if(!starthit) return kFALSE; Float_t starttime = starthit->getTime(); window = fabs(window); HStart2Cal* pCal = NULL; HCategory* fCatCal = HCategoryManager::getCategory(catStart2Cal,-1); if(!fCatCal) return kFALSE; Bool_t notfound = kTRUE; for (Int_t entry = 0; entry < fCatCal->getEntries(); ++entry) { if (NULL == (pCal = static_cast(fCatCal->getObject(entry)))) { return kFALSE; } if (pCal->getModule()>1) { continue; } for (Int_t i = 0; i < pCal->getMultiplicity() && i < pCal->getMaxMultiplicity(); ++i) { Float_t t = pCal->getTime(i + 1); if(fabs(starttime-t) maxStart ) continue; notfound = kFALSE; } } return notfound; } Bool_t HParticleTool::isGoodSTARTVETO(Float_t minStart,Float_t maxStart, Float_t window) { // return kTRUE is no starthit between minStart [ns] and maxStart [ns] // has been found which has no correlation with VETO inside window [ns] // +-window arround starthit startime is excluded. // requires catStart2Hit + catStart2Cal category. If the categories are // not available kFALSE is returned. HCategory* fCatCal = HCategoryManager::getCategory(catStart2Cal,-1); if(!fCatCal) return kFALSE; window=fabs(window); HStart2Hit* starthit = 0; starthit = HCategoryManager::getObject(starthit,catStart2Hit,0,kTRUE); if(!starthit) return kFALSE; Float_t starttime = starthit->getTime(); Bool_t kDontUseEvent = kFALSE; HStart2Cal* pCal = NULL; HStart2Cal* pVeto = NULL; //------------------------------------------------------- for (Int_t entry = 0; entry < fCatCal->getEntries(); ++entry) { // loop start if (NULL == (pCal = static_cast(fCatCal->getObject(entry)))) { return kFALSE; } if (pCal->getModule()>1) { continue; } // only module 0 and 1 (no veto or others) for (Int_t i = 0; i < pCal->getMultiplicity() && i < pCal->getMaxMultiplicity(); ++i) { // loop all hits of this strip Float_t t = pCal->getTime(i + 1); if(fabs(t-starttime) < window) continue; if(tmaxStart) continue; Bool_t isCorrelated = kFALSE; //------------------------------------------------------- for (Int_t entry2 = 0; entry2 < fCatCal->getEntries(); ++entry2) { // loop all veto if (NULL == (pVeto = static_cast(fCatCal->getObject(entry2)))) { return kFALSE; } if (pVeto->getModule()!=3) { continue; } // only veto for (Int_t j = 0; j < pVeto->getMultiplicity() && j < pVeto->getMaxMultiplicity(); ++j) { // loop all hits of this strip Float_t tveto = pVeto->getTime(j + 1); if(tvetomaxStart || (t<540&&t>520)) continue; // exclude reference time if( fabs(t - tveto ) < window ) { isCorrelated = kTRUE; } } } if(!isCorrelated) kDontUseEvent = kTRUE; // //------------------------------------------------------- } } //------------------------------------------------------- if(kDontUseEvent) return kFALSE; else return kTRUE; } Bool_t HParticleTool::isGoodSTARTMETA(Float_t minStart,Float_t maxStart,Int_t tresh,Float_t window, Float_t offset) { // search for metahits inside the event which are correlated with // late starthits. // return kTRUE if in range minStart [ns] to maxStart [ns] of START Hits // not more than tresh metahits inside correlation window [ns] are found. // offset [ns](apr12 = 7) has to be adjusted to move the fastest particles // of good triggered event to 0. // requires catRpcHit, catTofHit, catStart2Cal. if any of those is missing kFALSE // is returned. HCategory* rpcCat = HCategoryManager::getCategory(catRpcHit,-1); if(!rpcCat) return kFALSE; HCategory* tofCat = HCategoryManager::getCategory(catTofHit,-1); if(!tofCat) return kFALSE; HCategory* fCatCal = HCategoryManager::getCategory(catStart2Cal,-1); if(!fCatCal) return kFALSE; Bool_t kDontUseEvent = kFALSE; Int_t ctBad = 0; HStart2Cal* pCal = NULL; for(Int_t j=0;jgetEntries();j++){ // rpc loop HRpcHit* rpc = (HRpcHit*)rpcCat->getObject(j); if(!rpc) return kFALSE; //------------------------------------------------------- for (Int_t entry = 0; entry < fCatCal->getEntries(); ++entry) { if (NULL == (pCal = static_cast(fCatCal->getObject(entry)))) { return kFALSE; } if (pCal->getModule()>1) { continue; } // only module 0 and 1 (no veto or others) for (Int_t i = 0; i < pCal->getMultiplicity() && i < pCal->getMaxMultiplicity(); ++i) { // loop all hits of this strip Float_t t = pCal->getTime(i + 1); if(tmaxStart || (t<540&&t>520)) continue; // exclude reference time if( t - rpc->getTof() - offset < window ) { ctBad ++; kDontUseEvent = kTRUE;} } } //------------------------------------------------------- } for(Int_t j=0;jgetEntries();j++){ // tof loop HTofHit* tof = (HTofHit*)tofCat->getObject(j); if(!tof) return kFALSE; //------------------------------------------------------- for (Int_t entry = 0; entry < fCatCal->getEntries(); ++entry) { if (NULL == (pCal = static_cast(fCatCal->getObject(entry)))) { return kFALSE; } if (pCal->getModule()>1) { continue; } // only module 0 and 1 (no veto or others) for (Int_t i = 0; i < pCal->getMultiplicity() && i < pCal->getMaxMultiplicity(); ++i) { // loop all hits of this strip Float_t t = pCal->getTime(i + 1); if(tmaxStart || (t<540&&t>520)) continue; // exclude reference time if( t - tof->getTof() - offset < window ) { ctBad ++; kDontUseEvent = kTRUE;} } } //------------------------------------------------------- } if(kDontUseEvent&&ctBad>tresh) {return kFALSE;} else return kTRUE; } Bool_t HParticleTool::isNoSTARTPileUp() { // request existing start hit object and starthit->getMultiplicity()==1 HStart2Hit* starthit = 0; starthit = HCategoryManager::getObject(starthit,catStart2Hit,0,kTRUE); if(starthit && starthit->getMultiplicity()==1) return kTRUE; return kFALSE; } Bool_t HParticleTool::isNoMETAPileUp(Float_t ftimeTofCut,Int_t threshold) { // loops TofHit and RpcCluster and counds objects with tof < 0 // or tof > ftimeTofCut. if n > threshold (pileup) kFALSE is returned // if tofhit and rpccluster cateories are not available kTRUE is returned // if threshold < 0 gHades is used to retrieve the beamTimeID default values Int_t thresh = threshold; if(gHades && threshold < 0){ // try to get it from gHades if(gHades->getBeamTimeID()==Particle::kUnknownBeam) { ::Error("isNoMETAPileUp()","threshold < 0 but no beamTimeID set in gHades"); return kFALSE; } else if (gHades->getBeamTimeID()==Particle::kApr12) { thresh = 5*7; } else if (gHades->getBeamTimeID()==Particle::kJul14||gHades->getBeamTimeID()==Particle::kAug14) { thresh = 7; } } else if(!gHades && threshold < 0) { ::Error("isNoMETAPileUp()","threshold < 0 but gHades==NULL"); return kFALSE; } HCategory* tofhitCat = HCategoryManager::getCategory(catTofHit,-1); HCategory* rpcclustCat = HCategoryManager::getCategory(catRpcCluster,-1); if(tofhitCat&&rpcclustCat){ HTofHit* tofhit=0; HRpcCluster* rpccluster=0; Int_t N=0; Int_t n = tofhitCat ->getEntries(); for(Int_t i=0;igetTof() > ftimeTofCut && tofhit->getTof() < 0 ) N++; } n = rpcclustCat ->getEntries(); for(Int_t i=0;igetTof() > ftimeTofCut && rpccluster->getTof() < 0 ) N++; } if(N>thresh) return kFALSE; else return kTRUE; } else return kTRUE; // if there is nothing to check its ok :-) } Bool_t HParticleTool::isGoodTrigger(Int_t triggerbit) { // test if the triggerbit is set ( bit number is pt (physics trigger) + 10) // example : pt3 -> 13 // if triggerbit < 0 gHades is used to retrieve the beamTimeID default values if(triggerbit > 0 && triggerbit < 10){ ::Error("isGoodTrigger()","triggerbit has to be larger 10 (pt+10)"); return kFALSE; } Int_t tbit = triggerbit; if(gHades && triggerbit < 0){ // try to get it from gHades if(gHades->getBeamTimeID()==Particle::kUnknownBeam) { ::Error("isGoodTrigger()","triggerbit < 0 but no beamTimeID set in gHades"); return kFALSE; } else if (gHades->getBeamTimeID()==Particle::kApr12) { tbit = 13; } else if (gHades->getBeamTimeID()==Particle::kJul14||gHades->getBeamTimeID()==Particle::kAug14) { tbit = 11; } } else if(!gHades && triggerbit < 0) { ::Error("isGoodTrigger()","triggerbit < 0 but gHades==NULL"); return kFALSE; } if(gHades->getCurrentEvent()->getHeader()->isTBit(tbit)) return kTRUE; return kFALSE; } Int_t HParticleTool::getTofHitMult(Float_t minTof,Float_t maxTof, Int_t* sector) { // return the TofHit multiplicity between minTof and maxTof. if sector // (pointer to Int_t sector[6]) is not NULL sectorwise mult is filled. // the array is rest to 0 before. requires catTofHit HCategory* tofhitCat = HCategoryManager::getCategory(catTofHit); Int_t ct = 0; if(tofhitCat) { if(sector) for(Int_t i=0;i<6;i++) sector[i] =0; Int_t size = tofhitCat->getEntries(); HTofHit* hit=0; for(Int_t j = 0; j < size; j++) { hit = HCategoryManager::getObject(hit,tofhitCat,j); if(hit) { if(hit->getTof()>minTof&&hit->getTof()getSector()]++; } } } } return ct; } Int_t HParticleTool::getRpcHitMult(Float_t minTof,Float_t maxTof, Int_t* sector) { // return the RpcHit multiplicity between minTof and maxTof. if sector // (pointer to Int_t sector[6]) is not NULL sectorwise mult is filled. // the array is rest to 0 before. requires catRpcHit. HCategory* rpchitCat = HCategoryManager::getCategory(catRpcHit); Int_t ct = 0; if(rpchitCat) { if(sector) for(Int_t i=0;i<6;i++) sector[i] =0; Int_t size = rpchitCat->getEntries(); HRpcHit* hit=0; for(Int_t j = 0; j < size; j++) { hit = HCategoryManager::getObject(hit,rpchitCat,j); if(hit) { if(hit->getTof()>minTof&&hit->getTof()getSector()]++; } } } } return ct; }