diff --git a/analysis/data/v0_projection_2016_config.json b/analysis/data/v0_projection_2016_config.json new file mode 100644 index 000000000..fc24801bb --- /dev/null +++ b/analysis/data/v0_projection_2016_config.json @@ -0,0 +1,586 @@ +{ + "7629": { + "target_position": -4.3, + "rotated_mean_x": -0.24175124390756064, + "rotated_mean_y": -0.0841215989869986, + "rotated_sigma_x": 0.3209415790391567, + "rotated_sigma_y": 0.09259800976999152, + "rotation_angle_mrad": -50.36137375374897 + }, + "7630": { + "target_position": -4.3, + "rotated_mean_x": -0.10853070867807432, + "rotated_mean_y": -0.05681037797761715, + "rotated_sigma_x": 0.3184934749199056, + "rotated_sigma_y": 0.09275885589498911, + "rotation_angle_mrad": -46.23035117234371 + }, + "7636": { + "target_position": -4.3, + "rotated_mean_x": -0.22082700852206247, + "rotated_mean_y": -0.06668426373830501, + "rotated_sigma_x": 0.33086701088023096, + "rotated_sigma_y": 0.09363436968355388, + "rotation_angle_mrad": -53.320135209273246 + }, + "7644": { + "target_position": -4.3, + "rotated_mean_x": -0.23697406235152177, + "rotated_mean_y": -0.06651318077968818, + "rotated_sigma_x": 0.33757814486100035, + "rotated_sigma_y": 0.10379397336685675, + "rotation_angle_mrad": -38.73772989358739 + }, + "7653": { + "target_position": -4.3, + "rotated_mean_x": -0.2672865930009524, + "rotated_mean_y": -0.0783603982849895, + "rotated_sigma_x": 0.327967438536696, + "rotated_sigma_y": 0.10825572390387936, + "rotation_angle_mrad": -41.27264008425549 + }, + "7779": { + "target_position": -4.3, + "rotated_mean_x": -0.24380906864696564, + "rotated_mean_y": -0.09617371007194137, + "rotated_sigma_x": 0.35776447917545856, + "rotated_sigma_y": 0.09923847113497254, + "rotation_angle_mrad": -75.32732430212002 + }, + "7780": { + "target_position": -4.3, + "rotated_mean_x": -0.23723711199085368, + "rotated_mean_y": -0.07483333486097161, + "rotated_sigma_x": 0.35328279864972684, + "rotated_sigma_y": 0.10378725985857315, + "rotation_angle_mrad": -93.8241952328115 + }, + "7781": { + "target_position": -4.3, + "rotated_mean_x": -0.268892040569693, + "rotated_mean_y": -0.07389749686711232, + "rotated_sigma_x": 0.348671667104869, + "rotated_sigma_y": 0.09748216858981719, + "rotation_angle_mrad": -59.51154767830512 + }, + "7795": { + "target_position": -4.3, + "rotated_mean_x": -0.24226850752216486, + "rotated_mean_y": -0.03588097201649974, + "rotated_sigma_x": 0.36648235003371576, + "rotated_sigma_y": 0.11199775586677405, + "rotation_angle_mrad": -122.02331995866206 + }, + "7796": { + "target_position": -4.3, + "rotated_mean_x": -0.2423345878393665, + "rotated_mean_y": -0.03704949499316914, + "rotated_sigma_x": 0.3649364391153848, + "rotated_sigma_y": 0.11445780580489512, + "rotation_angle_mrad": -122.26673570532898 + }, + "7798": { + "target_position": -4.3, + "rotated_mean_x": -0.24470551057636938, + "rotated_mean_y": -0.05237861787564785, + "rotated_sigma_x": 0.35507213295707557, + "rotated_sigma_y": 0.10166519981120302, + "rotation_angle_mrad": -81.85556948547814 + }, + "7799": { + "target_position": -4.3, + "rotated_mean_x": -0.2518347753842008, + "rotated_mean_y": -0.052743138592004885, + "rotated_sigma_x": 0.3502637252987632, + "rotated_sigma_y": 0.10101105046495297, + "rotation_angle_mrad": -77.64193977501824 + }, + "7800": { + "target_position": -4.3, + "rotated_mean_x": -0.24891375757666734, + "rotated_mean_y": -0.05671012038956239, + "rotated_sigma_x": 0.35043148748833236, + "rotated_sigma_y": 0.09905029714140397, + "rotation_angle_mrad": -76.1354756262848 + }, + "7801": { + "target_position": -4.3, + "rotated_mean_x": -0.251434773460241, + "rotated_mean_y": -0.06099919245002469, + "rotated_sigma_x": 0.3480943234084851, + "rotated_sigma_y": 0.10096012304619151, + "rotation_angle_mrad": -75.40502474179134 + }, + "7803": { + "target_position": -4.3, + "rotated_mean_x": -0.24797182708752025, + "rotated_mean_y": -0.061326620157149285, + "rotated_sigma_x": 0.3298914215023227, + "rotated_sigma_y": 0.09407409065230352, + "rotation_angle_mrad": -50.95765333768275 + }, + "7804": { + "target_position": -4.3, + "rotated_mean_x": -0.2436064851860926, + "rotated_mean_y": -0.05395052806164744, + "rotated_sigma_x": 0.34412221378521685, + "rotated_sigma_y": 0.09772663762988872, + "rotation_angle_mrad": -57.28676904206268 + }, + "7805": { + "target_position": -4.3, + "rotated_mean_x": -0.2394473851290572, + "rotated_mean_y": -0.0590040451933412, + "rotated_sigma_x": 0.34600677878777497, + "rotated_sigma_y": 0.09591707989032405, + "rotation_angle_mrad": -52.26677381181967 + }, + "7807": { + "target_position": -4.3, + "rotated_mean_x": -0.24422937714665913, + "rotated_mean_y": -0.05828467429563133, + "rotated_sigma_x": 0.3432714788252632, + "rotated_sigma_y": 0.09840663987070958, + "rotation_angle_mrad": -61.00519136446077 + }, + "7947": { + "target_position": -4.3, + "rotated_mean_x": -0.1903308426581071, + "rotated_mean_y": -0.1089242312910535, + "rotated_sigma_x": 0.3635209656825039, + "rotated_sigma_y": 0.09356023937472986, + "rotation_angle_mrad": -39.87596092824682 + }, + "7948": { + "target_position": -4.3, + "rotated_mean_x": -0.15636580520466378, + "rotated_mean_y": -0.13148426563936258, + "rotated_sigma_x": 0.35320893511581075, + "rotated_sigma_y": 0.0908365556186038, + "rotation_angle_mrad": -40.1750923681627 + }, + "7949": { + "target_position": -4.3, + "rotated_mean_x": -0.14016448990153968, + "rotated_mean_y": -0.1441567576282953, + "rotated_sigma_x": 0.34798055523942617, + "rotated_sigma_y": 0.09237011117447004, + "rotation_angle_mrad": -42.87868623621624 + }, + "7962": { + "target_position": -4.3, + "rotated_mean_x": -0.18567028219050677, + "rotated_mean_y": -0.13195225127295007, + "rotated_sigma_x": 0.33700814119896694, + "rotated_sigma_y": 0.0889257280943231, + "rotation_angle_mrad": -54.64889950966569 + }, + "7963": { + "target_position": -4.3, + "rotated_mean_x": -0.15127703461686315, + "rotated_mean_y": -0.15448316457917208, + "rotated_sigma_x": 0.3387768673485481, + "rotated_sigma_y": 0.08884331459065495, + "rotation_angle_mrad": -47.385123919361774 + }, + "7964": { + "target_position": -4.3, + "rotated_mean_x": -0.15417283472555077, + "rotated_mean_y": -0.1497868401638537, + "rotated_sigma_x": 0.3454616820299382, + "rotated_sigma_y": 0.09086326762111298, + "rotation_angle_mrad": -43.35758205080592 + }, + "7965": { + "target_position": -4.3, + "rotated_mean_x": -0.16558633463586378, + "rotated_mean_y": -0.15576444083370428, + "rotated_sigma_x": 0.34479882872373635, + "rotated_sigma_y": 0.08982786312798145, + "rotation_angle_mrad": -44.830498263166774 + }, + "7966": { + "target_position": -4.3, + "rotated_mean_x": -0.1832527197951328, + "rotated_mean_y": -0.14596819580077122, + "rotated_sigma_x": 0.3420594060808304, + "rotated_sigma_y": 0.08971155196407041, + "rotation_angle_mrad": -37.12323829087083 + }, + "7968": { + "target_position": -4.3, + "rotated_mean_x": -0.17626353514005402, + "rotated_mean_y": -0.14212842930360925, + "rotated_sigma_x": 0.3477943387704159, + "rotated_sigma_y": 0.0921444336063571, + "rotation_angle_mrad": -38.830575766492856 + }, + "7970": { + "target_position": -4.3, + "rotated_mean_x": -0.19069571530441323, + "rotated_mean_y": -0.1291018212809842, + "rotated_sigma_x": 0.3480725230030799, + "rotated_sigma_y": 0.09381566106973012, + "rotation_angle_mrad": -33.52511338975475 + }, + "7972": { + "target_position": -4.3, + "rotated_mean_x": -0.05146698275350976, + "rotated_mean_y": -0.10367834196692258, + "rotated_sigma_x": 0.33964286691913687, + "rotated_sigma_y": 0.08904573026599312, + "rotation_angle_mrad": -34.13615441277002 + }, + "7976": { + "target_position": -4.3, + "rotated_mean_x": -0.05733850772891847, + "rotated_mean_y": -0.1031188264849115, + "rotated_sigma_x": 0.3165173747733235, + "rotated_sigma_y": 0.08694483844997235, + "rotation_angle_mrad": -42.92163705627083 + }, + "7983": { + "target_position": -4.3, + "rotated_mean_x": -0.05791057317283221, + "rotated_mean_y": -0.1093756160554428, + "rotated_sigma_x": 0.3320181229782187, + "rotated_sigma_y": 0.08980551606263029, + "rotation_angle_mrad": -41.57561029730484 + }, + "7984": { + "target_position": -4.3, + "rotated_mean_x": -0.018672077708074653, + "rotated_mean_y": -0.117278810001552, + "rotated_sigma_x": 0.3380214942514133, + "rotated_sigma_y": 0.09052009257286839, + "rotation_angle_mrad": -41.97373846575658 + }, + "7985": { + "target_position": -4.3, + "rotated_mean_x": -0.04662182705458965, + "rotated_mean_y": -0.11805472171925811, + "rotated_sigma_x": 0.3386689821250091, + "rotated_sigma_y": 0.08942650301390162, + "rotation_angle_mrad": -40.63085709672157 + }, + "7986": { + "target_position": -4.3, + "rotated_mean_x": -0.025832564163778573, + "rotated_mean_y": -0.13689720958850102, + "rotated_sigma_x": 0.33820794020844736, + "rotated_sigma_y": 0.0904331441754824, + "rotation_angle_mrad": -46.48847207379007 + }, + "7987": { + "target_position": -4.3, + "rotated_mean_x": -0.017044698484978287, + "rotated_mean_y": -0.14802365052656583, + "rotated_sigma_x": 0.35314163708409735, + "rotated_sigma_y": 0.09120403130843951, + "rotation_angle_mrad": -42.58629246929173 + }, + "7988": { + "target_position": -4.3, + "rotated_mean_x": -0.03484207000275738, + "rotated_mean_y": -0.14892503623068828, + "rotated_sigma_x": 0.3558627660414487, + "rotated_sigma_y": 0.09129908222479527, + "rotation_angle_mrad": -44.98119655706769 + }, + "8025": { + "target_position": -4.3, + "rotated_mean_x": -0.07914768250362011, + "rotated_mean_y": -0.12250924725940532, + "rotated_sigma_x": 0.32145959314279815, + "rotated_sigma_y": 0.1082440931482268, + "rotation_angle_mrad": -33.215254284915076 + }, + "8026": { + "target_position": -4.3, + "rotated_mean_x": -0.06710499989511251, + "rotated_mean_y": -0.11803814817984208, + "rotated_sigma_x": 0.32119241431862045, + "rotated_sigma_y": 0.10735435522642882, + "rotation_angle_mrad": -31.251728333683808 + }, + "8027": { + "target_position": -4.3, + "rotated_mean_x": -0.08055040042280137, + "rotated_mean_y": -0.11072191426556628, + "rotated_sigma_x": 0.32680388188871967, + "rotated_sigma_y": 0.10668727818974595, + "rotation_angle_mrad": -33.642699216608904 + }, + "8028": { + "target_position": -4.3, + "rotated_mean_x": -0.08379250372522418, + "rotated_mean_y": -0.10709915617681622, + "rotated_sigma_x": 0.3267926811986611, + "rotated_sigma_y": 0.10655937942189041, + "rotation_angle_mrad": -35.362349338138486 + }, + "8029": { + "target_position": -4.3, + "rotated_mean_x": -0.10741328359668914, + "rotated_mean_y": -0.10459672746049889, + "rotated_sigma_x": 0.3189686542441795, + "rotated_sigma_y": 0.10547311376523541, + "rotation_angle_mrad": -45.726959750343404 + }, + "8030": { + "target_position": -4.3, + "rotated_mean_x": -0.07431367563955617, + "rotated_mean_y": -0.11837201909245432, + "rotated_sigma_x": 0.32029902894225415, + "rotated_sigma_y": 0.10600300725291996, + "rotation_angle_mrad": -34.468349977479974 + }, + "8031": { + "target_position": -4.3, + "rotated_mean_x": -0.07173763193816436, + "rotated_mean_y": -0.12324723660310753, + "rotated_sigma_x": 0.3295047853819846, + "rotated_sigma_y": 0.10439887004293082, + "rotation_angle_mrad": -28.79156901837602 + }, + "8039": { + "target_position": -4.3, + "rotated_mean_x": -0.08105404568337934, + "rotated_mean_y": -0.11113614436122443, + "rotated_sigma_x": 0.32612498027143666, + "rotated_sigma_y": 0.10828069419220399, + "rotation_angle_mrad": -24.35638512333936 + }, + "8040": { + "target_position": -4.3, + "rotated_mean_x": -0.05794422420583392, + "rotated_mean_y": -0.10803323058561984, + "rotated_sigma_x": 0.3268389853764586, + "rotated_sigma_y": 0.10853338143285175, + "rotation_angle_mrad": -25.300245011554587 + }, + "8041": { + "target_position": -4.3, + "rotated_mean_x": -0.07940426958345932, + "rotated_mean_y": -0.11099875264706761, + "rotated_sigma_x": 0.31878724277266035, + "rotated_sigma_y": 0.10884887612245851, + "rotation_angle_mrad": -27.394273301326393 + }, + "8043": { + "target_position": -4.3, + "rotated_mean_x": -0.07121578822880463, + "rotated_mean_y": -0.11144265560676153, + "rotated_sigma_x": 0.32958298135991065, + "rotated_sigma_y": 0.11076153428777236, + "rotation_angle_mrad": -22.990078701077028 + }, + "8044": { + "target_position": -4.3, + "rotated_mean_x": -0.0660867048140294, + "rotated_mean_y": -0.10432017755622264, + "rotated_sigma_x": 0.3300276794683384, + "rotated_sigma_y": 0.11228201615894902, + "rotation_angle_mrad": -31.662603935862823 + }, + "8045": { + "target_position": -4.3, + "rotated_mean_x": -0.06812455047475856, + "rotated_mean_y": -0.12281299392014637, + "rotated_sigma_x": 0.3232108549421421, + "rotated_sigma_y": 0.11414520332891173, + "rotation_angle_mrad": -30.557996134325986 + }, + "8046": { + "target_position": -4.3, + "rotated_mean_x": -0.07661378200005259, + "rotated_mean_y": -0.1096015693759903, + "rotated_sigma_x": 0.3237081244960617, + "rotated_sigma_y": 0.10707049543137508, + "rotation_angle_mrad": -31.95529877076208 + }, + "8047": { + "target_position": -4.3, + "rotated_mean_x": -0.11592283461603257, + "rotated_mean_y": -0.09925226976610646, + "rotated_sigma_x": 0.3295543460335013, + "rotated_sigma_y": 0.1091174726419095, + "rotation_angle_mrad": -16.498728711918552 + }, + "8048": { + "target_position": -4.3, + "rotated_mean_x": -0.05720457511251374, + "rotated_mean_y": -0.11984760958654236, + "rotated_sigma_x": 0.31733928705062125, + "rotated_sigma_y": 0.10327148072533696, + "rotation_angle_mrad": -23.408964775753223 + }, + "8049": { + "target_position": -4.3, + "rotated_mean_x": -0.09199675543782215, + "rotated_mean_y": -0.12004241767818058, + "rotated_sigma_x": 0.31815679192292534, + "rotated_sigma_y": 0.1058961892847695, + "rotation_angle_mrad": -18.41517457263274 + }, + "8051": { + "target_position": -4.3, + "rotated_mean_x": -0.07776062304510763, + "rotated_mean_y": -0.11656064968415268, + "rotated_sigma_x": 0.3252610850846826, + "rotated_sigma_y": 0.10278021409394952, + "rotation_angle_mrad": -29.952650179449435 + }, + "8055": { + "target_position": -4.3, + "rotated_mean_x": -0.04635380869657864, + "rotated_mean_y": -0.10907374815645625, + "rotated_sigma_x": 0.3133639249338768, + "rotated_sigma_y": 0.10399941962693295, + "rotation_angle_mrad": -31.21893286459669 + }, + "8057": { + "target_position": -4.3, + "rotated_mean_x": -0.08968645970235978, + "rotated_mean_y": -0.12056333173375548, + "rotated_sigma_x": 0.32133549810692824, + "rotated_sigma_y": 0.1086222220462036, + "rotation_angle_mrad": -15.921200250507402 + }, + "8058": { + "target_position": -4.3, + "rotated_mean_x": -0.09350419895829268, + "rotated_mean_y": -0.12383291095379685, + "rotated_sigma_x": 0.32532409394232625, + "rotated_sigma_y": 0.10682707141177604, + "rotation_angle_mrad": -21.101803142542558 + }, + "8059": { + "target_position": -4.3, + "rotated_mean_x": -0.07302681965444255, + "rotated_mean_y": -0.12021232277329685, + "rotated_sigma_x": 0.32561575143160515, + "rotated_sigma_y": 0.1072202081178368, + "rotation_angle_mrad": -31.66538031847714 + }, + "8072": { + "target_position": -4.3, + "rotated_mean_x": 0.03584178814558944, + "rotated_mean_y": -0.15102169347085143, + "rotated_sigma_x": 0.33727901371636165, + "rotated_sigma_y": 0.11382033832940353, + "rotation_angle_mrad": 11.817782089230104 + }, + "8073": { + "target_position": -4.3, + "rotated_mean_x": -0.07048742357723133, + "rotated_mean_y": -0.1417540058769301, + "rotated_sigma_x": 0.3229933818719775, + "rotated_sigma_y": 0.10711870942719766, + "rotation_angle_mrad": -23.42132610080223 + }, + "8074": { + "target_position": -4.3, + "rotated_mean_x": -0.08608891772026173, + "rotated_mean_y": -0.13466263297716413, + "rotated_sigma_x": 0.31901591629126735, + "rotated_sigma_y": 0.10461330731410175, + "rotation_angle_mrad": -35.806372270567664 + }, + "8075": { + "target_position": -4.3, + "rotated_mean_x": -0.08770044201984302, + "rotated_mean_y": -0.11301566608808554, + "rotated_sigma_x": 0.3286592688509847, + "rotated_sigma_y": 0.10209773158158562, + "rotation_angle_mrad": -26.590766762034228 + }, + "8085": { + "target_position": -4.3, + "rotated_mean_x": -0.03724898699324563, + "rotated_mean_y": -0.09269450736674105, + "rotated_sigma_x": 0.31547895734708825, + "rotated_sigma_y": 0.10125580385041444, + "rotation_angle_mrad": -29.734555522794263 + }, + "8086": { + "target_position": -4.3, + "rotated_mean_x": -0.004112766704551729, + "rotated_mean_y": -0.0953411639114083, + "rotated_sigma_x": 0.33683396688923195, + "rotated_sigma_y": 0.09903690269475146, + "rotation_angle_mrad": -17.495037533969864 + }, + "8087": { + "target_position": -4.3, + "rotated_mean_x": -0.026597998644875557, + "rotated_mean_y": -0.09802393759860437, + "rotated_sigma_x": 0.34645011392506186, + "rotated_sigma_y": 0.11051706560888684, + "rotation_angle_mrad": -5.652117744631256 + }, + "8088": { + "target_position": -4.3, + "rotated_mean_x": -0.061809582007132755, + "rotated_mean_y": -0.07789610567434799, + "rotated_sigma_x": 0.3144850355359677, + "rotated_sigma_y": 0.1001023455939108, + "rotation_angle_mrad": -38.99943294001127 + }, + "8090": { + "target_position": -4.3, + "rotated_mean_x": -0.10803042190538667, + "rotated_mean_y": -0.07163180034305795, + "rotated_sigma_x": 0.3400015390963398, + "rotated_sigma_y": 0.10333046630601485, + "rotation_angle_mrad": -25.256258412965057 + }, + "8092": { + "target_position": -4.3, + "rotated_mean_x": -0.07417207043766089, + "rotated_mean_y": -0.07097249118686298, + "rotated_sigma_x": 0.33661454189584905, + "rotated_sigma_y": 0.10188613579741329, + "rotation_angle_mrad": -24.18392449240955 + }, + "8094": { + "target_position": -4.3, + "rotated_mean_x": -0.10424509356183052, + "rotated_mean_y": -0.07060455889901504, + "rotated_sigma_x": 0.334415533504094, + "rotated_sigma_y": 0.10070867061904525, + "rotation_angle_mrad": -21.542677827579375 + }, + "8095": { + "target_position": -4.3, + "rotated_mean_x": -0.1040980010323121, + "rotated_mean_y": -0.06909340699121978, + "rotated_sigma_x": 0.32482777089828757, + "rotated_sigma_y": 0.10196256736497583, + "rotation_angle_mrad": -35.722397548782865 + }, + "8096": { + "target_position": -4.3, + "rotated_mean_x": -0.08198016841687912, + "rotated_mean_y": -0.07152659058096561, + "rotated_sigma_x": 0.33832920932808114, + "rotated_sigma_y": 0.0992559505647756, + "rotation_angle_mrad": -8.165615761212278 + }, + "8098": { + "target_position": -4.3, + "rotated_mean_x": 0.023117572983090137, + "rotated_mean_y": -0.0948157871278651, + "rotated_sigma_x": 0.3277261218033586, + "rotated_sigma_y": 0.10108165645201359, + "rotation_angle_mrad": -15.458105438694497 + }, + "8099": { + "target_position": -4.3, + "rotated_mean_x": -0.07611336747642825, + "rotated_mean_y": -0.07629179626918209, + "rotated_sigma_x": 0.3400031965885583, + "rotated_sigma_y": 0.10336626796936195, + "rotation_angle_mrad": -8.46874684718849 + } +} diff --git a/analysis/data/v0_projection_2016_mc_7800_config.json b/analysis/data/v0_projection_2016_mc_7800_config.json new file mode 100644 index 000000000..cd6f5875c --- /dev/null +++ b/analysis/data/v0_projection_2016_mc_7800_config.json @@ -0,0 +1,10 @@ +{ + "7800": { + "target_position": -4.3, + "rotated_mean_x": -2.31164e-01, + "rotated_mean_y": -5.15632e-02, + "rotated_sigma_x": 2.01293e-01, + "rotated_sigma_y": 8.16385e-02, + "rotation_angle_mrad": -120.719 + } +} diff --git a/analysis/data/v0_projection_2016_mc_signal_config.json b/analysis/data/v0_projection_2016_mc_signal_config.json new file mode 100644 index 000000000..32f1a431a --- /dev/null +++ b/analysis/data/v0_projection_2016_mc_signal_config.json @@ -0,0 +1,10 @@ +{ + "7800": { + "target_position": -4.3, + "rotated_mean_x": 0.00032589410352742634, + "rotated_mean_y": 0.0005273697881972964, + "rotated_sigma_x": 0.19701452441142028, + "rotated_sigma_y": 0.09424252765235556, + "rotation_angle_mrad": -171.21036357655905 + } +} diff --git a/analysis/include/HistoManager.h b/analysis/include/HistoManager.h index 612b5775c..68badab03 100644 --- a/analysis/include/HistoManager.h +++ b/analysis/include/HistoManager.h @@ -182,7 +182,7 @@ class HistoManager { std::string xtitle, int nbinsX, float xmin, float xmax, std::string ytitle, int nbinsY, float ymin, float ymax, std::string ztitle, int nbinsZ, float zmin, float zmax); - + /** * @brief description * @@ -257,6 +257,20 @@ class HistoManager { */ void Fill2DHisto(const std::string& histoName, float valuex, float valuey, float weight=1.); + + /** + * @brief description + * + * @param histoName + * @param valuex + * @param valuey + * @param valuez + * @param weight + */ + void Fill3DHisto(const std::string& histoName, float valuex, float valuey, float valuez, float weight=1.); + + + /** * @brief Get histograms from input file * diff --git a/analysis/include/MutableTTree.h b/analysis/include/MutableTTree.h index f988850a1..bf2a3775a 100644 --- a/analysis/include/MutableTTree.h +++ b/analysis/include/MutableTTree.h @@ -87,6 +87,10 @@ class MutableTTree { */ bool variableExists(std::string variable); + virtual void addVariable(std::string variableName, double param)=0; + + void addVariableToTBranch(const std::string& variableName); + ~MutableTTree(); protected: diff --git a/analysis/include/SimpAnaTTree.h b/analysis/include/SimpAnaTTree.h index 92f8c3525..c213ace8a 100644 --- a/analysis/include/SimpAnaTTree.h +++ b/analysis/include/SimpAnaTTree.h @@ -64,9 +64,20 @@ class SimpAnaTTree : public MutableTTree { void addVariable_unc_vtx_abs_delta_z0tanlambda(); + //V0 target projection + //void vertex_target_projection(double target_pos); + //void vertex_target_projection_rotation(double angle); + //void addVariable_unc_vtx_proj_significance(); + + void addVariable(std::string variableName, double param) override; + + //Restructring add variables + void unc_vtx_ele_zalpha(double slope); + void unc_vtx_pos_zalpha(double slope); + void unc_vtx_deltaZ(); + //misc bool impactParameterCut2016Canonical(double mass); - bool testImpactParameterCut(); private: double skipCutVarValue_ = -9876543210.0;//=10" + }, + "posN2Dhits_gt" : { + "cut" : 10, + "id" : 2, + "info" : "p^{+} Track n2d Hits >=10" + }, "pSum_gt" : { "cut" : 0.5, - "id" : 1, + "id" : 3, "info" : "P_{e^{-}} + P_{e^{+}}>0.5 GeV" }, + "deltaZ_lt" : { + "cut" : 10.0, + "id" : 4, + "info" : "|#Delta Z| < 10" + }, "nVtxs_eq" : { "cut" : 1, - "id" : 2, + "id" : 5, "info" : "N_{vtx}=1" } - } diff --git a/analysis/selections/Tight_9hits_2021.json b/analysis/selections/Tight_9hits_2021.json new file mode 100644 index 000000000..c02e4e535 --- /dev/null +++ b/analysis/selections/Tight_9hits_2021.json @@ -0,0 +1,27 @@ +{ + "chi2unc_lt" : { + "cut" : 20.0, + "id" : 0, + "info" : "#chi^{2}_{unc}<20" + }, + "eleN2Dhits_gt" : { + "cut" : 9, + "id" : 1, + "info" : "e^{-} Track n2d Hits >=9" + }, + "posN2Dhits_gt" : { + "cut" : 9, + "id" : 2, + "info" : "p^{+} Track n2d Hits >=9" + }, + "pSum_gt" : { + "cut" : 0.5, + "id" : 3, + "info" : "P_{e^{-}} + P_{e^{+}}>0.5 GeV" + }, + "nVtxs_eq" : { + "cut" : 1, + "id" : 4, + "info" : "N_{vtx}=1" + } +} diff --git a/analysis/selections/Tight_L1L1_2021.json b/analysis/selections/Tight_L1L1_2021.json new file mode 100644 index 000000000..d9036eb89 --- /dev/null +++ b/analysis/selections/Tight_L1L1_2021.json @@ -0,0 +1,47 @@ +{ + "uncVtxX_lt" : { + "cut" : 1.0, + "id" : 0, + "info" : "|vtx x| < 1.0 mm" + }, + "uncVtxY_lt" : { + "cut" : 0.15, + "id" : 1, + "info" : "|vtx y| < 0.15 mm" + }, + "eleN2Dhits_gt" : { + "cut" : 10, + "id" : 2, + "info" : "e^{-} Track n2d Hits >=10" + }, + "posN2Dhits_gt" : { + "cut" : 10, + "id" : 3, + "info" : "p^{+} Track n2d Hits >=10" + }, + "chi2unc_lt" : { + "cut" : 20.0, + "id" : 4, + "info" : "#chi^{2}_{unc}<20" + }, + "L1Requirement_eq" : { + "cut" : 1, + "id" : 5, + "info" : "L1L1" + }, + "pSum_gt" : { + "cut" : 3.0, + "id" : 6, + "info" : "P_{e^{-}} + P_{e^{+}}>3.0 GeV" + }, + "deltaZ_lt" : { + "cut" : 10.0, + "id" : 7, + "info" : "|#Delta Z| < 10" + }, + "nVtxs_eq" : { + "cut" : 1, + "id" : 8, + "info" : "N_{vtx}=1" + } +} diff --git a/analysis/selections/Tight_pBot_2021.json b/analysis/selections/Tight_pBot_2021.json index 96e5e9ef0..6aa5ffd9c 100644 --- a/analysis/selections/Tight_pBot_2021.json +++ b/analysis/selections/Tight_pBot_2021.json @@ -1,22 +1,32 @@ { + "eleN2Dhits_gt" : { + "cut" : 10, + "id" : 0, + "info" : "e^{-} Track n2d Hits >=10" + }, + "posN2Dhits_gt" : { + "cut" : 10, + "id" : 1, + "info" : "p^{+} Track n2d Hits >=10" + }, "chi2unc_lt" : { "cut" : 20.0, - "id" : 0, + "id" : 2, "info" : "#chi^{2}_{unc}<20" }, "pSum_gt" : { "cut" : 0.5, - "id" : 1, + "id" : 3, "info" : "P_{e^{-}} + P_{e^{+}}>0.5 GeV" }, "volPos_bot" : { "cut" : 0.0, - "id" : 2, + "id" : 4, "info" : "e^{+} in Bot" }, "nVtxs_eq" : { "cut" : 1, - "id" : 3, + "id" : 5, "info" : "N_{vtx}=1" } diff --git a/analysis/selections/Tight_pTop_2021.json b/analysis/selections/Tight_pTop_2021.json index 46cbc1619..5150d6bd8 100644 --- a/analysis/selections/Tight_pTop_2021.json +++ b/analysis/selections/Tight_pTop_2021.json @@ -1,22 +1,32 @@ { + "eleN2Dhits_gt" : { + "cut" : 10, + "id" : 0, + "info" : "e^{-} Track n2d Hits >=10" + }, + "posN2Dhits_gt" : { + "cut" : 10, + "id" : 1, + "info" : "p^{+} Track n2d Hits >=10" + }, "chi2unc_lt" : { "cut" : 20.0, - "id" : 0, + "id" : 2, "info" : "#chi^{2}_{unc}<20" }, "pSum_gt" : { "cut" : 0.5, - "id" : 1, + "id" : 3, "info" : "P_{e^{-}} + P_{e^{+}}>0.5 GeV" }, "volPos_top" : { "cut" : 0.0, - "id" : 2, + "id" : 4, "info" : "e^{+} in Top" }, "nVtxs_eq" : { "cut" : 1, - "id" : 3, + "id" : 5, "info" : "N_{vtx}=1" } diff --git a/analysis/selections/simps/Tight_2016_simp_reach_CR.json b/analysis/selections/simps/Tight_2016_simp_reach_CR.json index 7c8885659..17448ab80 100644 --- a/analysis/selections/simps/Tight_2016_simp_reach_CR.json +++ b/analysis/selections/simps/Tight_2016_simp_reach_CR.json @@ -4,19 +4,14 @@ "id" : 0, "info" : "L1L1" }, - "pSum_lt" : { - "cut" : 2.4, - "id" : 1, - "info" : "P_{e^{-}} + P_{e^{+}} < 2.4 [GeV]" - }, "pSum_gt" : { "cut" : 1.9, - "id" : 2, + "id" : 1, "info" : "P_{e^{-}} + P_{e^{+}} > 1.9 GeV" }, "nVtxs_eq" : { "cut" : 1, - "id" : 3, + "id" : 2, "info" : "N_{vtx}=1" } } diff --git a/analysis/selections/simps/Tight_L1L1_nvtx1.json b/analysis/selections/simps/Tight_L1L1_nvtx1.json new file mode 100644 index 000000000..b09da7493 --- /dev/null +++ b/analysis/selections/simps/Tight_L1L1_nvtx1.json @@ -0,0 +1,12 @@ +{ + "L1Requirement_eq": { + "cut": 1, + "id": 0, + "info": "L1L1" + }, + "nVtxs_eq": { + "cut": 1, + "id": 1, + "info": "N_{vtx}=1" + } +} diff --git a/analysis/selections/simps/Tight_nocuts.json b/analysis/selections/simps/Tight_nocuts.json new file mode 100644 index 000000000..10dfc6dfa --- /dev/null +++ b/analysis/selections/simps/Tight_nocuts.json @@ -0,0 +1,7 @@ +{ + "chi2unc_lt" : { + "cut" : 20, + "id" : 0, + "info" : "#chi^{2}_{unc}<20" + } +} diff --git a/analysis/selections/simps/radMatchTight_2016_L1L1_nvtx1.json b/analysis/selections/simps/radMatchTight_2016_L1L1_nvtx1.json new file mode 100644 index 000000000..5eeb424b5 --- /dev/null +++ b/analysis/selections/simps/radMatchTight_2016_L1L1_nvtx1.json @@ -0,0 +1,17 @@ +{ + "L1Requirement_eq" : { + "cut" : 1, + "id" : 0, + "info" : "L1L1" + }, + "isRadEle_eq" : { + "cut" : 1, + "id" : 1, + "info" : "isRadEle" + }, + "nVtxs_eq" : { + "cut" : 1, + "id" : 2, + "info" : "N_{vtx}=1" + } +} diff --git a/analysis/selections/simps/radMatchTight_2016_simp_reach_CR.json b/analysis/selections/simps/radMatchTight_2016_simp_reach_CR.json index c9acaa8eb..825bd792b 100644 --- a/analysis/selections/simps/radMatchTight_2016_simp_reach_CR.json +++ b/analysis/selections/simps/radMatchTight_2016_simp_reach_CR.json @@ -4,24 +4,19 @@ "id" : 0, "info" : "L1L1" }, - "pSum_lt" : { - "cut" : 2.4, - "id" : 1, - "info" : "P_{e^{-}} + P_{e^{+}} < 2.4 [GeV]" - }, "pSum_gt" : { "cut" : 1.9, - "id" : 2, + "id" : 1, "info" : "P_{e^{-}} + P_{e^{+}} > 1.9 GeV" }, "isRadEle_eq" : { "cut" : 1, - "id" : 3, + "id" : 2, "info" : "isRadEle" }, "nVtxs_eq" : { "cut" : 1, - "id" : 4, + "id" : 3, "info" : "N_{vtx}=1" } } diff --git a/analysis/selections/simps/vertexSelection_2016_simp_nocuts.json b/analysis/selections/simps/vertexSelection_2016_simp_nocuts.json new file mode 100644 index 000000000..e36e63d9a --- /dev/null +++ b/analysis/selections/simps/vertexSelection_2016_simp_nocuts.json @@ -0,0 +1,7 @@ +{ + "Pair1_eq" : { + "cut" : 1, + "id" : 0, + "info" : "Trigger Pair1" + } +} diff --git a/analysis/selections/simps/vertexSelection_2016_simp_preselection.json b/analysis/selections/simps/vertexSelection_2016_simp_preselection.json new file mode 100644 index 000000000..d35526e6e --- /dev/null +++ b/analysis/selections/simps/vertexSelection_2016_simp_preselection.json @@ -0,0 +1,67 @@ +{ + "Pair1_eq" : { + "cut" : 1, + "id" : 0, + "info" : "Trigger Pair1" + }, + "eleTrkTime_lt" : { + "cut" : 6, + "id" : 1, + "info" : "e^{-} trkTime<6ns" + }, + "posTrkTime_lt" : { + "cut" : 6, + "id" : 2, + "info" : "e^{+} trkTime<6ns" + }, + "eleposCluTimeDiff_lt" : { + "cut" : 1.45, + "id" : 3, + "info" : "#Delta_{t}(clu_{e^{-}},clu_{e^{+}})<1.45ns" + }, + "eleTrkCluTimeDiff_lt" : { + "cut" : 4, + "id" : 4, + "info" : "e^{-} #Delta_{t}(trk,clu)<4ns" + }, + "posTrkCluTimeDiff_lt" : { + "cut" : 4, + "id" : 5, + "info" : "e^{+} #Delta_{t}(trk,clu)<4ns" + }, + "eleTrkChi2Ndf_lt" : { + "cut" : 20, + "id" : 6, + "info" : "e^{-} Track #chi^{2}/Ndf<20" + }, + "posTrkChi2Ndf_lt" : { + "cut" : 20, + "id" : 7, + "info" : "e^{+} Track #chi^{2}/Ndf<20" + }, + "eleMom_lt" : { + "cut" : 1.75, + "id" : 8, + "info" : "p e^{-}<1.75GeV" + }, + "eleN2Dhits_gt" : { + "cut" : 7, + "id" : 9, + "info" : "eleN2Dhits > 7" + }, + "posN2Dhits_gt" : { + "cut" : 7, + "id" : 10, + "info" : "posN2Dhits > 7" + }, + "chi2unc_lt" : { + "cut" : 20, + "id" : 11, + "info" : "#chi^{2}_{unc}<20" + }, + "maxVtxMom_lt" : { + "cut" : 2.4, + "id" : 12, + "info" : "maxVtxMom<2.4 [GeV]" + } +} diff --git a/analysis/selections/simps/vertexSelection_2016_simp_reach.json b/analysis/selections/simps/vertexSelection_2016_simp_reach.json new file mode 100644 index 000000000..e627733f4 --- /dev/null +++ b/analysis/selections/simps/vertexSelection_2016_simp_reach.json @@ -0,0 +1,72 @@ +{ + "Pair1_eq" : { + "cut" : 1, + "id" : 0, + "info" : "Trigger Pair1" + }, + "eleTrkTime_lt" : { + "cut" : 10, + "id" : 1, + "info" : "e^{-} trkTime<10ns" + }, + "posTrkTime_lt" : { + "cut" : 10, + "id" : 2, + "info" : "e^{+} trkTime<10ns" + }, + "eleposCluTimeDiff_lt" : { + "cut" : 1.45, + "id" : 3, + "info" : "#Delta_{t}(clu_{e^{-}},clu_{e^{+}})<1.45ns" + }, + "eleTrkCluTimeDiff_lt" : { + "cut" : 10, + "id" : 4, + "info" : "e^{-} #Delta_{t}(trk,clu)<10ns" + }, + "posTrkCluTimeDiff_lt" : { + "cut" : 10, + "id" : 5, + "info" : "e^{+} #Delta_{t}(trk,clu)<10ns" + }, + "eleTrkChi2Ndf_lt" : { + "cut" : 20, + "id" : 6, + "info" : "e^{-} Track #chi^{2}/Ndf<20" + }, + "posTrkChi2Ndf_lt" : { + "cut" : 20, + "id" : 7, + "info" : "e^{+} Track #chi^{2}/Ndf<20" + }, + "eleMom_lt" : { + "cut" : 1.75, + "id" : 8, + "info" : "p e^{-}<1.75GeV" + }, + "eleMom_gt" : { + "cut" : 0.4, + "id" : 9, + "info" : "p_e^{-}>0.4 [GeV]" + }, + "posMom_gt" : { + "cut" : 0.4, + "id" : 10, + "info" : "p_e^{+}>0.4 [GeV]" + }, + "eleN2Dhits_gt" : { + "cut" : 9, + "id" : 11, + "info" : "eleN2Dhits > 9" + }, + "posN2Dhits_gt" : { + "cut" : 9, + "id" : 12, + "info" : "posN2Dhits > 9" + }, + "chi2unc_lt" : { + "cut" : 20, + "id" : 13, + "info" : "#chi^{2}_{unc}<20" + } +} diff --git a/analysis/selections/vertexSelection_2021.json b/analysis/selections/vertexSelection_2021.json index c338eac8e..285e9ca8e 100644 --- a/analysis/selections/vertexSelection_2021.json +++ b/analysis/selections/vertexSelection_2021.json @@ -1,8 +1,8 @@ { - "eleMom_lt" : { - "cut" : 2.8, - "id" : 0, - "info" : "p e^{-} < 2.8 GeV" + "eleTrkTime_lt" : { + "cut" : 6.0, + "id" : 0, + "info" : "|e^{-} Track Time| < 6.0 ns" }, "eleTrkChi2Ndf_lt" : { "cut" : 30, @@ -14,34 +14,29 @@ "id" : 2, "info" : "e^{+} Track #chi^{2}/ndf < 30" }, - "eleMom_gt" : { - "cut" : 0.2, + "eleMom_lt" : { + "cut" : 2.9, "id" : 3, - "info" : "p_e^{-} > 0.2 [GeV]" + "info" : "p e^{-} < 2.9 GeV" }, - "posMom_gt" : { - "cut" : 0.2, + "eleMom_gt" : { + "cut" : 1.2, "id" : 4, - "info" : "p_e^{+} > 0.2 [GeV]" - }, - "eleN2Dhits_gt" : { - "cut" : 10, - "id" : 5, - "info" : "e^{-} Track n2d Hits >=10" + "info" : "p_e^{-} > 1.2 [GeV]" }, - "posN2Dhits_gt" : { - "cut" : 10, - "id" : 6, - "info" : "p^{+} Track n2d Hits >=10" + "posMom_gt" : { + "cut" : 1.2, + "id" : 5, + "info" : "p_e^{+} > 1.2 [GeV]" }, "chi2unc_lt" : { "cut" : 20, - "id" : 7, + "id" : 6, "info" : "#chi^{2}_{unc} < 20" }, "maxVtxMom_lt" : { - "cut" : 5.6, - "id" : 8, - "info" : "p_{vtx} < 5.6 [GeV]" + "cut" : 4.0, + "id" : 7, + "info" : "p_{vtx} < 4.0 [GeV]" } } diff --git a/analysis/src/AnaHelpers.cxx b/analysis/src/AnaHelpers.cxx index 04392b09e..74f6dee13 100644 --- a/analysis/src/AnaHelpers.cxx +++ b/analysis/src/AnaHelpers.cxx @@ -30,29 +30,29 @@ void AnaHelpers::InnermostLayerCheck(Track* trk, bool& foundL1, bool& foundL2) { bool hasL1 = false; bool hasL2 = false; bool hasL3 = false; - for (int ihit=0; ihitgetSvtHits().GetEntries();++ihit) { - TrackerHit* hit3d = (TrackerHit*) trk->getSvtHits().At(ihit); + for (int ihit=0; ihitgetHitLayers().size();++ihit) { + int hit3d = trk->getHitLayers().at(ihit); if(isKF){ - if (hit3d->getLayer() == 0 ) { + if (hit3d == 0 ) { innerCount++; } - if (hit3d->getLayer() == 1) { + if (hit3d == 1) { innerCount++; } - if (hit3d->getLayer() == 2) { + if (hit3d == 2) { innerCount++; hasL2 = true; } - if (hit3d->getLayer() == 3) { + if (hit3d == 3) { innerCount++; hasL3 = true; } } else{ - if (hit3d->getLayer() == 0 ) { + if (hit3d == 0 ) { innerCount++; } - if (hit3d->getLayer() == 1) { + if (hit3d == 1) { innerCount++; hasL1 = true; } diff --git a/analysis/src/HistoManager.cxx b/analysis/src/HistoManager.cxx index 3cee944c4..8dc77ec8d 100644 --- a/analysis/src/HistoManager.cxx +++ b/analysis/src/HistoManager.cxx @@ -85,7 +85,23 @@ void HistoManager::DefineHistos(){ hist.value().at("xtitle"),hist.value().at("binsX"),hist.value().at("minX"),hist.value().at("maxX"), hist.value().at("ytitle"),hist.value().at("binsY"),hist.value().at("minY"),hist.value().at("maxY")); } - + + //3D histogram + else if (extension == "hhh") { + histos3d[h_name] = plot3D(h_name, + hist.value().at("xtitle"), + hist.value().at("binsX"), + hist.value().at("minX"), + hist.value().at("maxX"), + hist.value().at("ytitle"), + hist.value().at("binsY"), + hist.value().at("minY"), + hist.value().at("maxY"), + hist.value().at("ztitle"), + hist.value().at("binsZ"), + hist.value().at("minZ"), + hist.value().at("maxZ")); + } }//loop on config } @@ -325,6 +341,24 @@ void HistoManager::Fill2DHisto(const std::string& histoName,float valuex, float } } +void HistoManager::Fill3DHisto(const std::string& histoName,float valuex, float valuey, float valuez, float weight) { + if (histos3d[m_name+"_"+histoName]) + histos3d[m_name+"_"+histoName]->Fill(valuex,valuey,valuez, weight); + else { + printWarnings_++; + if (doPrintWarnings_) { + if (printWarnings_ < maxWarnings_) + std::cout<<"ERROR::Fill3DHisto Histogram not found! "< *mcParts, std::string analysis, float weight) { + if(mcParts == nullptr) + std::cout << "MCPARTS IS NULL" << std::endl; int nParts = mcParts->size(); Fill1DHisto("numMCparts_h", (float)nParts, weight); int nMuons = 0; diff --git a/analysis/src/MutableTTree.cxx b/analysis/src/MutableTTree.cxx index 1d80f3953..0e71edcf6 100644 --- a/analysis/src/MutableTTree.cxx +++ b/analysis/src/MutableTTree.cxx @@ -109,3 +109,11 @@ void MutableTTree::initializeFlatTuple(TTree* tree, std::mapSetBranchAddress(varname.c_str(),tuple_map[varname]); } } + +void MutableTTree::addVariableToTBranch(const std::string& variableName){ + double* variable = new double{999.9}; + tuple_[variableName] = variable; + newtree_->Branch(variableName.c_str(), tuple_[variableName], (variableName+"/D").c_str()); + new_variables_[variableName] = variable; +} + diff --git a/analysis/src/SimpAnaTTree.cxx b/analysis/src/SimpAnaTTree.cxx index 9cdc62c18..5a4f886cc 100644 --- a/analysis/src/SimpAnaTTree.cxx +++ b/analysis/src/SimpAnaTTree.cxx @@ -1,31 +1,5 @@ #include -bool SimpAnaTTree::testImpactParameterCut(){ - double ele_z0 = getValue("unc_vtx_ele_track_z0"); - double pos_z0 = getValue("unc_vtx_pos_track_z0"); - double Z = getValue("unc_vtx_z"); - bool passCut = true; - if(ele_z0 > 0.0){ - if(ele_z0 < 0.029816*(Z-3.471875)) - passCut = false; - } - else{ - if(ele_z0 > -0.029530*(Z-3.471875)) - passCut = false; - } - - if(pos_z0 > 0.0){ - if(pos_z0 < 0.029816*(Z-3.471875)) - passCut = false; - } - else{ - if(pos_z0 > -0.029530*(Z-3.471875)) - passCut = false; - } - - return passCut; -} - bool SimpAnaTTree::impactParameterCut2016Canonical(double mass){ mass = mass/1000.0; double ele_z0 = getValue("unc_vtx_ele_track_z0"); @@ -59,6 +33,75 @@ bool SimpAnaTTree::impactParameterCut2016Canonical(double mass){ return passCut; } +void SimpAnaTTree::addVariable(std::string variableName, double param){ + if(variableName == "unc_vtx_ele_zalpha") + unc_vtx_ele_zalpha(param); + if(variableName == "unc_vtx_pos_zalpha") + unc_vtx_pos_zalpha(param); + if(variableName == "unc_vtx_deltaZ") + unc_vtx_deltaZ(); + +} +/* +void SimpAnaTTree::addVariableToTBranch(const std::string& variableName){ + double* variable = new double{999.9}; + tuple_[variableName] = variable; + newtree_->Branch(variableName, tuple_[variableName], variableName+"/D"); + new_variables_[variableName] = variable; +} +*/ +void SimpAnaTTree::unc_vtx_ele_zalpha(double slope){ + std::cout << "[SimpAnaTTree]::Adding variable unc_vtx_ele_zalpha with param " << slope << std::endl; + std::string variableName = "unc_vtx_ele_zalpha"; + //Create TBranch for variable + addVariableToTBranch(variableName); + + double& ele_z0 = *tuple_["unc_vtx_ele_track_z0"]; + double& recon_z = *tuple_["unc_vtx_z"]; + + std::function calculate_ele_zalpha = [&, slope]()->double{ + if(ele_z0 > 0.0) + return (recon_z - (ele_z0/slope)); + else + return (recon_z - (ele_z0/-slope)); + }; + functions_[variableName] = calculate_ele_zalpha; +} + +void SimpAnaTTree::unc_vtx_pos_zalpha(double slope){ + std::cout << "[SimpAnaTTree]::Adding variable unc_vtx_pos_zalpha with param " << slope << std::endl; + std::string variableName = "unc_vtx_pos_zalpha"; + //Create TBranch for variable + addVariableToTBranch(variableName); + + double& pos_z0 = *tuple_["unc_vtx_pos_track_z0"]; + double& recon_z = *tuple_["unc_vtx_z"]; + + std::function calculate_pos_zalpha = [&, slope]()->double{ + if(pos_z0 > 0.0) + return (recon_z - (pos_z0/slope)); + else + return (recon_z - (pos_z0/-slope)); + }; + functions_[variableName] = calculate_pos_zalpha; +} + +void SimpAnaTTree::unc_vtx_deltaZ(){ + std::cout << "[SimpAnaTTree]::Adding variable unc_vtx_deltaZ" << std::endl; + std::string variableName = "unc_vtx_deltaZ"; + //Create TBranch for variable + addVariableToTBranch(variableName); + + double& pos_z0 = *tuple_["unc_vtx_pos_track_z0"]; + double& ele_z0 = *tuple_["unc_vtx_ele_track_z0"]; + double& ele_tanlambda = *tuple_["unc_vtx_ele_track_tanLambda"]; + double& pos_tanlambda = *tuple_["unc_vtx_pos_track_tanLambda"]; + + std::function calculate_unc_vtx_deltaZ = [&]()->double{ + return std::abs((pos_z0/pos_tanlambda) - (ele_z0/ele_tanlambda)); + }; + functions_[variableName] = calculate_unc_vtx_deltaZ; +} void SimpAnaTTree::addVariable_unc_vtx_ele_zalpha(double slope){ std::cout << "[SimpAnaTTree]::Adding variable unc_vtx_ele_zalpha with slope " << slope << std::endl; @@ -325,6 +368,88 @@ void SimpAnaTTree::addVariable_unc_vtx_abs_delta_z0tanlambda(){ functions_["unc_vtx_abs_delta_z0tanlambda"] = calculate_abs_delta_z0tanlambda; } +/* +void SimpAnaTTree::vertex_target_projection(double target_pos){ + double* proj_x = new double{9999.9}; + tuple_["unc_vtx_proj_x"] = proj_x; + newtree_->Branch("unc_vtx_proj_x", tuple_["unc_vtx_proj_x"], "unc_vtx_proj_x/D"); + new_variables_["unc_vtx_proj_x"] = proj_x; + std::function calculate_vtx_proj_x = [&, target_pos]()->double{ + double* vtx_x = tuple_["unc_vtx_x"]; + double* vtx_z = tuple_["unc_vtx_z"]; + double* pz = tuple_["unc_vtx_pz"]; + double* px = tuple_["unc_vtx_px"]; + double projx = (*vtx_x - (*vtx_z - target_pos))*( *px / *pz); + return projx; + }; + functions_["unc_vtx_proj_x"] = calculate_vtx_proj_x; + + double* proj_y = new double{9999.9}; + tuple_["unc_vtx_proj_y"] = proj_y; + newtree_->Branch("unc_vtx_proj_y", tuple_["unc_vtx_proj_y"], "unc_vtx_proj_y/D"); + new_variables_["unc_vtx_proj_y"] = proj_y; + std::function calculate_vtx_proj_y = [&, target_pos]()->double{ + double* vtx_y = tuple_["unc_vtx_y"]; + double* vtx_z = tuple_["unc_vtx_z"]; + double* pz = tuple_["unc_vtx_pz"]; + double* py = tuple_["unc_vtx_py"]; + double projy = (*vtx_y - (*vtx_z - target_pos))*( *py / *pz); + return projy; + }; + functions_["unc_vtx_proj_y"] = calculate_vtx_proj_y; +} + +void SimpAnaTTree::vertex_target_projection_rotation(double angle){ + double* rot_x = new double{9999.9}; + tuple_["unc_vtx_proj_x_rot"] = rot_x; + newtree_->Branch("unc_vtx_proj_x_rot", tuple_["unc_vtx_proj_x_rot"], "unc_vtx_proj_x_rot/D"); + new_variables_["unc_vtx_proj_x_rot"] = rot_x; + std::function rotate_vtx_proj_x = [&, angle]()->double{ + double* x = tuple_["unc_vtx_proj_x"]; + double* y = tuple_["unc_vtx_proj_y"]; + return *x*std::cos(angle) - *y*std::sin(angle); + }; + functions_["unc_vtx_proj_x_rot"] = rotate_vtx_proj_x; + + double* rot_y = new double{9999.9}; + tuple_["unc_vtx_proj_y_rot"] = rot_y; + newtree_->Branch("unc_vtx_proj_y_rot", tuple_["unc_vtx_proj_y_rot"], "unc_vtx_proj_y_rot/D"); + new_variables_["unc_vtx_proj_y_rot"] = rot_y; + std::function rotate_vtx_proj_y = [&, angle]()->double{ + double* x = tuple_["unc_vtx_proj_x"]; + double* y = tuple_["unc_vtx_proj_y"]; + return *x*std::sin(angle) + *y*std::cos(angle); + }; + functions_["unc_vtx_proj_y_rot"] = rotate_vtx_proj_y; +} + +void SimpAnaTTree::addVariable_unc_vtx_proj_significance(){ + + //Create target projection values + vertex_target_projection(-4.3); + //Rotate projections + double angle = -0.060; + vertex_target_projection_rotation(angle); + + double* radius = new double{0.0}; + tuple_["unc_vtx_proj_significance"] = radius; + newtree_->Branch("unc_vtx_proj_significance", tuple_["unc_vtx_proj_significance"], "unc_vtx_proj_significance/D"); + new_variables_["unc_vtx_proj_significance"] = radius; + std::function calculate_vtx_proj_significance = [&]()->double{ + double xmean = -0.24; + double ymean = -0.06; + double xsigma = 0.38 ; + double ysigma = 0.1 ; + double* x = tuple_["unc_vtx_proj_x_rot"]; + double* y = tuple_["unc_vtx_proj_y_rot"]; + double x_sig = (*x-xmean)/xsigma; + double y_sig = (*y-ymean)/ysigma; + return std::sqrt( std::pow(x_sig,2) + std::pow(y_sig, 2) ); + }; + functions_["unc_vtx_proj_significance"] = calculate_vtx_proj_significance; +} +*/ + /* void SimpAnaTTree::addVariable_unc_vtx_zbravosumAlpha(double slope){ std::cout << "[SimpAnaTTree]::Adding variable unc_vtx_zbravosumAlpha" << std::endl; diff --git a/analysis/src/SimpEquations.cxx b/analysis/src/SimpEquations.cxx index 5f14b9047..bdf632afd 100644 --- a/analysis/src/SimpEquations.cxx +++ b/analysis/src/SimpEquations.cxx @@ -11,7 +11,7 @@ SimpEquations::SimpEquations(int year){ std::cout << "mass ratio Ap:V_D = " << mass_ratio_Ap_to_Vd_ << std::endl; std::cout << "ratio of dark pion mass to dark pion decay constant: " << ratio_mPi_to_fPi_ << std::endl; std::cout << "Dark Vector decay lepton mass: " << m_l_ << std::endl; - std::cout << "Alpha Dark: " << alpha_D_ << std::endl; + std::cout << "Alpha Dark: " << alpha_dark_ << std::endl; } @@ -24,7 +24,7 @@ SimpEquations::SimpEquations(int year, const std::string paramsConfigFile){ std::cout << "mass ratio Ap:V_D = " << mass_ratio_Ap_to_Vd_ << std::endl; std::cout << "ratio of dark pion mass to dark pion decay constant: " << ratio_mPi_to_fPi_ << std::endl; std::cout << "Dark Vector decay lepton mass: " << m_l_ << std::endl; - std::cout << "Alpha Dark: " << alpha_D_ << std::endl; + std::cout << "Alpha Dark: " << alpha_dark_ << std::endl; } @@ -45,7 +45,7 @@ void SimpEquations::loadParametersConfig(const std::string paramsConfigFile){ mass_ratio_Ap_to_Pid_ = param.value().at("mass_ratio_Ap_to_Pid"); ratio_mPi_to_fPi_ = param.value().at("ratio_mPi_to_fPi"); m_l_ = param.value().at("lepton_mass"); - alpha_D_ = param.value().at("alpha_dark"); + alpha_dark_ = param.value().at("alpha_dark"); } } catch(...) @@ -55,30 +55,93 @@ void SimpEquations::loadParametersConfig(const std::string paramsConfigFile){ } } -double SimpEquations::rate_2pi(double m_Ap, double m_pi, double m_V, double alpha_D){ - double coeff = (2.0*alpha_D/3.0) * m_Ap; +double SimpEquations::rate_2pi(double m_Ap, double m_pi, double m_V, double alpha_dark){ + double coeff = (2.0*alpha_dark/3.0) * m_Ap; double pow1 = std::pow((1-(4*m_pi*m_pi/(m_Ap*m_Ap))),3/2.); double pow2 = std::pow(((m_V*m_V)/((m_Ap*m_Ap)-(m_V*m_V))),2); return coeff*pow1*pow2; } -double SimpEquations::rate_Vpi(double m_Ap, double m_pi, double m_V, double alpha_D, double f_pi, bool rho, bool phi){ - double x = m_pi/m_Ap; - double y = m_V/m_Ap; - double coeff = alpha_D*Tv(rho, phi)/(192.*std::pow(M_PI,4)); +double SimpEquations::rate_Vrho_pi(double m_Ap, double m_pi, double m_V, + double alpha_dark, double f_pi){ + double x = m_pi / m_Ap; + double y = m_V / m_Ap; + double Tv = 3.0/4.0; + double coeff = alpha_dark*Tv/(192.*std::pow(M_PI,4)); return coeff * std::pow((m_Ap/m_pi),2) * std::pow(m_V/m_pi,2) * std::pow((m_pi/f_pi),4) * m_Ap*std::pow(Beta(x,y),3./2.); } -double SimpEquations::br_Vpi(double m_Ap, double m_pi, double m_V, double alpha_D, double f_pi, bool rho, bool phi){ - double rate = rate_Vpi(m_Ap,m_pi,m_V,alpha_D,f_pi,rho,phi) + rate_2pi(m_Ap,m_pi,m_V,alpha_D); - if(2*m_V < m_Ap) rate = rate_Vpi(m_Ap, m_pi, m_V, alpha_D, f_pi, rho,phi) + rate_2pi(m_Ap,m_pi,m_V,alpha_D) + rate_2V(m_Ap,m_V,alpha_D); - return rate_Vpi(m_Ap,m_pi,m_V,alpha_D,f_pi,rho,phi)/rate; +double SimpEquations::rate_Vphi_pi(double m_Ap, double m_pi, double m_V, + double alpha_dark, double f_pi){ + double x = m_pi / m_Ap; + double y = m_V / m_Ap; + double Tv = 3.0/2.0; + double coeff = alpha_dark*Tv/(192.*std::pow(M_PI,4)); + return coeff * std::pow((m_Ap/m_pi),2) * std::pow(m_V/m_pi,2) * std::pow((m_pi/f_pi),4) * m_Ap*std::pow(Beta(x,y),3./2.); +} +double SimpEquations::rate_Vcharged_pi(double m_Ap, double m_pi, double m_V, + double alpha_dark, double f_pi){ + double x = m_pi / m_Ap; + double y = m_V / m_Ap; + double Tv = 18.0-((3.0/2.0) +(3.0/4.0)) ; + double coeff = alpha_dark*Tv/(192.*std::pow(M_PI,4)); + return coeff * std::pow((m_Ap/m_pi),2) * std::pow(m_V/m_pi,2) * std::pow((m_pi/f_pi),4) * m_Ap*std::pow(Beta(x,y),3./2.); +} + +double SimpEquations::br_2pi(double m_Ap, double m_pi, double m_V, double alpha_dark, + double f_pi){ + double total_rate = rate_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + rate_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + rate_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + rate_2pi(m_Ap, m_pi, m_V, alpha_dark); + if(m_Ap > 2.0*m_V) + total_rate + rate_2V(m_Ap, m_V, alpha_dark); + return rate_2pi(m_Ap, m_pi, m_V, alpha_dark)/total_rate; +} + +double SimpEquations::br_Vrho_pi(double m_Ap, double m_pi, double m_V, double alpha_dark, + double f_pi){ + double total_rate = rate_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + rate_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + rate_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + rate_2pi(m_Ap, m_pi, m_V, alpha_dark); + if(m_Ap > 2.0*m_V) + total_rate + rate_2V(m_Ap, m_V, alpha_dark); + return rate_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi)/total_rate; } -double SimpEquations::br_2V(double m_Ap,double m_pi,double m_V,double alpha_D,double f_pi,double rho,double phi){ - if(2*m_V >= m_Ap) return 0.; - double rate = rate_Vpi(m_Ap,m_pi,m_V,alpha_D,f_pi,rho,phi) + rate_2pi(m_Ap,m_pi,m_V,alpha_D) + rate_2V(m_Ap,m_V,alpha_D); - return rate_2V(m_Ap,m_V,alpha_D)/rate; +double SimpEquations::br_Vphi_pi(double m_Ap, double m_pi, double m_V, double alpha_dark, + double f_pi){ + double total_rate = rate_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + rate_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + rate_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + rate_2pi(m_Ap, m_pi, m_V, alpha_dark); + if(m_Ap > 2.0*m_V) + total_rate + rate_2V(m_Ap, m_V, alpha_dark); + return rate_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi)/total_rate; +} + +double SimpEquations::br_Vcharged_pi(double m_Ap, double m_pi, double m_V, + double alpha_dark, double f_pi){ + double total_rate = rate_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + rate_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + rate_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + rate_2pi(m_Ap, m_pi, m_V, alpha_dark); + if(m_Ap > 2.0*m_V) + total_rate + rate_2V(m_Ap, m_V, alpha_dark); + return rate_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark,f_pi)/total_rate; +} + +double SimpEquations::br_2V(double m_Ap,double m_pi,double m_V,double alpha_dark,double f_pi,double rho,double phi){ + double total_rate = rate_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + rate_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + rate_Vcharged_pi(m_Ap, m_pi, m_V, alpha_dark, f_pi) + + rate_2pi(m_Ap, m_pi, m_V, alpha_dark); + if(m_Ap > 2.0*m_V){ + total_rate + rate_2V(m_Ap, m_V, alpha_dark); + return 0.0; + } + return rate_2V(m_Ap,m_V,alpha_dark)/total_rate; } double SimpEquations::Tv(bool rho,bool phi){ @@ -91,9 +154,9 @@ double SimpEquations::Beta(double x,double y){ return (1+std::pow(y,2)-std::pow(x,2)-2*y)*(1+std::pow(y,2)-std::pow(x,2)+2*y); } -double SimpEquations::rate_2V(double m_Ap,double m_V,double alpha_D){ +double SimpEquations::rate_2V(double m_Ap,double m_V,double alpha_dark){ double r = m_V/m_Ap; - return alpha_D/6. * m_Ap * f(r); + return alpha_dark/6. * m_Ap * f(r); } double SimpEquations::f(double r){ @@ -102,9 +165,9 @@ double SimpEquations::f(double r){ return num/den * std::pow((1-4*std::pow(r,2)),0.5); } -double SimpEquations::rate_2l(double m_Ap,double m_pi,double m_V,double eps,double alpha_D,double f_pi,double m_l,bool rho){ +double SimpEquations::rate_2l(double m_Ap,double m_pi,double m_V,double eps,double alpha_dark,double f_pi,double m_l,bool rho){ double alpha = 1./137.0; - double coeff = 16*M_PI*alpha_D*alpha*std::pow(eps,2)*std::pow(f_pi,2)/(3*std::pow(m_V,2)); + double coeff = 16*M_PI*alpha_dark*alpha*std::pow(eps,2)*std::pow(f_pi,2)/(3*std::pow(m_V,2)); double term1 = std::pow((std::pow(m_V,2)/(std::pow(m_Ap,2) - std::pow(m_V,2))),2); double term2 = std::pow((1-(4*std::pow(m_l,2)/std::pow(m_V,2))),0.5); double term3 = 1+(2*std::pow(m_l,2)/std::pow(m_V,2)); @@ -113,10 +176,10 @@ double SimpEquations::rate_2l(double m_Ap,double m_pi,double m_V,double eps,doub return coeff * term1 * term2 * term3 * m_V * constant; } -double SimpEquations::getCtau(double m_Ap,double m_pi,double m_V,double eps,double alpha_D,double f_pi,double m_l,bool rho){ +double SimpEquations::getCtau(double m_Ap,double m_pi,double m_V,double eps,double alpha_dark,double f_pi,double m_l,bool rho){ double c = 3.00e11; //mm/s double hbar = 6.58e-22; //MeV*sec - double rate = rate_2l(m_Ap,m_pi,m_V,eps,alpha_D,f_pi,m_l,rho); //MeV + double rate = rate_2l(m_Ap,m_pi,m_V,eps,alpha_dark,f_pi,m_l,rho); //MeV double tau = hbar/rate; double ctau = c*tau; return ctau; @@ -127,8 +190,8 @@ double SimpEquations::gamma(double m_V,double E_V){ return gamma; } -//Allow passage of radFrac, radAcc, and dNdm from python config -double SimpEquations::expectedSignalCalculation(double m_V, double eps, bool rho, bool phi, +//Calculate expected signal by passing radFrac, radAcc, and dNdm +double SimpEquations::expectedSignalCalculation(double m_V, double eps, bool rho, double E_V, TEfficiency* effCalc_h, double dNdm, double radFrac, double radAcc, double target_pos, double zcut){ //Signal mass dependent SIMP parameters @@ -137,7 +200,7 @@ double SimpEquations::expectedSignalCalculation(double m_V, double eps, bool rho double f_pi = m_pi/ratio_mPi_to_fPi_; //Mass in MeV - double ctau = getCtau(m_Ap,m_pi, m_V,eps,alpha_D_,f_pi,m_l_,rho); + double ctau = getCtau(m_Ap,m_pi, m_V,eps,alpha_dark_,f_pi,m_l_,rho); double gcTau = ctau * gamma(m_V/1000.0, E_V); //E_V in GeV std::cout << "gcTau: " << gcTau << std::endl; std::cout << "radFrac: " << radFrac << std::endl; @@ -154,13 +217,17 @@ double SimpEquations::expectedSignalCalculation(double m_V, double eps, bool rho effCalc_h->GetTotalHistogram()->GetBinWidth(zbin); } - //Total A' Production Rate double apProduction = (3.*137/2.)*3.14159*(m_Ap*eps*eps*radFrac*dNdm)/radAcc; //A' -> V+Pi Branching Ratio - double br_VPi = br_Vpi(m_Ap, m_pi, m_V, alpha_D_, f_pi, rho, phi); + double br_VPi = 0.0; + if(rho) + br_VPi = br_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark_, f_pi); + else + br_VPi = br_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark_, f_pi); + std::cout << "Branching ratio is " << br_VPi << std::endl; //Vector to e+e- BR = 1 double br_V_ee = 1.0; @@ -170,7 +237,7 @@ double SimpEquations::expectedSignalCalculation(double m_V, double eps, bool rho return expSignal; } -double SimpEquations::expectedSignalCalculation(double m_V, double eps, bool rho, bool phi, +double SimpEquations::expectedSignalCalculation(double m_V, double eps, bool rho, double E_V, TEfficiency* effCalc_h, double target_pos, double zcut){ //Signal mass dependent SIMP parameters @@ -179,7 +246,7 @@ double SimpEquations::expectedSignalCalculation(double m_V, double eps, bool rho double f_pi = m_pi/ratio_mPi_to_fPi_; //Mass in MeV - double ctau = getCtau(m_Ap,m_pi, m_V,eps,alpha_D_,f_pi,m_l_,rho); + double ctau = getCtau(m_Ap,m_pi, m_V,eps,alpha_dark_,f_pi,m_l_,rho); double gcTau = ctau * gamma(m_V/1000.0, E_V); //E_V in GeV //Calculate the Efficiency Vertex (Displaced VD Acceptance) @@ -198,7 +265,11 @@ double SimpEquations::expectedSignalCalculation(double m_V, double eps, bool rho /radiativeAcceptance(m_Ap); //A' -> V+Pi Branching Ratio - double br_VPi = br_Vpi(m_Ap, m_pi, m_V, alpha_D_, f_pi, rho, phi); + double br_VPi; + if(rho) + br_VPi = br_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark_, f_pi); + else + br_VPi = br_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark_, f_pi); //Vector to e+e- BR = 1 double br_V_ee = 1.0; @@ -209,11 +280,11 @@ double SimpEquations::expectedSignalCalculation(double m_V, double eps, bool rho return expSignal; } -double SimpEquations::expectedSignalCalculation(double m_Ap, double m_pi, double m_V, double eps, double alpha_D, - double f_pi, double m_l, bool rho, bool phi, double E_V, TEfficiency* effCalc_h, double target_pos, double zcut){ +double SimpEquations::expectedSignalCalculation(double m_Ap, double m_pi, double m_V, double eps, double alpha_dark, + double f_pi, double m_l, bool rho, double E_V, TEfficiency* effCalc_h, double target_pos, double zcut){ //Mass in MeV - double ctau = getCtau(m_Ap,m_pi, m_V,eps,alpha_D,f_pi,m_l,rho); + double ctau = getCtau(m_Ap,m_pi, m_V,eps,alpha_dark,f_pi,m_l,rho); double gcTau = ctau * gamma(m_V/1000.0, E_V); //E_V in GeV //Calculate the Efficiency Vertex (Displaced VD Acceptance) @@ -232,7 +303,11 @@ double SimpEquations::expectedSignalCalculation(double m_Ap, double m_pi, double /radiativeAcceptance(m_Ap); //A' -> V+Pi Branching Ratio - double br_VPi = br_Vpi(m_Ap, m_pi, m_V, alpha_D, f_pi, rho, phi); + double br_VPi; + if(rho) + br_VPi = br_Vrho_pi(m_Ap, m_pi, m_V, alpha_dark_, f_pi); + else + br_VPi = br_Vphi_pi(m_Ap, m_pi, m_V, alpha_dark_, f_pi); //Vector to e+e- BR = 1 double br_V_ee = 1.0; diff --git a/analysis/src/TrackHistos.cxx b/analysis/src/TrackHistos.cxx index 2f9ccfd1e..323635904 100644 --- a/analysis/src/TrackHistos.cxx +++ b/analysis/src/TrackHistos.cxx @@ -122,8 +122,32 @@ void TrackHistos::Fill1DVertex(Vertex* vtx, Fill1DHisto("ele_clusE_h",eleClus.getEnergy(),weight); Fill1DHisto("pos_clusE_h",posClus.getEnergy(),weight); Fill1DHisto("ele_EoP_h",eleClus.getEnergy()/p_ele.P(),weight); + Fill1DHisto("ele_trkClusX_h",ele_trk->getPositionAtEcal()[0] - eleClus.getPosition()[0],weight); + Fill1DHisto("ele_trkClusY_h",ele_trk->getPositionAtEcal()[1] - eleClus.getPosition()[1],weight); + Fill1DHisto("pos_trkClusX_h",pos_trk->getPositionAtEcal()[0] - posClus.getPosition()[0],weight); + Fill1DHisto("pos_trkClusY_h",pos_trk->getPositionAtEcal()[1] - posClus.getPosition()[1],weight); Fill1DHisto("pos_EoP_h",posClus.getEnergy()/p_pos.P(),weight); Fill2DHisto("EoP_hh", eleClus.getEnergy()/p_ele.P(), posClus.getEnergy()/p_pos.P(),weight); + Fill2DHisto("ele_trkClusX_p_hh", ele_trk->getP(),ele_trk->getPositionAtEcal()[0] - eleClus.getPosition()[0], weight); + Fill2DHisto("pos_trkClusX_p_hh", pos_trk->getP(),pos_trk->getPositionAtEcal()[0] - posClus.getPosition()[0], weight); + Fill2DHisto("ele_trkClusX_phi_hh", ele_trk->getPhi(),ele_trk->getPositionAtEcal()[0] - eleClus.getPosition()[0], weight); + Fill2DHisto("pos_trkClusX_phi_hh", pos_trk->getPhi(),pos_trk->getPositionAtEcal()[0] - posClus.getPosition()[0], weight); + Fill2DHisto("ele_trkClusX_z0_hh", ele_trk->getZ0(),ele_trk->getPositionAtEcal()[0] - eleClus.getPosition()[0], weight); + Fill2DHisto("pos_trkClusX_z0_hh", pos_trk->getZ0(),pos_trk->getPositionAtEcal()[0] - posClus.getPosition()[0], weight); + Fill2DHisto("ele_trkClusX_d0_hh", ele_trk->getD0(),ele_trk->getPositionAtEcal()[0] - eleClus.getPosition()[0], weight); + Fill2DHisto("pos_trkClusX_d0_hh", pos_trk->getD0(),pos_trk->getPositionAtEcal()[0] - posClus.getPosition()[0], weight); + Fill2DHisto("ele_trkClusX_tanL_hh", ele_trk->getTanLambda(),ele_trk->getPositionAtEcal()[0] - eleClus.getPosition()[0], weight); + Fill2DHisto("pos_trkClusX_tanL_hh", pos_trk->getTanLambda(), pos_trk->getPositionAtEcal()[0] - posClus.getPosition()[0], weight); + Fill2DHisto("ele_trkClusX_clusX_hh", eleClus.getPosition()[0],ele_trk->getPositionAtEcal()[0] - eleClus.getPosition()[0], weight); + Fill2DHisto("pos_trkClusX_clusX_hh", posClus.getPosition()[0], pos_trk->getPositionAtEcal()[0] - posClus.getPosition()[0], weight); + for(int i = 0; i < ele_trk->getHitLayers().size(); i++) + { + Fill2DHisto("ele_trkClusX_hitLay_hh", (float) ele_trk->getHitLayers().at(i), ele_trk->getPositionAtEcal()[0] - eleClus.getPosition()[0], weight); + } + for(int i = 0; i < pos_trk->getHitLayers().size(); i++) + { + Fill2DHisto("pos_trkClusX_hitLay_hh", (float) pos_trk->getHitLayers().at(i), pos_trk->getPositionAtEcal()[0] - posClus.getPosition()[0], weight); + } //Compute some extra variables @@ -183,7 +207,7 @@ void TrackHistos::Fill2DTrack(Track* track, float weight, const std::string& trk double d0 = track->getD0(); double z0 = track->getZ0(); - Fill2DHisto(trkname+"tanlambda_vs_phi0_hh",track->getPhi(),track->getTanLambda(), weight); + //Fill2DHisto(trkname+"tanlambda_vs_phi0_hh",track->getPhi(),track->getTanLambda(), weight); Fill2DHisto(trkname+"d0_vs_p_hh",track->getP(),d0,weight); Fill2DHisto(trkname+"d0_vs_phi0_hh",track->getPhi(),d0,weight); Fill2DHisto(trkname+"d0_vs_tanlambda_hh",track->getTanLambda(),d0,weight); @@ -192,6 +216,11 @@ void TrackHistos::Fill2DTrack(Track* track, float weight, const std::string& trk Fill2DHisto(trkname+"phi0_vs_p_hh",track->getP(),track->getPhi(),weight); Fill2DHisto(trkname+"z0_vs_phi0_hh",track->getPhi(),z0,weight); Fill2DHisto(trkname+"z0_vs_tanlambda_hh",track->getTanLambda(),z0,weight); + + Fill2DHisto(trkname+"TanLambda_vs_Phi_hh" ,track->getPhi() , track->getTanLambda() ,weight); + Fill2DHisto(trkname+"p_vs_Phi_hh" ,track->getPhi() , track->getP() ,weight); + Fill2DHisto(trkname+"p_vs_TanLambda_hh" ,track->getTanLambda() , track->getP() ,weight); + } } @@ -213,6 +242,7 @@ void TrackHistos::Fill1DTrack(Track* track, float weight, const std::string& trk Fill1DHisto(trkname+"invpT_h" ,-1*charge/track->getPt(),weight); Fill1DHisto(trkname+"TanLambda_h",track->getTanLambda() ,weight); Fill1DHisto(trkname+"Z0_h" ,track->getZ0() ,weight); + Fill1DHisto(trkname+"Z0oTanLambda_h",track->getZ0()/track->getTanLambda() ,weight); Fill1DHisto(trkname+"time_h" ,track->getTrackTime() ,weight); Fill1DHisto(trkname+"chi2_h" ,track->getChi2() ,weight); Fill1DHisto(trkname+"chi2ndf_h" ,track->getChi2Ndf() ,weight); @@ -222,8 +252,8 @@ void TrackHistos::Fill1DTrack(Track* track, float weight, const std::string& trk Fill1DHisto(trkname+"track_ypos_h",track->getPosition().at(1) ,weight); Fill1DHisto(trkname+"track_zpos_h",track->getPosition().at(2) ,weight); Fill1DHisto(trkname+"xpos_at_ecal_h",track->getPositionAtEcal().at(0) ,weight); - Fill1DHisto(trkname+"xpos_at_ecal_h",track->getPositionAtEcal().at(1) ,weight); - Fill1DHisto(trkname+"xpos_at_ecal_h",track->getPositionAtEcal().at(2) ,weight); + Fill1DHisto(trkname+"ypos_at_ecal_h",track->getPositionAtEcal().at(1) ,weight); + Fill1DHisto(trkname+"zpos_at_ecal_h",track->getPositionAtEcal().at(2) ,weight); //Top vs Bot if(track->getTanLambda() > 0.0) @@ -238,16 +268,12 @@ void TrackHistos::Fill1DTrack(Track* track, float weight, const std::string& trk Fill1DHisto(trkname+"TanLambda_err_h",track->getTanLambdaErr() ,weight); Fill1DHisto(trkname+"Z0_err_h" ,track->getZ0Err() ,weight); - for (int ihit=0; ihitgetSvtHits().GetEntries();++ihit) + for (int ihit=0; ihitgetHitLayers().size();++ihit) { - TrackerHit* hit2d = (TrackerHit*) track->getSvtHits().At(ihit); - Fill1DHisto(trkname+"hit_lay_h",hit2d->getLayer() ,weight); + int hit2d = track->getHitLayers().at(ihit); + Fill1DHisto(trkname+"hit_lay_h",(float) hit2d ,weight); } - //Fill 2D histos - Fill2DHisto(trkname+"TanLambda_vs_Phi_hh" ,track->getPhi() , track->getTanLambda() ,weight); - Fill2DHisto(trkname+"p_vs_Phi_hh" ,track->getPhi() , track->getP() ,weight); - Fill2DHisto(trkname+"p_vs_TanLambda_hh" ,track->getTanLambda() , track->getP() ,weight); - + //All Tracks Fill1DHisto(trkname+"sharingHits_h",0,weight); if (track->getNShared() == 0) diff --git a/analysis/src/ZBiHistos.cxx b/analysis/src/ZBiHistos.cxx index 1591a4f95..f182dd19b 100644 --- a/analysis/src/ZBiHistos.cxx +++ b/analysis/src/ZBiHistos.cxx @@ -395,7 +395,7 @@ TF1* ZBiHistos::fitExponentialPlusConst(std::string histogramName, double start_ return fitFunc; }*/ -TF1* ZBiHistos::fitExponentialPlusConst(std::string histogramName, double start_nevents){ +TF1* ZBiHistos::fitExponentialPlusExp(std::string histogramName, double start_nevents){ //Get z vertex distribution corresponding to Test Cut std::string histoname = m_name+"_"+histogramName+"_h"; @@ -460,11 +460,26 @@ TF1* ZBiHistos::fitExponentialPlusConst(std::string histogramName, double start_ double best_chi2 = 99999.9; double param_2; - TF1* fitFunc = new TF1("fitFunc", " (x<[2])*([0]*exp([1]*x)) + (x>=[2])*([0]*exp([1]*[2]))", xmin, xmax); + double best_seed_3; + double best_seed_4; + //TF1* fitFunc = new TF1("fitFunc", " (x<[2])*([0]*exp([1]*x)) + (x>=[2])*([3]*exp([4]*x))", xmin, xmax); + TF1* fitFunc = new TF1("fitFunc", " (x<[2])*([0]*exp([1]*x)) + (x>=[2])*(-[3]*x)", xmin, xmax); fitFunc->FixParameter(0, seed_0); fitFunc->FixParameter(1, seed_1); for(int i=0; i < 1000; i++){ fitFunc->FixParameter(2, (double)i /10.0); + //TFitResultPtr intermed_fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc_seed, "QS", "", (double)i /10.0 ,xmax); + //double seed_3 = intermed_fitResult->Parameter(0); + //double seed_4 = intermed_fitResult->Parameter(1); + //fitFunc->SetParameter(3,seed_3); + //fitFunc->SetParameter(4,seed_4); + + fitFunc_seed->FixParameter(0, seed_0); + fitFunc_seed->FixParameter(1, seed_1); + double seed_3 = fitFunc_seed->Eval(histos1d[histoname]->GetBinCenter(histos1d[histoname]->FindLastBinAbove(0.0))); + //std::cout << "LOOK FUCKER: " << seed_3 << std::endl; + //fitFunc->SetParameter(3,seed_3); + fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QSIM", "", xmin,xmax); if(fitResult->Ndf() <= 0) @@ -475,12 +490,212 @@ TF1* ZBiHistos::fitExponentialPlusConst(std::string histogramName, double start_ if(fitResult->Chi2()/fitResult->Ndf() < best_chi2){ best_chi2 = fitResult->Chi2()/fitResult->Ndf(); param_2 = fitResult->Parameter(2); - //std::cout << "set param_2: " << param_2 << std::endl; + best_seed_3 = seed_3; + //best_seed_4 = seed_4; } } //delete ran; + fitFunc->FixParameter(0, seed_0); + fitFunc->FixParameter(1, seed_1); + fitFunc->FixParameter(2, param_2); + fitFunc->FixParameter(3, best_seed_3); + //fitFunc->FixParameter(4, best_seed_4); + fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QSIM", "", xmin,xmax); + fitFunc->Draw(); + + //return histogram to original range + histos1d[histoname]->GetXaxis()->SetRangeUser(originalXmin,originalXmax); + + delete fitFunc_seed; + return fitFunc; +} + +TF1* ZBiHistos::fitExponentialTail(std::string histogramName, double start_nevents){ + + //Get z vertex distribution corresponding to Test Cut + std::string histoname = m_name+"_"+histogramName+"_h"; + + //If histogram is empty, return nullptr + if(histos1d[histoname]->GetEntries() < 1){ + std::cout << "[ZBiHistos]::WARNING: Background Model is NULL: " << + histogramName << " distribution is empty and could not be fit!" << std::endl; + return nullptr; + } + double originalXmin = histos1d[histoname]->GetXaxis()->GetXmin(); + double originalXmax = histos1d[histoname]->GetXaxis()->GetXmax(); + + //Start fit where start_nevents are left in tail + int lastbin = histos1d[histoname]->FindLastBinAbove(0.0); + if(histos1d[histoname]->Integral() < start_nevents) + start_nevents = histos1d[histoname]->GetBinContent(histos1d[histoname]->GetMaximumBin()); + int firstbin = lastbin - 1; + double test_integral = 0.0; + while(test_integral < start_nevents && histos1d[histoname]->GetBinLowEdge(firstbin) > 0.0){ + test_integral = histos1d[histoname]->Integral(firstbin, lastbin); + firstbin = firstbin - 1; + } + + double xmin = histos1d[histoname]->GetBinLowEdge(firstbin); + double xmax = histos1d[histoname]->GetBinLowEdge(lastbin+1); + histos1d[histoname]->GetXaxis()->SetRangeUser(xmin,xmax); + + //Initially fit tail with exponential to seed params + double best_seed_chi2 = 9999999.9; + double seed_0; + double seed_1; + + TF1* fitFunc = new TF1("fitFunc", "[0]*exp([1]*x)", xmin, xmax); + TFitResultPtr fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QSIM", "", xmin,xmax); + seed_0 = fitResult->Parameter(0); + seed_1 = fitResult->Parameter(1); + if(fitResult->Ndf() > 0){ + best_seed_chi2 = fitResult->Chi2()/fitResult->Ndf(); + } + + TRandom3* ran = new TRandom3(); + ran->SetSeed(0); + + //best_seed_chi2 = 999999.9; + for(int i=0; i < 30; i++){ + fitFunc->SetParameters(std::abs(ran->Gaus((double)start_nevents,3.0)), -std::abs(ran->Uniform(0.01,1.0))); + fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QSIM", "", xmin,xmax); + + if(fitResult->Ndf() <= 0) + continue; + if(fitFunc->GetProb() < 0.001) + continue; + + if(fitResult->Chi2()/fitResult->Ndf() < best_seed_chi2){ + if(fitResult->Parameter(0) < 1000.0) + continue; + best_seed_chi2 = fitResult->Chi2()/fitResult->Ndf(); + seed_0 = fitResult->Parameter(0); + seed_1 = fitResult->Parameter(1); + } + } + + double param_0 = seed_0; + double param_1 = seed_1; + for(int i=0; i < 30; i++){ + fitFunc->SetParameters(ran->Gaus(seed_0,1.0), -std::abs(ran->Gaus(seed_1,1.0))); + fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QSIM", "", xmin,xmax); + if(fitResult->Ndf() <= 0) + continue; + if(fitFunc->GetProb() < 0.001) + continue; + + if(fitResult->Chi2()/fitResult->Ndf() < best_seed_chi2){ + if(fitResult->Parameter(0) < 1000.0) + continue; + best_seed_chi2 = fitResult->Chi2()/fitResult->Ndf(); + param_0 = fitResult->Parameter(0); + param_1 = fitResult->Parameter(1); + } + } + + fitFunc->FixParameter(0, param_0); + fitFunc->FixParameter(1, param_1); + fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QSIM", "", xmin,xmax); + fitFunc->Draw(); + + //return histogram to original range + histos1d[histoname]->GetXaxis()->SetRangeUser(originalXmin,originalXmax); + + delete ran; + return fitFunc; +} + + + +TF1* ZBiHistos::fitExponentialPlusConst(std::string histogramName, double start_nevents){ + + //Get z vertex distribution corresponding to Test Cut + std::string histoname = m_name+"_"+histogramName+"_h"; + + //If histogram is empty, return nullptr + if(histos1d[histoname]->GetEntries() < 1){ + std::cout << "[ZBiHistos]::WARNING: Background Model is NULL: " << + histogramName << " distribution is empty and could not be fit!" << std::endl; + return nullptr; + } + double originalXmin = histos1d[histoname]->GetXaxis()->GetXmin(); + double originalXmax = histos1d[histoname]->GetXaxis()->GetXmax(); + + //Start fit where start_nevents are left in tail + int lastbin = histos1d[histoname]->FindLastBinAbove(0.0); + if(histos1d[histoname]->Integral() < start_nevents) + start_nevents = histos1d[histoname]->GetBinContent(histos1d[histoname]->GetMaximumBin()); + int firstbin = lastbin - 1; + double test_integral = 0.0; + while(test_integral < start_nevents && histos1d[histoname]->GetBinLowEdge(firstbin) > 0.0){ + test_integral = histos1d[histoname]->Integral(firstbin, lastbin); + firstbin = firstbin - 1; + } + + double xmin = histos1d[histoname]->GetBinLowEdge(firstbin); + double xmax = histos1d[histoname]->GetBinLowEdge(lastbin+1); + histos1d[histoname]->GetXaxis()->SetRangeUser(xmin,xmax); + + //Initially fit tail with exponential to seed params + double best_seed_chi2 = 9999999.9; + double seed_0; + double seed_1; + + TF1* fitFunc_seed = new TF1("fitFunc_seed", "[0]*exp([1]*x)", xmin, xmax); + //TF1* fitFunc_seed = new TF1("fitFunc_seed", "expo", xmin, xmax); + TFitResultPtr fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc_seed, "QSIM", "", xmin,xmax); + seed_0 = fitResult->Parameter(0); + seed_1 = fitResult->Parameter(1); + if(fitResult->Ndf() > 0 && fitFunc_seed->GetProb() > 0.001){ + best_seed_chi2 = fitResult->Chi2()/fitResult->Ndf(); + } + + TRandom3* ran = new TRandom3(); + ran->SetSeed(0); + + //best_seed_chi2 = 999999.9; + for(int i=0; i < 30; i++){ + fitFunc_seed->SetParameters(std::abs(ran->Gaus((double)start_nevents,3.0)), -std::abs(ran->Uniform(0.01,1.0))); + fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc_seed, "QSIM", "", xmin,xmax); + + if(fitResult->Ndf() <= 0) + continue; + if(fitFunc_seed->GetProb() < 0.001) + continue; + + if(fitResult->Chi2()/fitResult->Ndf() < best_seed_chi2){ + if(fitResult->Parameter(0) < 1000.0) + continue; + best_seed_chi2 = fitResult->Chi2()/fitResult->Ndf(); + seed_0 = fitResult->Parameter(0); + seed_1 = fitResult->Parameter(1); + } + } + + double best_chi2 = 99999.9; + double param_2; + TF1* fitFunc = new TF1("fitFunc", " (x<[2])*([0]*exp([1]*x)) + (x>=[2])*([0]*exp([1]*[2]))", xmin, xmax); + fitFunc->FixParameter(0, seed_0); + fitFunc->FixParameter(1, seed_1); + for(int i=0; i < 1000; i++){ + fitFunc->FixParameter(2, (double)i /10.0); + fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QSIM", "", xmin,xmax); + + if(fitResult->Ndf() <= 0) + continue; + if(fitFunc->GetProb() < 0.001) + continue; + + if(fitResult->Chi2()/fitResult->Ndf() < best_chi2){ + best_chi2 = fitResult->Chi2()/fitResult->Ndf(); + param_2 = fitResult->Parameter(2); + } + } + + delete ran; + fitFunc->FixParameter(0, seed_0); fitFunc->FixParameter(1, seed_1); fitFunc->FixParameter(2, param_2); @@ -494,6 +709,8 @@ TF1* ZBiHistos::fitExponentialPlusConst(std::string histogramName, double start_ return fitFunc; } +/* + TF1* ZBiHistos::fitExponentialTail(std::string histogramName, double start_nevents){ //Background Model Fit Function //TF1* fitFunc = new TF1("fitfunc","[0]*exp([1]*x)",0.0,100.0); @@ -572,6 +789,7 @@ TF1* ZBiHistos::fitExponentialTail(std::string histogramName, double start_neven return fitFunc; } +*/ /* TF1* ZBiHistos::fitExponentialTail(std::string histogramName, double start_nevents){ //Background Model Fit Function diff --git a/event/include/Track.h b/event/include/Track.h index 2df502e1d..685bfebf2 100644 --- a/event/include/Track.h +++ b/event/include/Track.h @@ -95,7 +95,12 @@ class Track : public TObject { void setCov(const std::vector& cov) {cov_ = cov;} std::vector getCov() {return cov_;} + + std::vector getHitLayers() {return hit_layers_;} + void addHitLayer(int layer) {hit_layers_.push_back(layer);} + void addMcpHit(int layer, int mcpID) {mcp_hits_.push_back(std::make_pair(layer,mcpID));} + std::vector> getMcpHits() {return mcp_hits_;} double getD0Err () const {return sqrt(cov_[0]);} double getPhiErr () const {return sqrt(cov_[2]);} @@ -329,6 +334,11 @@ class Track : public TObject { */ bool isBottomTrack() const { return track_volume_ ? true : false; }; + /** + * Set the number of tracker hits associated with this track. + */ + void setTrackerHitCount(int nHits) { n_hits_ = nHits; }; + /** * @return Number of tracker hits associated with this track. */ @@ -353,7 +363,7 @@ class Track : public TObject { private: /** Reference to the 3D hits associated with this track. */ - TRefArray tracker_hits_{}; + TRefArray tracker_hits_{TRefArray{}}; /** Reference to the reconstructed particle associated with this track. */ TRef particle_; @@ -372,7 +382,12 @@ class Track : public TObject { /** Cov matrix */ std::vector cov_; - + + /** hit layers */ + std::vector hit_layers_; + + /** truth mcp hits */ + std::vector> mcp_hits_; /** The distance of closest approach to the reference point. */ double d0_{-999.}; diff --git a/plotUtils/iterativeFitting.py b/plotUtils/iterativeFitting.py new file mode 100644 index 000000000..957ddaecd --- /dev/null +++ b/plotUtils/iterativeFitting.py @@ -0,0 +1,419 @@ +from ROOT import * +from array import array + +def MakeFit(histoGram, fitType, markerColor,fitrange=[-2e5,2e5],sigmarange=2): + + #make sure the styles are integers + markerColor = int(markerColor) + + #no Fit + if fitType=="noFit": + return None + elif fitType=="singleGausIterative": + fit = singleGausIterative(histoGram, sigmarange) + + fit.SetLineColor(markerColor) + + return fit + + +def ProfileYwithIterativeGaussFit(hist, mu_graph, sigma_graph, num_bins,fitrange=[-2e5,2e5],error_scaling=1.): + + if (num_bins < 1): + return + + minEntries = 50 + fDebug = False + + num_bins_x = hist.GetXaxis().GetNbins() + mu_graph.Rebin(num_bins) + sigma_graph.Rebin(num_bins) + + + f = (num_bins_x / num_bins + 2) + errs_mu = [0. for x in range (int(num_bins_x / num_bins + 2))] + errs_sigma = [0. for x in range (int(num_bins_x / num_bins + 2))] + + min_sigma = 0. + max_sigma = 0. + min_mu = 0. + max_mu = 0. + + num_skipped = 0 + + current_proj = None + + for i in range(1,num_bins_x+1, num_bins): + index = int(i / num_bins) + if (num_bins == 1): + index -= 1 + + + current_proj = hist.ProjectionY(hist.GetName()+"_"+str(index),i,i+num_bins-1) + + mu = 0. + mu_err = 0. + sigma = 0. + sigma_err = 0. + fit = None + + if (current_proj.GetEntries() < minEntries): + continue + else: + fit = singleGausIterative(current_proj,2,fitrange) + + mu = fit.GetParameter(1) + mu_err = fit.GetParError(1) + + sigma = fit.GetParameter(2) + sigma_err = fit.GetParError(2) + + + singlePlots=False + if (singlePlots): + c = TCanvas() + c.cd() + current_proj.Draw("p") + fit.Draw("same") + c.SaveAs(current_proj.GetName() + ".pdf") + + + if (sigma > max_sigma or max_sigma == 0): + max_sigma = sigma + if (sigma < min_sigma or min_sigma == 0): + min_sigma = sigma + + if (mu > max_mu or max_mu == 0) : + max_mu = mu + if (mu < min_mu or min_mu == 0) : + min_mu = mu + + value_x = (hist.GetXaxis().GetBinLowEdge(i) + hist.GetXaxis().GetBinUpEdge(i+num_bins-1))/2. + + #Important!! Use Fill to increment the graph with each iteration, or SetBinContent to replace contents... + if (sigma_graph!=None): + sigma_graph.Fill(value_x, sigma); + + if (mu_graph!=None): + mu_graph.Fill(value_x, mu); + + errs_mu[index+1] = mu_err * error_scaling + errs_sigma[index+1 ] = sigma_err * error_scaling + + + a_errs_mu = array("d", errs_mu) + a_errs_sigma = array("d",errs_sigma) + if (sigma_graph != None): + sigma_graph.SetError(a_errs_sigma) + sigma_graph.GetYaxis().SetTitleOffset(1.5) + sigma_graph.GetYaxis().SetTitle(hist.GetYaxis().GetTitle()) + sigma_graph.GetXaxis().SetTitle(hist.GetXaxis().GetTitle()) + sigma_graph.SetTitle("") + + if (mu_graph != None): + mu_graph.SetError(a_errs_mu) + mu_graph.GetYaxis().SetTitleOffset(1.5) + mu_graph.GetYaxis().SetTitle(hist.GetYaxis().GetTitle()) + mu_graph.GetXaxis().SetTitle(hist.GetXaxis().GetTitle()) + mu_graph.SetTitle("") + + if (fDebug and num_skipped): + print("Number of skipped bins: ", num_skipped) + + +#This fit requires some initial parameters guess to be reliable +def CrystalBallFit(hist, initialParameters, fitrange=[]): + debug = False + + if len(fitrange) ==0: + fitrange.append(histo.GetXaxis().GetXmin()) + fitrange.append(histo.GetXaxis().GetXmax()) + + CrystalFit = TF1("CrystalFit","[0] * ROOT::Math::crystalball_function(x,[1],[2],[3],[4])",fitrange[0],fitrange[1]) + + #https://en.wikipedia.org/wiki/Crystal_Ball_function#/media/File:CrystalBallFunction.svg + CrystalFit.SetParameter(0,initialParameters[0]) + CrystalFit.SetParameter(1,1) #alpha - 1 is a reasonable choice + CrystalFit.SetParameter(2,3) #n - 3 is a reasonable choice + CrystalFit.SetParameter(3,initialParameters[2]) #sigma + CrystalFit.SetParameter(4,initialParameters[1]) #mu + + + hist.Fit("CrystalFit","ORQN","same") + + + #c = TCanvas(hist.GetName()) + #hist.Draw() + #CrystalFit.Draw("same") + #c.SaveAs("./"+hist.GetName()+".pdf") + + return CrystalFit + + +def singleGausIterative(hist, sigmaRange,rangeFit=[]): + debug = False + # first perform a single Gaus fit across full range of histogram or in a specified range + + min = hist.GetBinLowEdge(1) + max = (hist.GetBinLowEdge(hist.GetNbinsX()))+hist.GetBinWidth(hist.GetNbinsX()) + userRange = False + minUser = -2e5 + maxUser = 2e5 + + if (len(rangeFit) != 0): + userRange = True + minUser = rangeFit[0] + maxUser = rangeFit[1] + min = minUser + max = maxUser + + fitA = TF1("fitA", "gaus", min,max) + hist.Fit("fitA","ORQN","same") + fitAMean = fitA.GetParameter(1) + fitASig = fitA.GetParameter(2) + + # performs a second fit with range determined by first fit + max = fitAMean + (fitASig*sigmaRange) + min = fitAMean - (fitASig*sigmaRange) + if (userRange): + if (max > maxUser): + max = maxUser + if (min < minUser): + min = minUser + + fitB = TF1("fitB", "gaus", min,max) + hist.Fit("fitB","ORQN","same") + fitMean = fitB.GetParameter(1) + fitSig = fitB.GetParameter(2) + + newFitSig = 99999 + newFitMean = 99999 + i = 0 + max = fitMean + (fitSig*sigmaRange) + min = fitMean - (fitSig*sigmaRange) + if (userRange): + if (max > maxUser): + max = maxUser + if (min < minUser): + min = minUser + + + fit = TF1("fit", "gaus", min,max) + + while abs(fitSig - newFitSig) > 0.0005 or abs(fitMean - newFitMean) > 0.0005: + + if(i > 0): + fitMean = newFitMean + fitSig = newFitSig + #print("i = ",i," fitMean = ",fitMean," fitSig = ",fitSig) + max = fitMean + (fitSig*sigmaRange) + min = fitMean - (fitSig*sigmaRange) + if (userRange): + if (max > maxUser): + max = maxUser + if (min < minUser): + min = minUser + + fit.SetRange(min,max) + hist.Fit("fit","ORQN","same") + newFitMean = fit.GetParameter(1) + newFitSig = fit.GetParameter(2) + #print("i = ",i," newFitMean = ", newFitMean, " newFitSig = ",newFitSig) + if(i > 50): + if debug: + print("WARNING terminate iterative gaus fit because of convergence problems") + print("final mean = ", newFitMean, ", previous iter mean = ", fitMean) + print("final sigma = ", newFitSig, ", previous iter sigma = ", fitSig) + break + + i = i + 1 + + + + if debug: + print("Final i = ",i," finalFitMean = ", fit.GetParameter(1), " finalFitSig = ",fit.GetParameter(2)) + + fit.SetLineWidth(2) + + return fit + + +#transform = 0: do nothing +#transform = 1: flip x +#transform = 2: flip y +#transform = 3: flip y and bin content +#transform = 4: flip x and y and bin content +def profileZwithCustomFit(hist,num_bins,fitrange=[-2e5,2e5],fittype="gaus",transform=0,name="") : + + + if ( not hist): + print("ProfileZwithIterativeGaussFit(): No histogram supplied!") + return + + num_bins_x = hist.GetXaxis().GetNbins(); + num_bins_y = hist.GetYaxis().GetNbins(); + xmin = hist.GetXaxis().GetXmin() + xmax = hist.GetXaxis().GetXmax() + ymin = hist.GetYaxis().GetXmin() + ymax = hist.GetYaxis().GetXmax() + + mu_name = hist.GetName() + "_mu_profiled_"+fittype+name + mu_graph = TH2D(mu_name,mu_name,num_bins_x, xmin, xmax,num_bins_y, ymin,ymax) + sigma_name = hist.GetName() + "_sigma_profiled_"+fittype+name + sigma_graph = TH2D(sigma_name,sigma_name,num_bins_x, xmin, xmax,num_bins_y, ymin,ymax) + mu_err_name = hist.GetName() + "_mu_err_profiled_"+fittype+name + mu_err_graph = TH2D(mu_err_name,mu_err_name,num_bins_x, xmin, xmax,num_bins_y, ymin,ymax) + sigma_err_name = hist.GetName() + "_sigma_err_profiled_"+fittype+name + sigma_err_graph = TH2D(sigma_err_name,sigma_err_name,num_bins_x, xmin, xmax,num_bins_y, ymin,ymax) + + + + mu_err = None + sigma_err = None + + minEntries = 50; + fDebug = 0; + + num_not_converged = 0; + num_skipped = 0; + + max_sigma = 0; + min_sigma = 0; + + max_mu = 0; + min_mu = 0; + + current_proj = None + + for i in range(1,num_bins_x+(num_bins==1),num_bins): + for j in range(1,num_bins_y+(num_bins==1),num_bins): + + index = i/num_bins; + index_y = j/num_bins; + + current_proj = hist.ProjectionZ(hist.GetName()+"_"+str(index)+"_"+str(index_y),i,i+num_bins-1,j,j+num_bins-1) + current_proj.SetTitle(hist.GetName()+"_"+str(index)+"_"+str(index_y)); + + current_mu = -999 + current_err_mu = -999 + current_sigma = -999 + current_err_sigma = -999 + + + if current_proj.GetEntries() < minEntries : + + current_mu = 0; + current_sigma = 0; + current_err_mu = 1; + current_err_sigma = 1; + + if (fDebug): + print("WARNING: Only "+current_proj.GetEntries()+" entries in bin "+index+","+index_y+ " in histogram " +hist.GetName()) + num_skipped+=1 + + + + else: + + + fit = singleGausIterative(current_proj,2,fitrange) + current_norm = fit.GetParameter(0) + current_mu = fit.GetParameter(1) + current_err_mu = fit.GetParError(1) + current_sigma = fit.GetParameter(2) + current_err_sigma = fit.GetParError(2) + + + #At thit point one can run a crystalBall + + if fittype=="cb": + crystalFit = CrystalBallFit(current_proj, + [current_norm, current_mu, current_sigma], + fitrange = [current_mu - 5*current_sigma, current_mu + 5*current_sigma]) + + current_mu = crystalFit.GetParameter(4) + current_err_mu = crystalFit.GetParError(4) + current_sigma = crystalFit.GetParameter(3) + current_err_sigma = crystalFit.GetParError(3) + + + if (current_sigma > max_sigma or max_sigma == 0): + max_sigma = current_sigma + if (current_sigma < min_sigma or min_sigma == 0): + min_sigma = current_sigma + if (current_mu > max_mu or max_mu == 0): + max_mu = current_mu + if (current_mu < min_mu or min_mu == 0): + min_mu = current_mu + + x_coord = (hist.GetXaxis().GetBinLowEdge(i) + hist.GetXaxis().GetBinUpEdge(i+num_bins-1))/2; + y_coord = (hist.GetYaxis().GetBinLowEdge(j) + hist.GetYaxis().GetBinUpEdge(j+num_bins-1))/2; + + if (sigma_graph): + sigma_graph.Fill(x_coord,y_coord,current_sigma) + if (mu_graph): + + if transform == 0: + mu_graph.Fill(x_coord,y_coord,current_mu) + elif transform == 1: + mu_graph.Fill(-x_coord,y_coord,current_mu) + elif transform == 2: + mu_graph.Fill(x_coord,-y_coord,current_mu) + elif transform == 3: + mu_graph.Fill(x_coord,-y_coord,-current_mu) + elif transform == 4: + mu_graph.Fill(-x_coord,-y_coord,-current_mu) + + #should probably be replace bin content, not fill? + if (sigma_err_graph): + sigma_err_graph.Fill(x_coord,y_coord,current_err_sigma); + if (mu_err_graph): + mu_err_graph.Fill(x_coord,y_coord,current_err_mu); + + + if (mu_graph) : + mu_graph.GetXaxis().SetTitle(hist.GetXaxis().GetTitle()) + mu_graph.GetYaxis().SetTitle(hist.GetYaxis().GetTitle()) + mu_graph.GetYaxis().SetTitleOffset(1) + mu_graph.GetZaxis().SetTitle(hist.GetZaxis().GetTitle()) + mu_graph.GetZaxis().SetTitleOffset(1.2) + mu_graph.SetTitle( "" ) + + if (sigma_graph) : + sigma_graph.GetXaxis().SetTitle(hist.GetXaxis().GetTitle()) + sigma_graph.GetYaxis().SetTitle(hist.GetYaxis().GetTitle()) + sigma_graph.GetYaxis().SetTitleOffset(1) + sigma_graph.GetZaxis().SetTitle(hist.GetZaxis().GetTitle()) + sigma_graph.GetZaxis().SetTitleOffset(1.2) + sigma_graph.SetTitle( "" ) + + + + if (mu_err_graph) : + mu_err_graph.GetXaxis().SetTitle(hist.GetXaxis().GetTitle()) + mu_err_graph.GetYaxis().SetTitle(hist.GetYaxis().GetTitle()) + mu_err_graph.GetYaxis().SetTitleOffset(1) + mu_err_graph.GetZaxis().SetTitle("Error of fit #mu:" + hist.GetZaxis().GetTitle()) + mu_err_graph.GetZaxis().SetTitleOffset(1.2) + mu_err_graph.SetTitle(hist.GetTitle()) + + + if (sigma_err_graph): + sigma_err_graph.GetXaxis().SetTitle(hist.GetXaxis().GetTitle()) + sigma_err_graph.GetYaxis().SetTitle(hist.GetYaxis().GetTitle()) + sigma_err_graph.GetYaxis().SetTitleOffset(1) + sigma_err_graph.GetZaxis().SetTitle("Error of fit #sigma: "+hist.GetZaxis().GetTitle()) + sigma_err_graph.GetZaxis().SetTitleOffset(1.2) + sigma_err_graph.SetTitle(hist.GetTitle()) + + return mu_graph,sigma_graph,mu_err_graph,sigma_err_graph + + + + + + + + + + diff --git a/processors/config/anaKalSimpTuple_Simps_cfg.py b/processors/config/anaKalSimpTuple_Simps_cfg.py index 0de147f9c..6f529c6b7 100644 --- a/processors/config/anaKalSimpTuple_Simps_cfg.py +++ b/processors/config/anaKalSimpTuple_Simps_cfg.py @@ -50,7 +50,7 @@ recoana_kf.parameters["mcColl"] = "MCParticle" recoana_kf.parameters["hitColl"] = "SiClusters" recoana_kf.parameters["ecalColl"] = "RecoEcalClusters" -recoana_kf.parameters["vtxSelectionjson"] = os.environ['HPSTR_BASE']+"/analysis/selections/simps/vertexSelection_2016_simp_reach.json" +recoana_kf.parameters["vtxSelectionjson"] = os.environ['HPSTR_BASE']+"/analysis/selections/simps/vertexSelection_2016_simp_preselection.json" recoana_kf.parameters["histoCfg"] = os.environ['HPSTR_BASE']+"/analysis/plotconfigs/tracking/simps/vtxAnalysis_2016_simp_reach.json" recoana_kf.parameters["mcHistoCfg"] = os.environ['HPSTR_BASE']+'/analysis/plotconfigs/mc/basicMC.json' ##### @@ -61,7 +61,14 @@ recoana_kf.parameters["isRadPDG"] = options.isRadPDG recoana_kf.parameters["makeFlatTuple"] = options.makeFlatTuple #recoana_kf.parameters["beamPosCfg"] = options.beamPosCorr -recoana_kf.parameters["beamPosCfg"] = os.environ['HPSTR_BASE']+'/analysis/data/beamspot_positions_2016.json' +#recoana_kf.parameters["beamPosCfg"] = os.environ['HPSTR_BASE']+'/analysis/data/beamspot_positions_2016.json' +if not options.isData: + recoana_kf.parameters["eleTrackTimeBias"] = -2.2 #MC + recoana_kf.parameters["posTrackTimeBias"] = -2.2 #MC +else: + recoana_kf.parameters["eleTrackTimeBias"] = -1.5 + recoana_kf.parameters["posTrackTimeBias"] = -1.5 + CalTimeOffset = -999 @@ -83,34 +90,9 @@ RegionPath+'Tight_2016_simp_reach_SR.json', RegionPath+'radMatchTight_2016_simp_reach_CR.json', RegionPath+'radMatchTight_2016_simp_reach_SR.json', - RegionPath+'../Tight_loose.json'] + RegionPath+'Tight_nocuts.json', + RegionPath+'Tight_L1L1_nvtx1.json'] -#RecoHitAna -recoana_gbl.parameters = recoana_kf.parameters.copy() -recoana_gbl.parameters["anaName"] = "vtxana_gbl" -recoana_gbl.parameters["vtxColl"] = "UnconstrainedV0Vertices" -recoana_gbl.parameters["tsColl"] = "TSData" -recoana_gbl.parameters["hitColl"] = "RotatedHelicalOnTrackHits" -recoana_gbl.parameters["trkColl"] = "GBLTracks" -recoana_gbl.parameters["mcColl"] = "MCParticle" -recoana_gbl.parameters["ecalColl"] = "RecoEcalClusters" -recoana_gbl.parameters["vtxSelectionjson"] = os.environ['HPSTR_BASE']+"/analysis/selections/vertexSelection.json" -recoana_gbl.parameters["histoCfg"] = os.environ['HPSTR_BASE']+"/analysis/plotconfigs/tracking/vtxAnalysis.json" -recoana_gbl.parameters["mcHistoCfg"] = os.environ['HPSTR_BASE']+'/analysis/plotconfigs/mc/basicMC.json' -##### -recoana_gbl.parameters["beamE"] = base.beamE[str(options.year)] -recoana_gbl.parameters["isData"] = options.isData -recoana_gbl.parameters["analysis"] = options.analysis -recoana_gbl.parameters["debug"] = 0 -recoana_gbl.parameters["isRadPDG"] = options.isRadPDG -recoana_gbl.parameters["makeFlatTuple"] = options.makeFlatTuple -recoana_gbl.parameters["beamPosCfg"] = options.beamPosCorr - -recoana_gbl.parameters["CalTimeOffset"] = CalTimeOffset -#Region definitions -RegionPath = os.environ['HPSTR_BASE']+"/analysis/selections/" -recoana_gbl.parameters["regionDefinitions"] = [RegionPath+'Tight.json', - RegionPath+'radMatchTight.json'] #MCParticleAna mcana.parameters["debug"] = 0 mcana.parameters["anaName"] = "mcAna" @@ -119,19 +101,12 @@ mcana.parameters["ecalHitColl"] = "EcalHits" mcana.parameters["histCfg"] = os.environ['HPSTR_BASE']+'/analysis/plotconfigs/mc/basicMC.json' -# -# RegionPath+'ESumCR.json', -# RegionPath+'TightNoSharedL0.json', -# RegionPath+'TightNoShared.json', # Sequence which the processors will run. #p.sequence = [recoana_kf,recoana_gbl] if (options.tracking == "KF"): print("Run KalmanFullTracks analysis") p.sequence = [recoana_kf] # ,mcana] -elif (options.tracking == "GBL"): - print("Run GBL analysis") - p.sequence = [recoana_gbl] # ,mcana] else: print("ERROR::Need to specify which tracks KF or GBL") exit(1) diff --git a/processors/config/anaNewVtxTuple_cfg.py b/processors/config/anaNewVtxTuple_cfg.py new file mode 100644 index 000000000..c78600065 --- /dev/null +++ b/processors/config/anaNewVtxTuple_cfg.py @@ -0,0 +1,83 @@ +import HpstrConf +import sys +import os +import baseConfig as base + +base.parser.add_argument("-f", "--makeFlatTuple", type=int, dest="makeFlatTuple", help="Make True to make vertex ana flat tuple", metavar="makeFlatTuple", default=0) + +options = base.parser.parse_args() + + +# Use the input file to set the output file name +infile = options.inFilename +outfile = options.outFilename + +print('Input file: %s' % infile) +print('Output file: %s' % outfile) + +p = HpstrConf.Process() + +p.run_mode = 1 +p.skip_events = options.skip_events +p.max_events = options.nevents + +#p.max_events = 1000 + +# Library containing processors +p.add_library("libprocessors") + +############################### +# Processors # +############################### + +vtxana = HpstrConf.Processor('vtxana', 'NewVertexAnaProcessor') + +############################### +# Processor Configuration # +############################### +#Vertex Analysis +vtxana.parameters["debug"] = 0 +vtxana.parameters["anaName"] = "vtxana" +vtxana.parameters["tsColl"] = "TSBank" +vtxana.parameters["hitColl"] = "SiClustersOnTrackOnPartOnUVtx" +vtxana.parameters["vtxColl"] = "UnconstrainedV0Vertices_KF" +vtxana.parameters["mcColl"] = "MCParticle" +vtxana.parameters["analysis"] = "vertex" +vtxana.parameters["vtxSelectionjson"] = os.environ['HPSTR_BASE']+'/analysis/selections/vertexSelection_2021.json' +vtxana.parameters["mcHistoCfg"] = os.environ['HPSTR_BASE']+'/analysis/plotconfigs/mc/basicMC.json' +vtxana.parameters["histoCfg"] = os.environ['HPSTR_BASE']+"/analysis/plotconfigs/tracking/vtxAnalysis_2021.json" +vtxana.parameters["beamE"] = base.beamE[str(options.year)] +vtxana.parameters["isData"] = options.isData +vtxana.parameters["isRadPDG"] = 622 +vtxana.parameters["makeFlatTuple"] = options.makeFlatTuple + +CalTimeOffset = -999 + +if (options.isData == 1): + CalTimeOffset = 56. + print("Running on data file: Setting CalTimeOffset %d" % CalTimeOffset) + +elif (options.isData == 0): + CalTimeOffset = 43. + print("Running on MC file: Setting CalTimeOffset %d" % CalTimeOffset) +else: + print("Specify which type of ntuple you are running on: -t 1 [for Data] / -t 0 [for MC]") + + +vtxana.parameters["CalTimeOffset"] = CalTimeOffset + +#Region definitions + +RegionPath = os.environ['HPSTR_BASE']+"/analysis/selections/" +if (options.year == 2019): + vtxana.parameters["regionDefinitions"] = [RegionPath+'Tight_2019.json', RegionPath+'Tight_pTop_2019.json', RegionPath+'Tight_pBot_2019.json'] +if (options.year == 2021): + vtxana.parameters["regionDefinitions"] = [RegionPath+'Tight_2021.json', RegionPath+'Tight_pTop_2021.json', RegionPath+'Tight_pBot_2021.json'] + +# Sequence which the processors will run. +p.sequence = [vtxana] + +p.input_files = infile +p.output_files = [outfile] + +p.printProcess() diff --git a/processors/config/anaSimps_2016_cfg.py b/processors/config/anaSimps_2016_cfg.py new file mode 100644 index 000000000..a673267d5 --- /dev/null +++ b/processors/config/anaSimps_2016_cfg.py @@ -0,0 +1,113 @@ +import HpstrConf +import sys +import os +import baseConfig as base + +base.parser.add_argument("-f", "--makeFlatTuple", type=int, dest="makeFlatTuple", + help="Make True to make vertex ana flat tuple", metavar="makeFlatTuple", default=1) +base.parser.add_argument("-r", "--isRadPDG", type=int, dest="isRadPDG", + help="Set radiative trident PDG ID", metavar="isRadPDG", default=622) +options = base.parser.parse_args() + +# Use the input file to set the output file name +infile = options.inFilename +outfile = options.outFilename + +outfile = outfile.split(".root")[0]+".root" + +print('Input file: %s' % infile) +print('Output file: %s' % outfile) + +p = HpstrConf.Process() + +p.run_mode = 1 +#p.max_events = 1000 + +# Library containing processors +p.add_library("libprocessors") + +############################### +# Processors # +############################### +vtxana = HpstrConf.Processor('vtxana', 'NewVertexAnaProcessor') +mcana = HpstrConf.Processor('mcpartana', 'MCAnaProcessor') +############################### +# Processor Configuration # +############################### +#RecoHitAna +vtxana.parameters["debug"] = 0 +vtxana.parameters["anaName"] = "vtxana" +vtxana.parameters["tsColl"] = "TSBank" +vtxana.parameters["vtxColl"] = "UnconstrainedV0Vertices_KF" +vtxana.parameters["mcColl"] = "MCParticle" +vtxana.parameters["hitColl"] = "SiClustersOnTrackOnPartOnUVtx" +vtxana.parameters["analysis"] = "vertex" +vtxana.parameters["vtxSelectionjson"] = os.environ['HPSTR_BASE']+"/analysis/selections/simps/vertexSelection_2016_simp_preselection.json" +vtxana.parameters["histoCfg"] = os.environ['HPSTR_BASE']+"/analysis/plotconfigs/tracking/simps/vtxAnalysis_2016_simp_reach_light.json" +vtxana.parameters["mcHistoCfg"] = os.environ['HPSTR_BASE']+'/analysis/plotconfigs/mc/basicMC.json' +##### +vtxana.parameters["beamE"] = base.beamE[str(options.year)] +vtxana.parameters["isData"] = options.isData +vtxana.parameters["isRadPDG"] = options.isRadPDG +vtxana.parameters["makeFlatTuple"] = options.makeFlatTuple +vtxana.parameters["beamPosCfg"] = "" +if options.isData: + vtxana.parameters["v0ProjectionFitsCfg"] = os.environ['HPSTR_BASE']+'/analysis/data/v0_projection_2016_config.json' +else: + vtxana.parameters["v0ProjectionFitsCfg"] = os.environ['HPSTR_BASE']+'/analysis/data/v0_projection_2016_mc_7800_config.json' #For tritrig and wab mc + #vtxana.parameters["v0ProjectionFitsCfg"] = os.environ['HPSTR_BASE']+'/analysis/data/v0_projection_2016_mc_signal_config.json' #For signal (accidentally gen with bspt=(0,0) + +if options.isData: + vtxana.parameters["eleTrackTimeBias"] = -1.5 + vtxana.parameters["posTrackTimeBias"] = -1.5 +else: + vtxana.parameters["eleTrackTimeBias"] = -2.2 #MC + vtxana.parameters["posTrackTimeBias"] = -2.2 #MC + #vtxana.parameters["eleTrackTimeBias"] = -5.5 #MC For TTs new smearing samples...due to readout bug + #vtxana.parameters["posTrackTimeBias"] = -5.5 #MC For TTs new smearing samples...due to readout bug + + +CalTimeOffset = -999 + +if (options.isData == 1): + CalTimeOffset = 56. + print("Running on data file: Setting CalTimeOffset %d" % CalTimeOffset) + +elif (options.isData == 0): + CalTimeOffset = 43. + print("Running on MC file: Setting CalTimeOffset %d" % CalTimeOffset) +else: + print("Specify which type of ntuple you are running on: -t 1 [for Data] / -t 0 [for MC]") + +vtxana.parameters["CalTimeOffset"] = CalTimeOffset +#Region definitions +RegionPath = os.environ['HPSTR_BASE']+"/analysis/selections/simps/" +if (options.year == 2016): + RegionDefinitions = [RegionPath+'Tight_2016_simp_reach_CR.json', + RegionPath+'Tight_2016_simp_reach_SR.json', + RegionPath+'Tight_nocuts.json', + RegionPath+'Tight_L1L1_nvtx1.json'] + if(options.isData != 1): + RegionDefinitions.extend([RegionPath+'radMatchTight_2016_simp_reach_CR.json', + RegionPath+'radMatchTight_2016_simp_reach_SR.json']) + + vtxana.parameters["regionDefinitions"] = RegionDefinitions + +#MCParticleAna +mcana.parameters["debug"] = 0 +mcana.parameters["anaName"] = "mcAna" +mcana.parameters["partColl"] = "MCParticle" +mcana.parameters["trkrHitColl"] = "TrackerHits" +mcana.parameters["ecalHitColl"] = "EcalHits" +mcana.parameters["histCfg"] = os.environ['HPSTR_BASE']+'/analysis/plotconfigs/mc/basicMC.json' + +# Sequence which the processors will run. +p.sequence = [vtxana] # ,mcana] + +p.skip_events = options.skip_events +p.max_events = options.nevents + +p.input_files = infile +p.output_files = [outfile] + +p.printProcess() diff --git a/processors/config/baseConfig.py b/processors/config/baseConfig.py index 4fdaf9f0c..23a344a58 100644 --- a/processors/config/baseConfig.py +++ b/processors/config/baseConfig.py @@ -32,8 +32,10 @@ beamE["2016"] = 2.3 beamE["2019"] = 4.55 beamE["2021"] = 3.74 +beamE["2022"] = 1.92 bfield = {} bfield["2016"] = 0.52 bfield["2019"] = 1.034 bfield["2021"] = 0.85 +bfield["2022"] = 0.437 diff --git a/processors/config/noRawHitTuple_cfg.py b/processors/config/noRawHitTuple_cfg.py new file mode 100644 index 000000000..3902199ec --- /dev/null +++ b/processors/config/noRawHitTuple_cfg.py @@ -0,0 +1,89 @@ +import HpstrConf +import sys + +import baseConfig as base +from baseConfig import bfield + +base.parser.add_argument("-w", "--tracking", type=str, dest="tracking", + help="Which tracking to use to make plots", metavar="tracking", default="KF") +base.parser.add_argument("-s", "--truthHits", type=int, dest="truthHits", + help="Get svt truth hits: 1=yes", metavar="truthHits", default=0) +base.parser.add_argument("-r", "--rawHits", type=int, dest="rawHits", + help="Keep raw svt hits: 1=yes", metavar="rawHits", default=0) +base.parser.add_argument("-TS", "--trackstate", type=str, dest="trackstate", + help="Specify Track State | 'AtECal' or 'AtTarget'. Default is origin (AtIP)", metavar="trackstate", default="") + + +options = base.parser.parse_args() + +# Use the input file to set the output file name +lcio_file = options.inFilename +root_file = options.outFilename + +print('LCIO file: %s' % lcio_file) +print('Root file: %s' % root_file) + +p = HpstrConf.Process() + +# p.max_events = 1000 +p.run_mode = 0 +p.skip_events = options.skip_events +p.max_events = options.nevents + +# Library containing processors +p.add_library("libprocessors") + +############################### +# Processors # +############################### +header = HpstrConf.Processor('header', 'EventProcessor') +vtx = HpstrConf.Processor('vtx', 'VertexProcessor') +mcpart = HpstrConf.Processor('mcpart', 'MCParticleProcessor') + +############################### +# Processor Configuration # +############################### +# Event +header.parameters["debug"] = 0 +header.parameters["headCollRoot"] = "EventHeader" +header.parameters["trigCollLcio"] = "TriggerBank" +header.parameters["rfCollLcio"] = "RFHits" +header.parameters["vtpCollLcio"] = "VTPBank" +header.parameters["vtpCollRoot"] = "VTPBank" +header.parameters["tsCollLcio"] = "TSBank" +header.parameters["tsCollRoot"] = "TSBank" + +# Vertex +vtx.parameters["debug"] = 0 +vtx.parameters["vtxCollLcio"] = 'UnconstrainedV0Vertices_KF' +vtx.parameters["vtxCollRoot"] = 'UnconstrainedV0Vertices_KF' +vtx.parameters["partCollRoot"] = 'ParticlesOnUVertices_KF' +vtx.parameters["kinkRelCollLcio"] = '' +vtx.parameters["trkRelCollLcio"] = 'KFTrackDataRelations' +vtx.parameters["trkhitCollRoot"] = '' +vtx.parameters["hitFitsCollLcio"] = 'SVTFittedRawTrackerHits' +vtx.parameters["rawhitCollRoot"] = '' +vtx.parameters["trackStateLocation"] = options.trackstate +vtx.parameters["mcPartRelLcio"] = 'SVTTrueHitRelations' +if options.trackstate == "": + vtx.parameters["bfield"] = bfield[str(options.year)] + +# MCParticle +mcpart.parameters["debug"] = 0 +mcpart.parameters["mcPartCollLcio"] = 'MCParticle' +mcpart.parameters["mcPartCollRoot"] = 'MCParticle' + + +# Sequence which the processors will run. + +sequence = [header, vtx] + +if (not options.isData): + sequence.append(mcpart) + +p.sequence = sequence + +p.input_files = lcio_file +p.output_files = [root_file] + +p.printProcess() diff --git a/processors/config/recoTuple_cfg.py b/processors/config/recoTuple_cfg.py index cfc957bc0..40ba376e3 100644 --- a/processors/config/recoTuple_cfg.py +++ b/processors/config/recoTuple_cfg.py @@ -9,7 +9,7 @@ base.parser.add_argument("-s", "--truthHits", type=int, dest="truthHits", help="Get svt truth hits: 1=yes", metavar="truthHits", default=0) base.parser.add_argument("-r", "--rawHits", type=int, dest="rawHits", - help="Keep raw svt hits: 1=yes", metavar="rawHits", default=0) + help="Keep raw svt hits: 1=yes", metavar="rawHits", default=0) options = base.parser.parse_args() diff --git a/processors/config/simpTuple_2016_cfg.py b/processors/config/simpTuple_2016_cfg.py new file mode 100644 index 000000000..b1d500b37 --- /dev/null +++ b/processors/config/simpTuple_2016_cfg.py @@ -0,0 +1,79 @@ +import HpstrConf +import sys +import baseConfig as base +from baseConfig import bfield + +base.parser.add_argument("-TS", "--trackstate", type=str, dest="trackstate", + help="Specify Track State | 'AtECal', 'AtTarget'. Default is origin ", metavar="trackstate", default="AtTarget") + +options = base.parser.parse_args() + +# Use the input file to set the output file name +lcio_file = options.inFilename +root_file = options.outFilename + +print('LCIO file: %s' % lcio_file) +print('Root file: %s' % root_file) + +p = HpstrConf.Process() + +# p.max_events = 1000 +p.run_mode = 0 +p.skip_events = options.skip_events +p.max_events = options.nevents + +# Library containing processors +p.add_library("libprocessors") + +############################### +# Processors # +############################### +header = HpstrConf.Processor('header', 'EventProcessor') +vtx = HpstrConf.Processor('vtx', 'VertexProcessor') +mcpart = HpstrConf.Processor('mcpart', 'MCParticleProcessor') + +############################### +# Processor Configuration # +############################### +# Event +header.parameters["debug"] = 0 +header.parameters["headCollRoot"] = "EventHeader" +header.parameters["trigCollLcio"] = "TriggerBank" +header.parameters["rfCollLcio"] = "RFHits" +header.parameters["vtpCollLcio"] = "VTPBank" +header.parameters["vtpCollRoot"] = "VTPBank" +header.parameters["tsCollLcio"] = "TSBank" +header.parameters["tsCollRoot"] = "TSBank" + +# Vertex +vtx.parameters["debug"] = 0 +vtx.parameters["vtxCollLcio"] = 'UnconstrainedV0Vertices_KF' +vtx.parameters["vtxCollRoot"] = 'UnconstrainedV0Vertices_KF' +vtx.parameters["partCollRoot"] = 'ParticlesOnUVertices_KF' +vtx.parameters["kinkRelCollLcio"] = '' +vtx.parameters["trkRelCollLcio"] = 'KFTrackDataRelations' +vtx.parameters["trkhitCollRoot"] = '' +vtx.parameters["hitFitsCollLcio"] = 'SVTFittedRawTrackerHits' +vtx.parameters["rawhitCollRoot"] = '' +vtx.parameters["trackStateLocation"] = options.trackstate +if options.trackstate == "": + vtx.parameters["bfield"] = bfield[str(options.year)] +vtx.parameters["mcPartRelLcio"] = 'SVTTrueHitRelations' + +# MCParticle +mcpart.parameters["debug"] = 0 +mcpart.parameters["mcPartCollLcio"] = 'MCParticle' +mcpart.parameters["mcPartCollRoot"] = 'MCParticle' + +sequence = [header, vtx] + +# If MC, get MCParticles +if (not options.isData): + sequence.append(mcpart) + +p.sequence = sequence + +p.input_files = lcio_file +p.output_files = [root_file] + +p.printProcess() diff --git a/processors/config/simp_zbi_config.py b/processors/config/simp_zbi_config.py index 10bbbbfa6..c94694795 100644 --- a/processors/config/simp_zbi_config.py +++ b/processors/config/simp_zbi_config.py @@ -20,7 +20,7 @@ help="Define Zcut based on n background events in background fit", metavar="ztail_nevents", default=0.5) base.parser.add_argument("-z", "--scan_zcut", type=int, dest="scan_zcut", - help="Choose best ZBi using Zcut Scan (1=yes, 0 = No)", metavar="scan_zcut", default=1) + help="Choose best ZBi using Zcut Scan (1=yes, 0 = No)", metavar="scan_zcut", default=0) base.parser.add_argument("-m", "--step_size", type=float, dest="step_size", help="Cut % of signal with each iteration", metavar="step_size", default=0.01) @@ -30,6 +30,8 @@ options = base.parser.parse_args() p = HpstrConf.Process() + +######################################### Expected Signal Calculation ##################################################### #Calculated in 'makeRadFrac.py' def radiativeFraction(mass): radF = -1.04206e-01 + 9.92547e-03*mass + -1.99437e-04*pow(mass,2) + 1.83534e-06*pow(mass,3) + -7.93138e-9*pow(mass,4) + 1.30456e-11*pow(mass,5) #alic 2016 simps kf 11/15/22 @@ -40,7 +42,13 @@ def radiativeAcceptance(mass): acc = ( -7.35934e-01 + 9.75402e-02*mass + -5.22599e-03*pow(mass,2) + 1.47226e-04*pow(mass,3) + -2.41435e-06*pow(mass,4) + 2.45015e-08*pow(mass,5) + -1.56938e-10*pow(mass,6) + 6.19494e-13*pow(mass,7) + -1.37780e-15*pow(mass,8) + 1.32155e-18*pow(mass,9) ) #alic 2016 simps kf 11/15/22 return acc -#lcio_file = options.inFilename[0] +############################################################################################################### + +def chooseIterativeCutVariables(zbi, cut_variables = [], new_variables = [], variable_params = []): + zbi.parameters['cutVariables'] = cut_variables + zbi.parameters['add_new_variaibles'] = new_variables + zbi.parameters['new_variable_params'] = variable_params + in_file = 'dummy.root' out_file = options.outFilename @@ -53,68 +61,57 @@ def radiativeAcceptance(mass): ############################### # Processors # ################################ -#Get expected signal calculation values -#dNdm_VdMass = {55.0:602.e3 , 60.0:324.e3 , 75.0:93200 , 90.0:26247 , 110.0:4133} -#dNdm = dNdm_VdMass[simp_vd_mass] -#logeps2_VdMass = {55.0:-6.3 , 60.0:-6.4 , 75.0:-6.6, 90.0:-6.7 , 110.0:-6.8} -logeps2 = -6.5 +#Initialize processor zbi = HpstrConf.Processor('zbi','SimpZBiOptimizationProcessor') -#basic config -zbi.parameters['max_iteration'] = 5 + +#Configure basic settings +zbi.parameters['max_iteration'] = 25 zbi.parameters['year'] = 2016 zbi.parameters['debug'] = 1 -#zbi.parameters['debug'] = options.debug zbi.parameters['outFileName'] = options.outFilename zbi.parameters['scan_zcut'] = options.scan_zcut #1 will calculate ZBi as function of zcut position zbi.parameters['step_size'] = options.step_size #Specify %variable in signal to cut with each iteration zbi.parameters['ztail_events'] = options.ztail_nevents # 0.5 is the minimum allowed. ZBi calc breaks if 0.0 -#json configs -zbi.parameters['variableHistCfgFilename'] = '/sdf/group/hps/users/alspellm/src/test/hpstr/analysis/plotconfigs/tracking/zbiCutVariables.json' -zbi.parameters['cuts_cfgFile'] = '/sdf/group/hps/users/alspellm/run/simp_zbi_optimization/delta_z0_tanlambda/iterativeCuts.json' -eq_cfg_file = '/sdf/group/hps/users/alspellm/run/simp_zbi_optimization/delta_z0_tanlambda/simp_parameters.json' +#Histogram Config +zbi.parameters['variableHistCfgFilename'] = '/sdf/group/hps/users/alspellm/src/test/hpstr/analysis/plotconfigs/tracking/zbi_optimization_histograms.json' +#Config SIMP model +eq_cfg_file = '/sdf/group/hps/users/alspellm/src/hpstr/analysis/selections/simps/simp_parameters.json' zbi.parameters['eq_cfgFile'] = eq_cfg_file - -#abs delta z0/tanlambda -zbi.parameters['cutVariables'] = ["unc_vtx_abs_delta_z0tanlambda"] -zbi.parameters['add_new_variables'] = ["unc_vtx_abs_delta_z0tanlambda"] -zbi.parameters['new_variable_params'] = [0.0] - -#special config -zbi.parameters['testSpecialCut'] = 0 - -#background config -zbi.parameters['bkgVtxAnaFilename'] = '/sdf/group/hps/users/alspellm/projects/THESIS/data/2016/hps_7800/20230726_recon/20230726_tuples/20230727_ana/ana_files/hadd_hps_7800_KF_24nsClusterWindow_ana_367_files.root' -zbi.parameters['bkgVtxAnaTreename'] = 'vtxana_kf_Tight_2016_simp_reach_SR' +#Choose initial set of cuts +zbi.parameters['cuts_cfgFile'] = '/sdf/group/hps/users/alspellm/projects/THESIS/analysis/tight_selection_studies/v0_projection/zbi_opt/iterativeCuts.json' +#Choose cut variables to tighten iteratively. Must be present in json file above^ +chooseIterativeCutVariables(zbi, ["unc_vtx_proj_sig"]) + +#Configure Background +zbi.parameters['bkgVtxAnaFilename'] = '/sdf/group/hps/users/alspellm/projects/THESIS/data/2016/BLPass4c_20231006/ana_20231019/full_hadd_blpass4c_ana.root' +zbi.parameters['bkgVtxAnaTreename'] = 'vtxana_Tight_2016_simp_reach_SR' zbi.parameters['background_sf'] = 10.0 -#mc signal config -#zbi.parameters['signal_sf'] = 1.0 +#Configure Signal zbi.parameters['signal_sf'] = 1.0 zbi.parameters['signal_mass'] = options.mass zbi.parameters['mass_window_nsigma'] = 2.8 -zbi.parameters['signalVtxAnaFilename'] = '/sdf/group/hps/users/alspellm/projects/THESIS/mc/2016/simps/signal_beam/20230713_slic/20230713_readout/20230713_recon/20230721_tuples/20230724_ana/ana_files/hadd_simp_%s_MeV_egsv6_HPS-PhysicsRun2016-Pass2_recon_targetTracks_NObeamPosCorr_ana.root'%(int(options.mass)) -zbi.parameters['signalVtxAnaTreename'] = 'vtxana_kf_radMatchTight_2016_simp_reach_SR' -zbi.parameters['signalVtxMCSelection'] = 'vtxana_kf_mc_radMatchTight_2016_simp_reach_SR' +zbi.parameters['signalVtxAnaFilename'] = '/sdf/group/hps/users/alspellm/projects/THESIS/mc/2016/simps/signal_beam/20230713_slic/20230713_readout/hps-java_v5pt2pt1/pass4/recon_20231009/ana_20231020/hadd_simp_signal_%s_MeV_beam_ana.root'%(int(options.mass)) +zbi.parameters['signalVtxAnaTreename'] = 'vtxana_radMatchTight_2016_simp_reach_SR' +zbi.parameters['signalVtxMCSelection'] = 'vtxana_mc_radMatchTight_2016_simp_reach_SR' zbi.parameters['signalMCAnaFilename'] = '/sdf/group/hps/users/alspellm/projects/THESIS/mc/2016/simps/slic/20230713_slic/20230724_slic_ana/ana_files/hadd_simp_%s_MeV_rot_slic_mcana.root'%(int(options.mass)) zbi.parameters['signal_pdgid'] = '625' zbi.parameters['logEps2'] = options.logeps2 -#Total A' config +#Configure Expected Signal Calculation #Read simp params from json to get A' mass eq_cfg = open(eq_cfg_file) eq_data = json.load(eq_cfg) for key, val in eq_data['SIMP Parameters'].items(): mass_ratio_Ap_to_Vd = eq_data['SIMP Parameters']['mass_ratio_Ap_to_Vd'] -print("LOOK HERE FUCKER", mass_ratio_Ap_to_Vd) m_Ap = options.mass*mass_ratio_Ap_to_Vd -print("here too idiot!!!!!!!! ", m_Ap) zbi.parameters['radFrac'] = radiativeFraction(m_Ap) zbi.parameters['radAcc'] = radiativeAcceptance(m_Ap) -zbi.parameters['bkgControlRegionFilename'] = '/sdf/group/hps/users/alspellm/projects/THESIS/data/2016/BLPass4_rereco_20230823/ana_20230823/output/final_hadd_data_BLPass4_ReRecon_20230828.root' -zbi.parameters['bkgControlRegionTreename'] = 'vtxana_kf_Tight_2016_simp_reach_CR' -zbi.parameters['dNdm_sf'] = 1.0 +zbi.parameters['bkgControlRegionFilename'] = '/sdf/group/hps/users/alspellm/projects/THESIS/data/2016/BLPass4c_20231006/ana_20231019/full_hadd_blpass4c_ana.root' +zbi.parameters['bkgControlRegionTreename'] = 'vtxana_Tight_2016_simp_reach_CR' +zbi.parameters['dNdm_sf'] = 10.0 # Sequence which the processors will run. p.sequence = [zbi] diff --git a/processors/include/FinalStateParticleProcessor.h b/processors/include/FinalStateParticleProcessor.h index 3e2c9a57d..c94082b8b 100644 --- a/processors/include/FinalStateParticleProcessor.h +++ b/processors/include/FinalStateParticleProcessor.h @@ -100,6 +100,8 @@ class FinalStateParticleProcessor : public Processor { std::string kinkRelCollLcio_{"GBLKinkDataRelations"}; //!< description std::string trkRelCollLcio_{"TrackDataRelations"}; //!< description std::string hitFitsCollLcio_{"SVTFittedRawTrackerHits"}; + + double bfield_{-1.}; //!< magnetic field int debug_{0}; //!< Debug Level diff --git a/processors/include/NewVertexAnaProcessor.h b/processors/include/NewVertexAnaProcessor.h new file mode 100644 index 000000000..111a0d1ba --- /dev/null +++ b/processors/include/NewVertexAnaProcessor.h @@ -0,0 +1,151 @@ +#ifndef __NEWVERTEX_ANAPROCESSOR_H__ +#define __NEWVERTEX_ANAPROCESSOR_H__ + +// HPSTR +#include "HpsEvent.h" +#include "Collections.h" +#include "EventHeader.h" +#include "TSData.h" +#include "Vertex.h" +#include "Track.h" +#include "TrackerHit.h" +#include "MCParticle.h" +#include "Particle.h" +#include "Processor.h" +#include "BaseSelector.h" +#include "TrackHistos.h" +#include "MCAnaHistos.h" +#include "utilities.h" + +#include "FlatTupleMaker.h" +#include "AnaHelpers.h" + +// ROOT +#include "TFile.h" +#include "TTree.h" +#include "TRefArray.h" +#include "TBranch.h" +#include "TVector3.h" +#include "TLorentzVector.h" + +// C++ +#include + +struct char_cmp { + bool operator () (const char *a,const char *b) const + { + return strcmp(a,b)<0; + } +}; + +/** + * @brief Insert description here. + * more details + */ +class NewVertexAnaProcessor : public Processor { + + public: + /** + * @brief Constructor + * + * @param name + * @param process + */ + NewVertexAnaProcessor(const std::string& name, Process& process); + + ~NewVertexAnaProcessor(); + + /** + * @brief description + * + * @param ievent + * @return true + * @return false + */ + virtual bool process(IEvent* ievent); + + /** + * @brief description + * + * @param tree + */ + virtual void initialize(TTree* tree); + + /** + * @brief description + * + */ + virtual void finalize(); + + /** + * @brief description + * + * @param parameters + */ + virtual void configure(const ParameterSet& parameters); + + private: + std::shared_ptr vtxSelector; //!< description + std::vector regionSelections_; //!< description + + std::string selectionCfg_; //!< description + std::map brMap_; //!< description + TBranch* bts_{nullptr}; //!< description + TBranch* bvtxs_{nullptr}; //!< description + TBranch* bhits_{nullptr}; //!< description + TBranch* bmcParts_{nullptr}; //!< description + TBranch* bevth_{nullptr}; //!< description + TBranch* becal_{nullptr}; //!< description + + EventHeader* evth_{nullptr}; //!< description + TSData* ts_{nullptr}; //!< description + std::vector* ecal_{}; //!< description + std::vector* vtxs_{}; //!< description + std::vector* hits_{}; //!< description + std::vector* mcParts_{}; //!< description + + std::string anaName_{"vtxAna"}; //!< description + std::string tsColl_{"TSBank"}; //!< description + std::string vtxColl_{"Vertices"}; //!< description + std::string hitColl_{"RotatedHelicalTrackHits"}; //!< description + std::string ecalColl_{"RecoEcalClusters"}; //!< description + std::string mcColl_{"MCParticle"}; //!< description + int isRadPDG_{622}; //!< description + int makeFlatTuple_{0}; //!< make true in config to save flat tuple + TTree* tree_{nullptr}; //!< description + + std::shared_ptr _vtx_histos; //!< description + std::shared_ptr _mc_vtx_histos; //!< description + + /** \todo Duplicate.. We can make a single class.. ? */ + std::map> _reg_vtx_selectors; //!< description + std::map> _reg_vtx_histos; //!< description + std::map> _reg_mc_vtx_histos; //!< description + std::map> _reg_tuples; //!< description + + std::vector _regions; //!< description + + typedef std::map>::iterator reg_it; //!< description + typedef std::map>::iterator reg_mc_it; //!< description + + std::string histoCfg_{""}; //!< description + std::string mcHistoCfg_{""}; //!< description + double timeOffset_{-999}; //!< description + double beamE_{2.3}; //!< In GeV. Default is 2016 value; + int isData_{0}; //!< description + std::string analysis_{"vertex"}; //!< description + + std::shared_ptr _ah; //!< description + + int debug_{0}; //!< Debug level + std::string beamPosCfg_{""}; //!< json containing run dep beamspot positions + json bpc_configs_; //!< json object + std::vector beamPosCorrections_ = {0.0,0.0,0.0}; //!< holds beam position corrections + std::string v0ProjectionFitsCfg_{""};//!< json file w run dependent v0 projection fits + json v0proj_fits_;//!< json object v0proj + double eleTrackTimeBias_ = 0.0; + double posTrackTimeBias_ = 0.0; + int current_run_number_{-999}; //!< track current run number +}; + +#endif diff --git a/processors/include/SimpZBiOptimizationProcessor.h b/processors/include/SimpZBiOptimizationProcessor.h index b0c5e9e0c..ffab91e5d 100644 --- a/processors/include/SimpZBiOptimizationProcessor.h +++ b/processors/include/SimpZBiOptimizationProcessor.h @@ -192,8 +192,6 @@ class SimpZBiOptimizationProcessor : public Processor { double highMass_; //* tracks_{}; TBranch* btracks_{nullptr}; //!< description - + + /** Event header branch. */ + TBranch* bevth_{nullptr}; //! + + /** Clusters */ + TBranch* becal_{nullptr}; //! + + // Event Header + EventHeader* evth_{nullptr}; //! + std::vector* ecal_{}; //!< + std::string trkCollName_; //!< Track Collection name + std::string ecalCollName_{"RecoEcalClusters"}; //!< Cluster Collection name // Track Selector configuration std::string selectionCfg_; std::shared_ptr trkSelector_; //!< description - + std::vector regionSelections_; //!< track selections + std::map> reg_selectors_; //!< description + std::map> reg_histos_; //!< description + typedef std::map>::iterator reg_it; //!< description + // Containers to hold histogrammer info std::string histCfgFilename_; //!< description std::string truthHistCfgFilename_; //!< description TrackHistos* trkHistos_{nullptr}; //!< description TrackHistos* truthHistos_{nullptr}; //!< description + + std::vector regions_; //! + bool doTruth_{false}; //!< description - + int isData_{1}; //! is data int debug_{0}; //!< debug level + float time_offset_{0}; //! time offset + + std::shared_ptr smearingTool_; + }; // TrackingAnaProcessor diff --git a/processors/include/TrackingProcessor.h b/processors/include/TrackingProcessor.h index 3eacbf574..664f6eecc 100644 --- a/processors/include/TrackingProcessor.h +++ b/processors/include/TrackingProcessor.h @@ -15,6 +15,7 @@ #include #include #include +#include #include #include #include @@ -99,7 +100,7 @@ class TrackingProcessor : public Processor { std::string kinkRelCollLcio_{"GBLKinkDataRelations"}; //!< collection name std::string trkRelCollLcio_{"TrackDataRelations"}; //!< collection name std::string trkCollRoot_{"GBLTracks"}; //!< collection name - + /** Container to hold all raw hits objecs. */ std::vector rawhits_{}; std::string hitFitsCollLcio_{"SVTFittedRawTrackerHits"}; //!< collection name diff --git a/processors/include/VertexAnaProcessor.h b/processors/include/VertexAnaProcessor.h index 33ea9320f..97ea2c7b4 100644 --- a/processors/include/VertexAnaProcessor.h +++ b/processors/include/VertexAnaProcessor.h @@ -144,9 +144,13 @@ class VertexAnaProcessor : public Processor { std::string beamPosCfg_{""}; //!< json containing run dep beamspot positions json bpc_configs_; //!< json object std::vector beamPosCorrections_ = {0.0,0.0,0.0}; //!< holds beam position corrections + std::string v0ProjectionFitsCfg_{""};//!< json file w run dependent v0 projection fits + json v0proj_fits_;//!< json object v0proj double eleTrackTimeBias_ = 0.0; double posTrackTimeBias_ = 0.0; int current_run_number_{-999}; //!< track current run number + + bool mc_reg_on_ = false; }; #endif diff --git a/processors/include/VertexProcessor.h b/processors/include/VertexProcessor.h index 537bf21a2..e1cc8e320 100644 --- a/processors/include/VertexProcessor.h +++ b/processors/include/VertexProcessor.h @@ -12,9 +12,13 @@ // LCIO // //----------// #include +#include +#include #include #include #include +#include +#include //-----------// // hpstr // @@ -23,6 +27,8 @@ #include "Vertex.h" #include "Particle.h" #include "Event.h" +#include "TrackerHit.h" +#include "RawSvtHit.h" // Forward declarations class TTree; @@ -77,6 +83,10 @@ class VertexProcessor : public Processor { private: /** Containers to hold all TrackerHit objects. */ + std::vector hits_{}; + std::string trkhitCollRoot_{"fspOnTrackHits"}; + std::vector rawhits_{}; + std::string rawhitCollRoot_{"fspOnTrackRawHits"}; std::vector vtxs_{}; //!< description std::vector parts_{}; //!< description std::string vtxCollLcio_{"UnconstrainedV0Vertices"}; //!< description @@ -84,8 +94,12 @@ class VertexProcessor : public Processor { std::string partCollRoot_{"ParticlesOnVertices"}; //!< description std::string kinkRelCollLcio_{"GBLKinkDataRelations"}; //!< description std::string trkRelCollLcio_{"TrackDataRelations"}; //!< description + std::string hitFitsCollLcio_{"SVTFittedRawTrackerHits"}; + std::string mcPartRelLcio_{"RotatedHelicalTrackMCRelations"}; std::string trackStateLocation_{""}; //!< select track state for tracks DEFAULT AtIP + double bfield_{-1.}; //!< magnetic field + int debug_{0}; //!< Debug Level }; // VertexProcessor diff --git a/processors/include/utilities.h b/processors/include/utilities.h index 7803c8464..6da65c5a6 100644 --- a/processors/include/utilities.h +++ b/processors/include/utilities.h @@ -18,6 +18,10 @@ #include #include +#include +#include +#include +using json = nlohmann::json; //-----------// // hpstr // @@ -184,7 +188,16 @@ namespace utils { * * \todo extern? */ - void get2016KFMCTruthHitCodes(Track* ele_trk, Track* pos_trk, std::vector* hits, int& L1L2hitCode, int& L1hitCode, int& L2hitCode); + void get2016KFMCTruthHitCodes(Track* ele_trk, Track* pos_trk, int& L1L2hitCode, int& L1hitCode, int& L2hitCode); + + /** + * @brief description + * + * \todo extern? + */ + double v0_projection_to_target_significance(json v0proj_fits, int run, double &vtx_proj_x, double &vtx_proj_y, + double &vtx_proj_x_signif, double &vtx_proj_y_signif, double vtx_x, double vtx_y, double vtx_z, + double vtx_px, double vtx_py, double vtx_pz); } #endif //UTILITIES diff --git a/processors/src/FinalStateParticleProcessor.cxx b/processors/src/FinalStateParticleProcessor.cxx index 1c73b0f99..f3920cbe9 100644 --- a/processors/src/FinalStateParticleProcessor.cxx +++ b/processors/src/FinalStateParticleProcessor.cxx @@ -26,6 +26,7 @@ void FinalStateParticleProcessor::configure(const ParameterSet& parameters) { hitFitsCollLcio_ = parameters.getString("hitFitsCollLcio", hitFitsCollLcio_); trkhitCollRoot_ = parameters.getString("trkhitCollRoot",trkhitCollRoot_); rawhitCollRoot_ = parameters.getString("rawhitCollRoot",rawhitCollRoot_); + bfield_ = parameters.getDouble("bfield",bfield_); } catch (std::runtime_error& error) { @@ -129,6 +130,7 @@ bool FinalStateParticleProcessor::process(IEvent* ievent) { if (lc_fsp->getTracks().size()>0){ EVENT::Track* lc_track = static_cast(lc_fsp->getTracks()[0]); Track* track = utils::buildTrack(lc_track,"",gbl_kink_data,track_data); + if (bfield_ > 0.0) track->setMomentum(bfield_); if (track->isKalmanTrack()) hitType = 1; //SiClusters EVENT::TrackerHitVec lc_tracker_hits = lc_track->getTrackerHits(); for (auto lc_tracker_hit : lc_tracker_hits) { diff --git a/processors/src/NewVertexAnaProcessor.cxx b/processors/src/NewVertexAnaProcessor.cxx new file mode 100644 index 000000000..8b379beee --- /dev/null +++ b/processors/src/NewVertexAnaProcessor.cxx @@ -0,0 +1,1405 @@ +/** + *@file NewVertexAnaProcessor.cxx + *@brief Main vertex analysis processor + *@author PF, SLAC + */ + +#include "NewVertexAnaProcessor.h" +#include +#include +#include + +NewVertexAnaProcessor::NewVertexAnaProcessor(const std::string& name, Process& process) : Processor(name,process) { + +} + +//TODO Check this destructor + +NewVertexAnaProcessor::~NewVertexAnaProcessor(){} + +void NewVertexAnaProcessor::configure(const ParameterSet& parameters) { + std::cout << "Configuring NewVertexAnaProcessor" <(); + + vtxSelector = std::make_shared(anaName_+"_"+"vtxSelection",selectionCfg_); + vtxSelector->setDebug(debug_); + vtxSelector->LoadSelection(); + + _vtx_histos = std::make_shared(anaName_+"_"+"vtxSelection"); + _vtx_histos->loadHistoConfig(histoCfg_); + _vtx_histos->DefineHistos(); + + if(!isData_){ + _mc_vtx_histos = std::make_shared(anaName_+"_mc_"+"vtxSelection"); + _mc_vtx_histos->loadHistoConfig(mcHistoCfg_); + _mc_vtx_histos->DefineHistos(); + _mc_vtx_histos->Define2DHistos(); + } + + //Load Run Dependent V0 target projection fits from json + if(!v0ProjectionFitsCfg_.empty()){ + std::ifstream v0proj_file(v0ProjectionFitsCfg_); + v0proj_file >> v0proj_fits_; + v0proj_file.close(); + } + + //Run Dependent Corrections + //Beam Position + if(!beamPosCfg_.empty()){ + std::ifstream bpc_file(beamPosCfg_); + bpc_file >> bpc_configs_; + bpc_file.close(); + } + // histos = new MCAnaHistos(anaName_); + //histos->loadHistoConfig(histCfgFilename_) + //histos->DefineHistos(); + //histos->Define2DHistos(); + + + //For each region initialize plots + + for (unsigned int i_reg = 0; i_reg < regionSelections_.size(); i_reg++) { + std::string regname = AnaHelpers::getFileName(regionSelections_[i_reg],false); + std::cout<<"Setting up region:: " << regname <(anaName_+"_"+regname, regionSelections_[i_reg]); + _reg_vtx_selectors[regname]->setDebug(debug_); + _reg_vtx_selectors[regname]->LoadSelection(); + + _reg_vtx_histos[regname] = std::make_shared(anaName_+"_"+regname); + _reg_vtx_histos[regname]->loadHistoConfig(histoCfg_); + _reg_vtx_histos[regname]->DefineHistos(); + + + if(!isData_){ + _reg_mc_vtx_histos[regname] = std::make_shared(anaName_+"_mc_"+regname); + _reg_mc_vtx_histos[regname]->loadHistoConfig(mcHistoCfg_); + _reg_mc_vtx_histos[regname]->DefineHistos(); + } + + //Build a flat tuple for vertex and track params + if (makeFlatTuple_){ + _reg_tuples[regname] = std::make_shared(anaName_+"_"+regname+"_tree"); + + //vtx vars + _reg_tuples[regname]->addVariable("unc_vtx_mass"); + _reg_tuples[regname]->addVariable("unc_vtx_z"); + _reg_tuples[regname]->addVariable("unc_vtx_chi2"); + _reg_tuples[regname]->addVariable("unc_vtx_psum"); + _reg_tuples[regname]->addVariable("unc_vtx_px"); + _reg_tuples[regname]->addVariable("unc_vtx_py"); + _reg_tuples[regname]->addVariable("unc_vtx_pz"); + _reg_tuples[regname]->addVariable("unc_vtx_x"); + _reg_tuples[regname]->addVariable("unc_vtx_y"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_pos_clus_dt"); + _reg_tuples[regname]->addVariable("run_number"); + _reg_tuples[regname]->addVariable("unc_vtx_cxx"); + _reg_tuples[regname]->addVariable("unc_vtx_cyy"); + _reg_tuples[regname]->addVariable("unc_vtx_czz"); + _reg_tuples[regname]->addVariable("unc_vtx_cyx"); + _reg_tuples[regname]->addVariable("unc_vtx_czy"); + _reg_tuples[regname]->addVariable("unc_vtx_czx"); + _reg_tuples[regname]->addVariable("unc_vtx_proj_x"); + _reg_tuples[regname]->addVariable("unc_vtx_proj_y"); + _reg_tuples[regname]->addVariable("unc_vtx_proj_x_sig"); + _reg_tuples[regname]->addVariable("unc_vtx_proj_y_sig"); + _reg_tuples[regname]->addVariable("unc_vtx_proj_sig"); + _reg_tuples[regname]->addVariable("unc_vtx_deltaZ"); + + + //track vars + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_p"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_t"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_d0"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_phi0"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_omega"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_tanLambda"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_z0"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_chi2ndf"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_clust_dt"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_z0Err"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_d0Err"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_tanLambdaErr"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_PhiErr"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_OmegaErr"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_L1_isolation"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_nhits"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_lastlayer"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_si0"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_si1"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_ecal_x"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_ecal_y"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_z"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_px"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_py"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_track_pz"); + + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_clust_dt"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_p"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_t"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_d0"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_phi0"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_omega"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_tanLambda"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_z0"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_chi2ndf"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_z0Err"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_d0Err"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_tanLambdaErr"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_PhiErr"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_OmegaErr"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_L1_isolation"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_nhits"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_lastlayer"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_si0"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_si1"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_ecal_x"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_ecal_y"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_z"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_px"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_py"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_track_pz"); + + //clust vars + _reg_tuples[regname]->addVariable("unc_vtx_ele_clust_E"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_clust_x"); + _reg_tuples[regname]->addVariable("unc_vtx_ele_clust_corr_t"); + + _reg_tuples[regname]->addVariable("unc_vtx_pos_clust_E"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_clust_x"); + _reg_tuples[regname]->addVariable("unc_vtx_pos_clust_corr_t"); + + if(!isData_) + { + _reg_tuples[regname]->addVariable("true_vtx_z"); + _reg_tuples[regname]->addVariable("true_vtx_mass"); + _reg_tuples[regname]->addVariable("ap_true_vtx_z"); + _reg_tuples[regname]->addVariable("ap_true_vtx_mass"); + _reg_tuples[regname]->addVariable("ap_true_vtx_energy"); + _reg_tuples[regname]->addVariable("vd_true_vtx_z"); + _reg_tuples[regname]->addVariable("vd_true_vtx_mass"); + _reg_tuples[regname]->addVariable("vd_true_vtx_energy"); + _reg_tuples[regname]->addVariable("hitCode"); + _reg_tuples[regname]->addVariable("L1hitCode"); + _reg_tuples[regname]->addVariable("L2hitCode"); + } + } + + _regions.push_back(regname); + } + + // Get list of branches in tree to help protect accessing them + int nBr = tree_->GetListOfBranches()->GetEntries(); + if (debug_) std::cout << "Tree has " << nBr << " branches" << std::endl; + for(int iBr = 0; iBr < nBr; iBr++) + { + TBranch *br = dynamic_cast(tree_->GetListOfBranches()->At(iBr)); + brMap_.insert(std::map::value_type(br->GetName(), 1)); + if (debug_) std::cout << br->GetName() << ": " << brMap_[br->GetName()] << std::endl; + } + + //init Reading Tree + tree_->SetBranchAddress("EventHeader", &evth_ , &bevth_); + if (brMap_.find(tsColl_.c_str()) != brMap_.end()) tree_->SetBranchAddress(tsColl_.c_str(), &ts_ , &bts_); + tree_->SetBranchAddress(vtxColl_.c_str(), &vtxs_ , &bvtxs_); + //tree_->SetBranchAddress(hitColl_.c_str(), &hits_ , &bhits_); + if (brMap_.find(hitColl_.c_str()) != brMap_.end()) tree_->SetBranchAddress(hitColl_.c_str(), &hits_ , &bhits_); + if(!isData_ && !mcColl_.empty()) tree_->SetBranchAddress(mcColl_.c_str() , &mcParts_, &bmcParts_); +} + +bool NewVertexAnaProcessor::process(IEvent* ievent) { + if(debug_) { + std:: cout << "----------------- Event " << evth_->getEventNumber() << " -----------------" << std::endl; + } + HpsEvent* hps_evt = (HpsEvent*) ievent; + double weight = 1.; + int run_number = evth_->getRunNumber(); + int closest_run; + if (debug_) std::cout << "Check pbc_configs" << std::endl; + if(!bpc_configs_.empty()){ + for(auto run : bpc_configs_.items()){ + int check_run = std::stoi(run.key()); + if(check_run > run_number) + break; + else{ + closest_run = check_run; + } + } + beamPosCorrections_ = {bpc_configs_[std::to_string(closest_run)]["beamspot_x"], + bpc_configs_[std::to_string(closest_run)]["beamspot_y"], + bpc_configs_[std::to_string(closest_run)]["beamspot_z"]}; + } + + + //Get "true" values + //AP + double apMass = -0.9; + double apZ = -0.9; + double apEnergy = -0.9; + //Simp + double vdMass = -0.9; + double vdZ = -0.9; + double vdEnergy = -0.9; + + if (debug_) std::cout << "plot trigger info" << std::endl; + //Plot info about which trigger bits are present in the event + if (ts_ != nullptr) + { + _vtx_histos->Fill2DHisto("trig_count_hh", + ((int)ts_->prescaled.Single_3_Top)+((int)ts_->prescaled.Single_3_Bot), + ((int)ts_->prescaled.Single_2_Top)+((int)ts_->prescaled.Single_2_Bot)); + } + if (debug_) std::cout << "plot vtx N" << std::endl; + _vtx_histos->Fill1DHisto("n_vtx_h", vtxs_->size()); + + if (mcParts_) { + for(int i = 0; i < mcParts_->size(); i++) + { + if(mcParts_->at(i)->getPDG() == 622) + { + apMass = mcParts_->at(i)->getMass(); + apZ = mcParts_->at(i)->getVertexPosition().at(2); + apEnergy = mcParts_->at(i)->getEnergy(); + } + if(mcParts_->at(i)->getPDG() == 625) + { + vdMass = mcParts_->at(i)->getMass(); + vdZ = mcParts_->at(i)->getVertexPosition().at(2); + vdEnergy = mcParts_->at(i)->getEnergy(); + } + } + + if (!isData_) _mc_vtx_histos->FillMCParticles(mcParts_, analysis_); + } + //Store processed number of events + std::vector selected_vtxs; + bool passVtxPresel = false; + + if(debug_){ + std::cout<<"Number of vertices found in event: "<< vtxs_->size()<size(); i_vtx++ ) { + vtxSelector->getCutFlowHisto()->Fill(0.,weight); + + Vertex* vtx = vtxs_->at(i_vtx); + Particle* ele = nullptr; + Track* ele_trk = nullptr; + Particle* pos = nullptr; + Track* pos_trk = nullptr; + + //Trigger requirement - *really hate* having to do it here for each vertex. + + if (isData_) { + if (!vtxSelector->passCutEq("Pair1_eq",(int)evth_->isPair1Trigger(),weight)) + break; + } + + bool foundParts = _ah->GetParticlesFromVtx(vtx,ele,pos); + if (!foundParts) { + if(debug_) std::cout<<"NewVertexAnaProcessor::WARNING::Found vtx without ele/pos. Skip."<getTrack().Clone(); + pos_trk = (Track*) pos->getTrack().Clone(); + + //Beam Position Corrections + ele_trk->applyCorrection("z0", beamPosCorrections_.at(1)); + pos_trk->applyCorrection("z0", beamPosCorrections_.at(1)); + //Track Time Corrections + ele_trk->applyCorrection("track_time",eleTrackTimeBias_); + pos_trk->applyCorrection("track_time", posTrackTimeBias_); + + //Add the momenta to the tracks - do not do that + //ele_trk->setMomentum(ele->getMomentum()[0],ele->getMomentum()[1],ele->getMomentum()[2]); + //pos_trk->setMomentum(pos->getMomentum()[0],pos->getMomentum()[1],pos->getMomentum()[2]); + + double ele_E = ele->getEnergy(); + double pos_E = pos->getEnergy(); + if (debug_) std::cout << "got tracks" << std::endl; + + CalCluster eleClus = ele->getCluster(); + CalCluster posClus = pos->getCluster(); + + + //Compute analysis variables here. + TLorentzVector p_ele; + p_ele.SetPxPyPzE(ele_trk->getMomentum()[0],ele_trk->getMomentum()[1],ele_trk->getMomentum()[2],ele->getEnergy()); + TLorentzVector p_pos; + p_pos.SetPxPyPzE(pos_trk->getMomentum()[0],pos_trk->getMomentum()[1],pos_trk->getMomentum()[2],ele->getEnergy()); + + //Tracks in opposite volumes - useless + //if (!vtxSelector->passCutLt("eleposTanLambaProd_lt",ele_trk->getTanLambda() * pos_trk->getTanLambda(),weight)) + // continue; + + if (debug_) std::cout << "start selection" << std::endl; + //Ele Track Time + if (!vtxSelector->passCutLt("eleTrkTime_lt",fabs(ele_trk->getTrackTime()),weight)) + continue; + + //Pos Track Time + if (!vtxSelector->passCutLt("posTrkTime_lt",fabs(pos_trk->getTrackTime()),weight)) + continue; + + //Ele Track-cluster match + if (!vtxSelector->passCutLt("eleTrkCluMatch_lt",ele->getGoodnessOfPID(),weight)) + continue; + + //Pos Track-cluster match + if (!vtxSelector->passCutLt("posTrkCluMatch_lt",pos->getGoodnessOfPID(),weight)) + continue; + + //Require Positron Cluster exists + if (!vtxSelector->passCutGt("posClusE_gt",posClus.getEnergy(),weight)) + continue; + + //Require Positron Cluster does NOT exists + if (!vtxSelector->passCutLt("posClusE_lt",posClus.getEnergy(),weight)) + continue; + + + double corr_eleClusterTime = ele->getCluster().getTime() - timeOffset_; + double corr_posClusterTime = pos->getCluster().getTime() - timeOffset_; + + double botClusTime = 0.0; + if(ele->getCluster().getPosition().at(1) < 0.0) botClusTime = ele->getCluster().getTime(); + else botClusTime = pos->getCluster().getTime(); + + //Bottom Cluster Time + if (!vtxSelector->passCutLt("botCluTime_lt", botClusTime, weight)) + continue; + + if (!vtxSelector->passCutGt("botCluTime_gt", botClusTime, weight)) + continue; + + //Ele Pos Cluster Time Difference + if (!vtxSelector->passCutLt("eleposCluTimeDiff_lt",fabs(corr_eleClusterTime - corr_posClusterTime),weight)) + continue; + + //Ele Track-Cluster Time Difference + if (!vtxSelector->passCutLt("eleTrkCluTimeDiff_lt",fabs(ele_trk->getTrackTime() - corr_eleClusterTime),weight)) + continue; + + //Pos Track-Cluster Time Difference + if (!vtxSelector->passCutLt("posTrkCluTimeDiff_lt",fabs(pos_trk->getTrackTime() - corr_posClusterTime),weight)) + continue; + + TVector3 ele_mom; + //ele_mom.SetX(ele->getMomentum()[0]); + //ele_mom.SetY(ele->getMomentum()[1]); + //ele_mom.SetZ(ele->getMomentum()[2]); + ele_mom.SetX(ele_trk->getMomentum()[0]); + ele_mom.SetY(ele_trk->getMomentum()[1]); + ele_mom.SetZ(ele_trk->getMomentum()[2]); + + + TVector3 pos_mom; + //pos_mom.SetX(pos->getMomentum()[0]); + //pos_mom.SetY(pos->getMomentum()[1]); + //pos_mom.SetZ(pos->getMomentum()[2]); + pos_mom.SetX(pos_trk->getMomentum()[0]); + pos_mom.SetY(pos_trk->getMomentum()[1]); + pos_mom.SetZ(pos_trk->getMomentum()[2]); + + + + //Ele Track Quality - Chi2 + if (!vtxSelector->passCutLt("eleTrkChi2_lt",ele_trk->getChi2(),weight)) + continue; + + //Pos Track Quality - Chi2 + if (!vtxSelector->passCutLt("posTrkChi2_lt",pos_trk->getChi2(),weight)) + continue; + + //Ele Track Quality - Chi2Ndf + if (!vtxSelector->passCutLt("eleTrkChi2Ndf_lt",ele_trk->getChi2Ndf(),weight)) + continue; + + //Pos Track Quality - Chi2Ndf + if (!vtxSelector->passCutLt("posTrkChi2Ndf_lt",pos_trk->getChi2Ndf(),weight)) + continue; + + //Beam Electron cut + if (!vtxSelector->passCutLt("eleMom_lt",ele_mom.Mag(),weight)) + continue; + + //Ele min momentum cut + if (!vtxSelector->passCutGt("eleMom_gt",ele_mom.Mag(),weight)) + continue; + + //Pos min momentum cut + if (!vtxSelector->passCutGt("posMom_gt",pos_mom.Mag(),weight)) + continue; + + //Ele nHits + int ele2dHits = ele_trk->getTrackerHitCount(); + if (!ele_trk->isKalmanTrack()) + ele2dHits*=2; + + if (!vtxSelector->passCutGt("eleN2Dhits_gt",ele2dHits,weight)) { + continue; + } + + //Pos nHits + int pos2dHits = pos_trk->getTrackerHitCount(); + if (!pos_trk->isKalmanTrack()) + pos2dHits*=2; + + if (!vtxSelector->passCutGt("posN2Dhits_gt",pos2dHits,weight)) { + continue; + } + + //Less than 4 shared hits for ele/pos track + if (!vtxSelector->passCutLt("eleNshared_lt",ele_trk->getNShared(),weight)) { + continue; + } + + if (!vtxSelector->passCutLt("posNshared_lt",pos_trk->getNShared(),weight)) { + continue; + } + + + //Vertex Quality + if (!vtxSelector->passCutLt("chi2unc_lt",vtx->getChi2(),weight)) + continue; + + //Max vtx momentum + if (!vtxSelector->passCutLt("maxVtxMom_lt",(ele_mom+pos_mom).Mag(),weight)) + continue; + + //Min vtx momentum + + if (!vtxSelector->passCutGt("minVtxMom_gt",(ele_mom+pos_mom).Mag(),weight)) + continue; + + if (debug_) std::cout << "fill 1D Vertex" << std::endl; + _vtx_histos->Fill1DVertex(vtx, + ele, + pos, + ele_trk, + pos_trk, + weight); + + if (debug_) std::cout << "fill track histos" << std::endl; + double ele_pos_dt = corr_eleClusterTime - corr_posClusterTime; + double psum = ele_mom.Mag()+pos_mom.Mag(); + + _vtx_histos->Fill1DTrack(ele_trk,weight, "ele_"); + _vtx_histos->Fill1DTrack(pos_trk,weight, "pos_"); + _vtx_histos->Fill1DHisto("ele_track_n2dhits_h", ele2dHits, weight); + _vtx_histos->Fill1DHisto("pos_track_n2dhits_h", pos2dHits, weight); + _vtx_histos->Fill1DHisto("vtx_Psum_h", p_ele.P()+p_pos.P(), weight); + _vtx_histos->Fill1DHisto("vtx_Esum_h", ele_E + pos_E, weight); + _vtx_histos->Fill1DHisto("ele_pos_clusTimeDiff_h", (corr_eleClusterTime - corr_posClusterTime), weight); + _vtx_histos->Fill2DHisto("ele_vtxZ_iso_hh", TMath::Min(ele_trk->getIsolation(0), ele_trk->getIsolation(1)), vtx->getZ(), weight); + _vtx_histos->Fill2DHisto("pos_vtxZ_iso_hh", TMath::Min(pos_trk->getIsolation(0), pos_trk->getIsolation(1)), vtx->getZ(), weight); + _vtx_histos->Fill2DHistograms(vtx,weight); + _vtx_histos->Fill2DTrack(ele_trk,weight,"ele_"); + _vtx_histos->Fill2DTrack(pos_trk,weight,"pos_"); + _vtx_histos->Fill1DHisto("mcMass622_h",apMass); + _vtx_histos->Fill1DHisto("mcZ622_h",apZ); + + //New SIMP histos for developing loose preselection cuts + //2d histos + _vtx_histos->Fill2DHisto("ele_clusT_v_ele_trackT_hh", ele_trk->getTrackTime(), corr_eleClusterTime, weight); + _vtx_histos->Fill2DHisto("pos_clusT_v_pos_trackT_hh", pos_trk->getTrackTime(), corr_posClusterTime, weight); + _vtx_histos->Fill2DHisto("ele_track_time_v_P_hh", ele_trk->getP(), ele_trk->getTrackTime(), weight); + _vtx_histos->Fill2DHisto("pos_track_time_v_P_hh", pos_trk->getP(), pos_trk->getTrackTime(), weight); + _vtx_histos->Fill2DHisto("ele_pos_clusTimeDiff_v_pSum_hh",ele_mom.Mag()+pos_mom.Mag(), ele_pos_dt, weight); + _vtx_histos->Fill2DHisto("ele_cluster_energy_v_track_p_hh",ele_trk->getP(), eleClus.getEnergy(), weight); + _vtx_histos->Fill2DHisto("pos_cluster_energy_v_track_p_hh",pos_trk->getP(), posClus.getEnergy(), weight); + _vtx_histos->Fill2DHisto("ele_track_cluster_dt_v_EoverP_hh",eleClus.getEnergy()/ele_trk->getP(), ele_trk->getTrackTime() - corr_eleClusterTime, weight); + _vtx_histos->Fill2DHisto("pos_track_cluster_dt_v_EoverP_hh",posClus.getEnergy()/pos_trk->getP(), pos_trk->getTrackTime() - corr_posClusterTime, weight); + _vtx_histos->Fill2DHisto("ele_track_clus_dt_v_p_hh",ele_trk->getP(), ele_trk->getTrackTime() - corr_eleClusterTime, weight); + _vtx_histos->Fill2DHisto("pos_track_clus_dt_v_p_hh",pos_trk->getP(), pos_trk->getTrackTime() - corr_posClusterTime, weight); + //chi2 2d plots + _vtx_histos->Fill2DHisto("ele_track_chi2ndf_v_time_hh", ele_trk->getTrackTime(), ele_trk->getChi2Ndf(), weight); + _vtx_histos->Fill2DHisto("ele_track_chi2ndf_v_p_hh", ele_trk->getP(), ele_trk->getChi2Ndf(), weight); + _vtx_histos->Fill2DHisto("ele_track_chi2ndf_v_tanlambda_hh", ele_trk->getTanLambda(), ele_trk->getChi2Ndf(), weight); + _vtx_histos->Fill2DHisto("ele_track_chi2ndf_v_n2dhits_hh", ele2dHits, ele_trk->getChi2Ndf(), weight); + + _vtx_histos->Fill2DHisto("pos_track_chi2ndf_v_time_hh", pos_trk->getTrackTime(), pos_trk->getChi2Ndf(), weight); + _vtx_histos->Fill2DHisto("pos_track_chi2ndf_v_p_hh", pos_trk->getP(), pos_trk->getChi2Ndf(), weight); + _vtx_histos->Fill2DHisto("pos_track_chi2ndf_v_tanlambda_hh", pos_trk->getTanLambda(), pos_trk->getChi2Ndf(), weight); + _vtx_histos->Fill2DHisto("pos_track_chi2ndf_v_n2dhits_hh", pos2dHits, pos_trk->getChi2Ndf(), weight); + + _vtx_histos->Fill2DHisto("ele_track_p_v_tanlambda_hh", ele_trk->getTanLambda(), ele_trk->getP(), weight); + _vtx_histos->Fill2DHisto("pos_track_p_v_tanlambda_hh", pos_trk->getTanLambda(), pos_trk->getP(), weight); + + + //1d histos + _vtx_histos->Fill1DHisto("ele_track_clus_dt_h", ele_trk->getTrackTime() - corr_eleClusterTime, weight); + _vtx_histos->Fill1DHisto("pos_track_clus_dt_h", pos_trk->getTrackTime() - corr_posClusterTime, weight); + + passVtxPresel = true; + + selected_vtxs.push_back(vtx); + vtxSelector->clearSelector(); + } + + // std::cout << "Number of selected vtxs: " << selected_vtxs.size() << std::endl; + + _vtx_histos->Fill1DHisto("n_vertices_h",selected_vtxs.size()); + + + //not working atm + //hps_evt->addVertexCollection("selected_vtxs", selected_vtxs); + + //Make Plots for each region: loop on each region. Check if the region has the cut and apply it + //TODO Clean this up => Cuts should be implemented in each region? + //TODO Bring the preselection out of this stupid loop + + + if (debug_) std::cout << "start regions" << std::endl; + //TODO add yields. => Quite terrible way to loop. + for (auto region : _regions ) { + + int nGoodVtx = 0; + Vertex* goodVtx = nullptr; + std::vector goodVtxs; + + float truePsum = -1; + float trueEsum = -1; + + for ( auto vtx : selected_vtxs) { + + //No cuts. + _reg_vtx_selectors[region]->getCutFlowHisto()->Fill(0.,weight); + + + Particle* ele = nullptr; + Particle* pos = nullptr; + + _ah->GetParticlesFromVtx(vtx,ele,pos); + + CalCluster eleClus = ele->getCluster(); + CalCluster posClus = pos->getCluster(); + //vtx X position + if (!_reg_vtx_selectors[region]->passCutLt("uncVtxX_lt",fabs(vtx->getX()),weight)) + continue; + + //vtx Y position + if (!_reg_vtx_selectors[region]->passCutLt("uncVtxY_lt",fabs(vtx->getY()),weight)) + continue; + + //vtx Z position + if (!_reg_vtx_selectors[region]->passCutGt("uncVtxZ_gt",vtx->getZ(),weight)) + continue; + + double ele_E = ele->getEnergy(); + double pos_E = pos->getEnergy(); + + //Compute analysis variables here. + + Track * ele_trk = nullptr; + ele_trk = (Track*) ele->getTrack().Clone(); + + Track * pos_trk = nullptr; + pos_trk = (Track*) pos->getTrack().Clone(); + + //Beam Position Corrections + ele_trk->applyCorrection("z0", beamPosCorrections_.at(1)); + pos_trk->applyCorrection("z0", beamPosCorrections_.at(1)); + //Track Time Corrections + ele_trk->applyCorrection("track_time",eleTrackTimeBias_); + pos_trk->applyCorrection("track_time", posTrackTimeBias_); + + //Add the momenta to the tracks + //ele_trk->setMomentum(ele->getMomentum()[0],ele->getMomentum()[1],ele->getMomentum()[2]); + //pos_trk->setMomentum(pos->getMomentum()[0],pos->getMomentum()[1],pos->getMomentum()[2]); + TVector3 recEleP(ele->getMomentum()[0],ele->getMomentum()[1],ele->getMomentum()[2]); + TLorentzVector p_ele; + p_ele.SetPxPyPzE(ele_trk->getMomentum()[0],ele_trk->getMomentum()[1],ele_trk->getMomentum()[2], ele_E); + TLorentzVector p_pos; + p_pos.SetPxPyPzE(pos_trk->getMomentum()[0],pos_trk->getMomentum()[1],pos_trk->getMomentum()[2], pos_E); + + //Get the layers hit on each track + std::vector ele_hit_layers = ele_trk->getHitLayers(); + int ele_Si0 = 0; + int ele_Si1 = 0; + int ele_lastlayer = 0; + for(int i=0; i pos_hit_layers = pos_trk->getHitLayers(); + int pos_Si0 = 0; + int pos_Si1 = 0; + int pos_lastlayer = 0; + for(int i=0; igetTrackerHitCount()<InnermostLayerCheck(ele_trk, foundL1ele, foundL2ele); + + + if (debug_) { + std::cout<<"Check on pos_Track"<getTrackerHitCount()<InnermostLayerCheck(pos_trk, foundL1pos, foundL2pos); + + if (debug_) { + std::cout<<"Check on pos_Track"<passCutEq("Pair1_eq",(int)evth_->isPair1Trigger(),weight)) + break; + } + //Ele Track Time + if (!_reg_vtx_selectors[region]->passCutLt("eleTrkTime_lt",fabs(ele_trk->getTrackTime()),weight)) + continue; + + //Pos Track Time + if (!_reg_vtx_selectors[region]->passCutLt("posTrkTime_lt",fabs(pos_trk->getTrackTime()),weight)) + continue; + + //Ele Track-cluster match + if (!_reg_vtx_selectors[region]->passCutLt("eleTrkCluMatch_lt",ele->getGoodnessOfPID(),weight)) + continue; + + //Pos Track-cluster match + if (!_reg_vtx_selectors[region]->passCutLt("posTrkCluMatch_lt",pos->getGoodnessOfPID(),weight)) + continue; + + //Require Positron Cluster exists + if (!_reg_vtx_selectors[region]->passCutGt("posClusE_gt",posClus.getEnergy(),weight)) + continue; + + //Require Positron Cluster does NOT exists + if (!_reg_vtx_selectors[region]->passCutLt("posClusE_lt",posClus.getEnergy(),weight)) + continue; + + double corr_eleClusterTime = ele->getCluster().getTime() - timeOffset_; + double corr_posClusterTime = pos->getCluster().getTime() - timeOffset_; + + double botClusTime = 0.0; + if(ele->getCluster().getPosition().at(1) < 0.0) botClusTime = ele->getCluster().getTime(); + else botClusTime = pos->getCluster().getTime(); + + //Bottom Cluster Time + if (!_reg_vtx_selectors[region]->passCutLt("botCluTime_lt", botClusTime, weight)) + continue; + + if (!_reg_vtx_selectors[region]->passCutGt("botCluTime_gt", botClusTime, weight)) + continue; + + //Ele Pos Cluster Time Difference + if (!_reg_vtx_selectors[region]->passCutLt("eleposCluTimeDiff_lt",fabs(corr_eleClusterTime - corr_posClusterTime),weight)) + continue; + + //Ele Track-Cluster Time Difference + if (!_reg_vtx_selectors[region]->passCutLt("eleTrkCluTimeDiff_lt",fabs(ele_trk->getTrackTime() - corr_eleClusterTime),weight)) + continue; + + //Pos Track-Cluster Time Difference + if (!_reg_vtx_selectors[region]->passCutLt("posTrkCluTimeDiff_lt",fabs(pos_trk->getTrackTime() - corr_posClusterTime),weight)) + continue; + + TVector3 ele_mom; + ele_mom.SetX(ele_trk->getMomentum()[0]); + ele_mom.SetY(ele_trk->getMomentum()[1]); + ele_mom.SetZ(ele_trk->getMomentum()[2]); + + + TVector3 pos_mom; + pos_mom.SetX(pos_trk->getMomentum()[0]); + pos_mom.SetY(pos_trk->getMomentum()[1]); + pos_mom.SetZ(pos_trk->getMomentum()[2]); + + //Ele Track Quality - Chi2 + if (!_reg_vtx_selectors[region]->passCutLt("eleTrkChi2_lt",ele_trk->getChi2(),weight)) + continue; + + //Pos Track Quality - Chi2 + if (!_reg_vtx_selectors[region]->passCutLt("posTrkChi2_lt",pos_trk->getChi2(),weight)) + continue; + + //Ele Track Quality - Chi2Ndf + if (!_reg_vtx_selectors[region]->passCutLt("eleTrkChi2Ndf_lt",ele_trk->getChi2Ndf(),weight)) + continue; + + //Pos Track Quality - Chi2Ndf + if (!_reg_vtx_selectors[region]->passCutLt("posTrkChi2Ndf_lt",pos_trk->getChi2Ndf(),weight)) + continue; + + //Beam Electron cut + if (!_reg_vtx_selectors[region]->passCutLt("eleMom_lt",ele_mom.Mag(),weight)) + continue; + + //Ele min momentum cut + if (!_reg_vtx_selectors[region]->passCutGt("eleMom_gt",ele_mom.Mag(),weight)) + continue; + + //Pos min momentum cut + if (!_reg_vtx_selectors[region]->passCutGt("posMom_gt",pos_mom.Mag(),weight)) + continue; + + //Ele nHits + int ele2dHits = ele_trk->getTrackerHitCount(); + if (!ele_trk->isKalmanTrack()) + ele2dHits*=2; + + if (!_reg_vtx_selectors[region]->passCutGt("eleN2Dhits_gt",ele2dHits,weight)) { + continue; + } + + //Pos nHits + int pos2dHits = pos_trk->getTrackerHitCount(); + if (!pos_trk->isKalmanTrack()) + pos2dHits*=2; + + if (!_reg_vtx_selectors[region]->passCutGt("posN2Dhits_gt",pos2dHits,weight)) { + continue; + } + + //Less than 4 shared hits for ele/pos track + if (!_reg_vtx_selectors[region]->passCutLt("eleNshared_lt",ele_trk->getNShared(),weight)) { + continue; + } + + if (!_reg_vtx_selectors[region]->passCutLt("posNshared_lt",pos_trk->getNShared(),weight)) { + continue; + } + + //Vertex Quality + if (!_reg_vtx_selectors[region]->passCutLt("chi2unc_lt",vtx->getChi2(),weight)) + continue; + + //Max vtx momentum + if (!_reg_vtx_selectors[region]->passCutLt("maxVtxMom_lt",(ele_mom+pos_mom).Mag(),weight)) + continue; + + //Min vtx momentum + if (!_reg_vtx_selectors[region]->passCutGt("minVtxMom_gt",(ele_mom+pos_mom).Mag(),weight)) + continue; + + //END PRESELECTION CUTS + + //L1 requirement + if (!_reg_vtx_selectors[region]->passCutEq("L1Requirement_eq",(int)(foundL1ele&&foundL1pos),weight)) + continue; + + //L2 requirement + if (!_reg_vtx_selectors[region]->passCutEq("L2Requirement_eq",(int)(foundL2ele&&foundL2pos),weight)) + continue; + + //L1 requirement for positron + if (!_reg_vtx_selectors[region]->passCutEq("L1PosReq_eq",(int)(foundL1pos),weight)) + continue; + + //ESum low cut + if (!_reg_vtx_selectors[region]->passCutLt("eSum_lt",(ele_E+pos_E),weight)) + continue; + + //ESum high cut + if (!_reg_vtx_selectors[region]->passCutGt("eSum_gt",(ele_E+pos_E),weight)) + continue; + + //PSum low cut + if (!_reg_vtx_selectors[region]->passCutLt("pSum_lt",(p_ele.P()+p_pos.P()),weight)) + continue; + + //PSum high cut + if (!_reg_vtx_selectors[region]->passCutGt("pSum_gt",(p_ele.P()+p_pos.P()),weight)) + continue; + + //Require Electron Cluster exists + if (!_reg_vtx_selectors[region]->passCutGt("eleClusE_gt",eleClus.getEnergy(),weight)) + continue; + + + //Require Electron Cluster does NOT exists + if (!_reg_vtx_selectors[region]->passCutLt("eleClusE_lt",eleClus.getEnergy(),weight)) + continue; + + //No shared hits requirement + if (!_reg_vtx_selectors[region]->passCutEq("ele_sharedL0_eq",(int)ele_trk->getSharedLy0(),weight)) + continue; + if (!_reg_vtx_selectors[region]->passCutEq("pos_sharedL0_eq",(int)pos_trk->getSharedLy0(),weight)) + continue; + if (!_reg_vtx_selectors[region]->passCutEq("ele_sharedL1_eq",(int)ele_trk->getSharedLy1(),weight)) + continue; + if (!_reg_vtx_selectors[region]->passCutEq("pos_sharedL1_eq",(int)pos_trk->getSharedLy1(),weight)) + continue; + + //Min vtx Y pos + if (!_reg_vtx_selectors[region]->passCutGt("VtxYPos_gt", vtx->getY(), weight)) + continue; + + //Max vtx Y pos + if (!_reg_vtx_selectors[region]->passCutLt("VtxYPos_lt", vtx->getY(), weight)) + continue; + + //Tracking Volume for positron + if (!_reg_vtx_selectors[region]->passCutGt("volPos_top", p_pos.Py(), weight)) + continue; + + if (!_reg_vtx_selectors[region]->passCutLt("volPos_bot", p_pos.Py(), weight)) + continue; + + if (!_reg_vtx_selectors[region]->passCutLt("deltaZ_lt", std::abs((ele_trk->getZ0()/ele_trk->getTanLambda()) - (pos_trk->getZ0()/pos_trk->getTanLambda())), weight)) + continue; + + //If this is MC check if MCParticle matched to the electron track is from rad or recoil + if(!isData_) + { + //Fill MC plots after all selections + _reg_mc_vtx_histos[region]->FillMCParticles(mcParts_, analysis_); + + //Count the number of hits per part on the ele track + std::map nHits4part; + for(int i =0; i < ele_trk->getMcpHits().size(); i++) + { + int partID = ele_trk->getMcpHits().at(i).second; + if ( nHits4part.find(partID) == nHits4part.end() ) + { + // not found + nHits4part[partID] = 1; + } + else + { + // found + nHits4part[partID]++; + } + } + + //Determine the MC part with the most hits on the track + int maxNHits = 0; + int maxID = 0; + for (std::map::iterator it=nHits4part.begin(); it!=nHits4part.end(); ++it) + { + if(it->second > maxNHits) + { + maxNHits = it->second; + maxID = it->first; + } + } + + //Find the correct mc part and grab mother id + int isRadEle = -999; + int isRecEle = -999; + + + trueEleP.SetXYZ(-999,-999,-999); + truePosP.SetXYZ(-999,-999,-999); + if (mcParts_) { + float trueEleE = -1; + float truePosE = -1; + for(int i = 0; i < mcParts_->size(); i++) + { + int momPDG = mcParts_->at(i)->getMomPDG(); + if(mcParts_->at(i)->getPDG() == 11 && momPDG == isRadPDG_) + { + std::vector lP = mcParts_->at(i)->getMomentum(); + trueEleP.SetXYZ(lP[0],lP[1],lP[2]); + trueEleE = mcParts_->at(i)->getEnergy(); + + } + if(mcParts_->at(i)->getPDG() == -11 && momPDG == isRadPDG_) + { + std::vector lP = mcParts_->at(i)->getMomentum(); + truePosP.SetXYZ(lP[0],lP[1],lP[2]); + truePosE = mcParts_->at(i)->getEnergy(); + + } + if(trueEleP.X() != -999 && truePosP.X() != -999){ + truePsum = trueEleP.Mag() + trueEleP.Mag(); + trueEsum = trueEleE + truePosE; + } + + if(mcParts_->at(i)->getID() != maxID) continue; + //Default isRadPDG = 622 + if(momPDG == isRadPDG_) isRadEle = 1; + if(momPDG == 623) isRecEle = 1; + } + } + double momRatio = recEleP.Mag() / trueEleP.Mag(); + double momAngle = trueEleP.Angle(recEleP) * TMath::RadToDeg(); + if (!_reg_vtx_selectors[region]->passCutLt("momRatio_lt", momRatio, weight)) continue; + if (!_reg_vtx_selectors[region]->passCutGt("momRatio_gt", momRatio, weight)) continue; + if (!_reg_vtx_selectors[region]->passCutLt("momAngle_lt", momAngle, weight)) continue; + + if (!_reg_vtx_selectors[region]->passCutEq("isRadEle_eq", isRadEle, weight)) continue; + if (!_reg_vtx_selectors[region]->passCutEq("isNotRadEle_eq", isRadEle, weight)) continue; + if (!_reg_vtx_selectors[region]->passCutEq("isRecEle_eq", isRecEle, weight)) continue; + } + + goodVtx = vtx; + nGoodVtx++; + goodVtxs.push_back(vtx); + } // selected vertices + + //N selected vertices - this is quite a silly cut to make at the end. But okay. that's how we decided atm. + if (!_reg_vtx_selectors[region]->passCutEq("nVtxs_eq", nGoodVtx, weight)) + continue; + //Move to after N vertices cut (was filled before) + _reg_vtx_histos[region]->Fill1DHisto("n_vertices_h", nGoodVtx, weight); + + //Loop over all selected vertices in the region + for(std::vector::iterator it = goodVtxs.begin(); it != goodVtxs.end(); it++){ + + Vertex* vtx = *it; + + Particle* ele = nullptr; + Particle* pos = nullptr; + + if (!vtx || !_ah->GetParticlesFromVtx(vtx,ele,pos)) + continue; + + CalCluster eleClus = ele->getCluster(); + CalCluster posClus = pos->getCluster(); + + double corr_eleClusterTime = ele->getCluster().getTime() - timeOffset_; + double corr_posClusterTime = pos->getCluster().getTime() - timeOffset_; + + double ele_E = ele->getEnergy(); + double pos_E = pos->getEnergy(); + + //Compute analysis variables here. + Track * ele_trk = (Track*) ele->getTrack().Clone(); + Track * pos_trk = (Track*) pos->getTrack().Clone(); + //Get the shared info - TODO change and improve + + //Track Time Corrections + ele_trk->applyCorrection("track_time",eleTrackTimeBias_); + pos_trk->applyCorrection("track_time", posTrackTimeBias_); + + //Get the layers hit on each track + std::vector ele_hit_layers = ele_trk->getHitLayers(); + int ele_Si0 = 0; + int ele_Si1 = 0; + int ele_lastlayer = 0; + for(int i=0; i pos_hit_layers = pos_trk->getHitLayers(); + int pos_Si0 = 0; + int pos_Si1 = 0; + int pos_lastlayer = 0; + for(int i=0; i vtx_cov = vtx->getCovariance(); + float cxx = vtx_cov.at(0); + float cyx = vtx_cov.at(1); + float cyy = vtx_cov.at(2); + float czx = vtx_cov.at(3); + float czy = vtx_cov.at(4); + float czz = vtx_cov.at(5); + + + //MC Truth hits in first 4 sensors + int L1L2hitCode = 0; //hit code '1111' means truth ax+ster hits in L1_ele, L1_pos, L2_ele, L2_pos + int L1hitCode = 0; //hit code '1111' means truth in L1_ele_ax, L1_ele_ster, L1_pos_ax, L1_pos_ster + int L2hitCode = 0; // hit code '1111' means truth in L2_ele_ax, L2_ele_ster, L2_pos_ax, L2_pos_ster + if(!isData_){ + //Get hit codes. Only sure this works for 2016 KF as is. + utils::get2016KFMCTruthHitCodes(ele_trk, pos_trk, L1L2hitCode, L1hitCode, L2hitCode); + //L1L2 truth hit selection + if (!_reg_vtx_selectors[region]->passCutLt("hitCode_lt",((double)L1L2hitCode)-0.5, weight)) continue; + if (!_reg_vtx_selectors[region]->passCutGt("hitCode_gt",((double)L1L2hitCode)+0.5, weight)) continue; + //Fil hitcodes + _reg_vtx_histos[region]->Fill1DHisto("hitCode_h", L1L2hitCode,weight); + _reg_vtx_histos[region]->Fill1DHisto("L1hitCode_h", L1hitCode,weight); + _reg_vtx_histos[region]->Fill1DHisto("L2hitCode_h", L2hitCode,weight); + } + + //track isolations + //Only calculate isolations if both track L1 and L2 hits exist + bool hasL1ele = false; + bool hasL2ele = false; + _ah->InnermostLayerCheck(ele_trk, hasL1ele, hasL2ele); + + bool hasL1pos = false; + bool hasL2pos = false; + _ah->InnermostLayerCheck(pos_trk, hasL1pos, hasL2pos); + + TVector3 ele_mom; + //ele_mom.SetX(ele->getMomentum()[0]); + //ele_mom.SetY(ele->getMomentum()[1]); + //ele_mom.SetZ(ele->getMomentum()[2]); + ele_mom.SetX(ele_trk->getMomentum()[0]); + ele_mom.SetY(ele_trk->getMomentum()[1]); + ele_mom.SetZ(ele_trk->getMomentum()[2]); + + + TVector3 pos_mom; + //pos_mom.SetX(pos->getMomentum()[0]); + //pos_mom.SetY(pos->getMomentum()[1]); + //pos_mom.SetZ(pos->getMomentum()[2]); + pos_mom.SetX(pos_trk->getMomentum()[0]); + pos_mom.SetY(pos_trk->getMomentum()[1]); + pos_mom.SetZ(pos_trk->getMomentum()[2]); + + double ele_pos_dt = corr_eleClusterTime - corr_posClusterTime; + double psum = ele_mom.Mag()+pos_mom.Mag(); + + //Ele nHits + int ele2dHits = ele_trk->getTrackerHitCount(); + if (!ele_trk->isKalmanTrack()) + ele2dHits*=2; + + //pos nHits + int pos2dHits = pos_trk->getTrackerHitCount(); + + if(ts_ != nullptr) + { + _reg_vtx_histos[region]->Fill2DHisto("trig_count_hh", + ((int)ts_->prescaled.Single_3_Top)+((int)ts_->prescaled.Single_3_Bot), + ((int)ts_->prescaled.Single_2_Top)+((int)ts_->prescaled.Single_2_Bot)); + } + _reg_vtx_histos[region]->Fill1DHisto("n_vtx_h", vtxs_->size()); + + //Add the momenta to the tracks + //ele_trk->setMomentum(ele->getMomentum()[0],ele->getMomentum()[1],ele->getMomentum()[2]); + //pos_trk->setMomentum(pos->getMomentum()[0],pos->getMomentum()[1],pos->getMomentum()[2]); + TVector3 recEleP(ele->getMomentum()[0],ele->getMomentum()[1],ele->getMomentum()[2]); + TLorentzVector p_ele; + p_ele.SetPxPyPzE(ele_trk->getMomentum()[0],ele_trk->getMomentum()[1],ele_trk->getMomentum()[2], ele_E); + TLorentzVector p_pos; + p_pos.SetPxPyPzE(pos_trk->getMomentum()[0],pos_trk->getMomentum()[1],pos_trk->getMomentum()[2], pos_E); + + + _reg_vtx_histos[region]->Fill2DHistograms(vtx,weight); + _reg_vtx_histos[region]->Fill1DVertex(vtx, + ele, + pos, + ele_trk, + pos_trk, + weight); + + _reg_vtx_histos[region]->Fill1DHisto("ele_pos_clusTimeDiff_h", (corr_eleClusterTime - corr_posClusterTime), weight); + _reg_vtx_histos[region]->Fill1DHisto("ele_track_n2dhits_h", ele2dHits, weight); + _reg_vtx_histos[region]->Fill1DHisto("pos_track_n2dhits_h", pos2dHits, weight); + _reg_vtx_histos[region]->Fill1DHisto("vtx_Psum_h", p_ele.P()+p_pos.P(), weight); + _reg_vtx_histos[region]->Fill1DHisto("vtx_Esum_h", eleClus.getEnergy()+posClus.getEnergy(), weight); + _reg_vtx_histos[region]->Fill2DHisto("ele_vtxZ_iso_hh", TMath::Min(ele_trk->getIsolation(0), ele_trk->getIsolation(1)), vtx->getZ(), weight); + _reg_vtx_histos[region]->Fill2DHisto("pos_vtxZ_iso_hh", TMath::Min(pos_trk->getIsolation(0), pos_trk->getIsolation(1)), vtx->getZ(), weight); + _reg_vtx_histos[region]->Fill2DTrack(ele_trk,weight,"ele_"); + _reg_vtx_histos[region]->Fill2DTrack(pos_trk,weight,"pos_"); + _reg_vtx_histos[region]->Fill1DHisto("mcMass622_h",apMass); + _reg_vtx_histos[region]->Fill1DHisto("mcZ622_h",apZ); + _reg_vtx_histos[region]->Fill1DHisto("mcMass625_h",vdMass); + _reg_vtx_histos[region]->Fill1DHisto("mcZ625_h",vdZ); + + + //Just for the selected vertex + if(!isData_) + { + _reg_vtx_histos[region]->Fill2DHisto("vtx_Esum_vs_true_Esum_hh",eleClus.getEnergy()+posClus.getEnergy(), trueEsum, weight); + _reg_vtx_histos[region]->Fill2DHisto("vtx_Psum_vs_true_Psum_hh",p_ele.P()+p_pos.P(), truePsum, weight); + _reg_vtx_histos[region]->Fill1DHisto("true_vtx_psum_h",truePsum,weight); + } + + double reconz = vtx->getZ(); + double ele_trk_z0 = ele_trk->getZ0(); + double ele_trk_z0err = ele_trk->getZ0Err(); + double pos_trk_z0 = pos_trk->getZ0(); + double pos_trk_z0err = pos_trk->getZ0Err(); + + //DeltaZ + double deltaZ = std::abs( (ele_trk_z0/ele_trk->getTanLambda()) - (pos_trk_z0/pos_trk->getTanLambda()) ); + + //Project vertex to target + double vtx_proj_x = -999.9; + double vtx_proj_y = -999.9; + double vtx_proj_x_sig = -999.9; + double vtx_proj_y_sig = -999.9; + double vtx_proj_sig = -999.9; + if(!v0ProjectionFitsCfg_.empty()) + vtx_proj_sig = utils::v0_projection_to_target_significance(v0proj_fits_, evth_->getRunNumber(), + vtx_proj_x, vtx_proj_y, vtx_proj_x_sig, vtx_proj_y_sig, vtx->getX(), vtx->getY(), + reconz, vtx->getP().X(), vtx->getP().Y(), vtx->getP().Z()); + + _reg_vtx_histos[region]->Fill2DHisto("unc_vtx_x_v_unc_vtx_y_hh", vtx->getX(), vtx->getY()); + _reg_vtx_histos[region]->Fill2DHisto("unc_vtx_proj_x_v_unc_vtx_proj_y_hh", vtx_proj_x, vtx_proj_y); + _reg_vtx_histos[region]->Fill2DHisto("unc_vtx_proj_x_y_significance_hh", vtx_proj_x_sig, vtx_proj_y_sig); + _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_vtx_proj_significance_hh", vtx_proj_sig, reconz); + _reg_vtx_histos[region]->Fill2DHisto("vtx_track_reconz_v_Z0err_hh", ele_trk_z0err, reconz); + _reg_vtx_histos[region]->Fill2DHisto("vtx_track_reconz_v_Z0err_hh", pos_trk_z0err, reconz); + _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_z0_hh", ele_trk_z0, reconz); + _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_z0_hh", pos_trk_z0, reconz); + _reg_vtx_histos[region]->Fill2DHisto("vtx_track_recon_z_v_ABSdz0tanlambda_hh", std::abs((ele_trk_z0/ele_trk->getTanLambda()) - (pos_trk_z0/pos_trk->getTanLambda())), reconz); + _reg_vtx_histos[region]->Fill2DHisto("vtx_track_recon_z_v_dz0tanlambda_hh", ((ele_trk_z0/ele_trk->getTanLambda()) - (pos_trk_z0/pos_trk->getTanLambda())), reconz); + + _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_cxx_hh", cxx, reconz); + _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_cyy_hh", cyy, reconz); + _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_czz_hh", czz, reconz); + _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_cyx_hh", cyx, reconz); + _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_czx_hh", czx, reconz); + _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_czy_hh", czy, reconz); + _reg_vtx_histos[region]->Fill1DHisto("cxx_h", cxx); + _reg_vtx_histos[region]->Fill1DHisto("cyy_h", cyy); + _reg_vtx_histos[region]->Fill1DHisto("czz_h", czz); + _reg_vtx_histos[region]->Fill1DHisto("cyx_h", cyx); + _reg_vtx_histos[region]->Fill1DHisto("czx_h", czx); + _reg_vtx_histos[region]->Fill1DHisto("czy_h", czy); + _reg_vtx_histos[region]->Fill2DHisto("vtx_track_recon_z_v_z0tanlambda_hh", ele_trk_z0/ele_trk->getTanLambda(), reconz); + _reg_vtx_histos[region]->Fill2DHisto("vtx_track_recon_z_v_z0tanlambda_hh", pos_trk_z0/pos_trk->getTanLambda(), reconz); + _reg_vtx_histos[region]->Fill2DHisto("ele_clusT_v_ele_trackT_hh", ele_trk->getTrackTime(), corr_eleClusterTime, weight); + _reg_vtx_histos[region]->Fill2DHisto("pos_clusT_v_pos_trackT_hh", pos_trk->getTrackTime(), corr_posClusterTime, weight); + _reg_vtx_histos[region]->Fill2DHisto("ele_track_time_v_P_hh", ele_trk->getP(), ele_trk->getTrackTime(), weight); + _reg_vtx_histos[region]->Fill2DHisto("pos_track_time_v_P_hh", pos_trk->getP(), pos_trk->getTrackTime(), weight); + _reg_vtx_histos[region]->Fill2DHisto("ele_pos_clusTimeDiff_v_pSum_hh",ele_mom.Mag()+pos_mom.Mag(), ele_pos_dt, weight); + _reg_vtx_histos[region]->Fill2DHisto("ele_cluster_energy_v_track_p_hh",ele_trk->getP(), eleClus.getEnergy(), weight); + _reg_vtx_histos[region]->Fill2DHisto("pos_cluster_energy_v_track_p_hh",pos_trk->getP(), posClus.getEnergy(), weight); + _reg_vtx_histos[region]->Fill2DHisto("ele_track_cluster_dt_v_EoverP_hh",eleClus.getEnergy()/ele_trk->getP(), ele_trk->getTrackTime() - corr_eleClusterTime, weight); + _reg_vtx_histos[region]->Fill2DHisto("pos_track_cluster_dt_v_EoverP_hh",posClus.getEnergy()/pos_trk->getP(), pos_trk->getTrackTime() - corr_posClusterTime, weight); + _reg_vtx_histos[region]->Fill2DHisto("ele_track_clus_dt_v_p_hh",ele_trk->getP(), ele_trk->getTrackTime() - corr_eleClusterTime, weight); + _reg_vtx_histos[region]->Fill2DHisto("pos_track_clus_dt_v_p_hh",pos_trk->getP(), pos_trk->getTrackTime() - corr_posClusterTime, weight); + _reg_vtx_histos[region]->Fill2DHisto("ele_z0_vs_pos_z0_hh",ele_trk->getZ0(), pos_trk->getZ0(), weight); + + //chi2 2d plots + _reg_vtx_histos[region]->Fill2DHisto("ele_track_chi2ndf_v_time_hh", ele_trk->getTrackTime(), ele_trk->getChi2Ndf(), weight); + _reg_vtx_histos[region]->Fill2DHisto("ele_track_chi2ndf_v_p_hh", ele_trk->getP(), ele_trk->getChi2Ndf(), weight); + _reg_vtx_histos[region]->Fill2DHisto("ele_track_chi2ndf_v_tanlambda_hh", ele_trk->getTanLambda(), ele_trk->getChi2Ndf(), weight); + _reg_vtx_histos[region]->Fill2DHisto("ele_track_chi2ndf_v_n2dhits_hh", ele2dHits, ele_trk->getChi2Ndf(), weight); + + _reg_vtx_histos[region]->Fill2DHisto("pos_track_chi2ndf_v_time_hh", pos_trk->getTrackTime(), pos_trk->getChi2Ndf(), weight); + _reg_vtx_histos[region]->Fill2DHisto("pos_track_chi2ndf_v_p_hh", pos_trk->getP(), pos_trk->getChi2Ndf(), weight); + _reg_vtx_histos[region]->Fill2DHisto("pos_track_chi2ndf_v_tanlambda_hh", pos_trk->getTanLambda(), pos_trk->getChi2Ndf(), weight); + _reg_vtx_histos[region]->Fill2DHisto("pos_track_chi2ndf_v_n2dhits_hh", pos2dHits, pos_trk->getChi2Ndf(), weight); + + _reg_vtx_histos[region]->Fill2DHisto("ele_track_p_v_tanlambda_hh", ele_trk->getTanLambda(), ele_trk->getP(), weight); + _reg_vtx_histos[region]->Fill2DHisto("pos_track_p_v_tanlambda_hh", pos_trk->getTanLambda(), pos_trk->getP(), weight); + + + //1d histos + _reg_vtx_histos[region]->Fill1DHisto("ele_track_clus_dt_h", ele_trk->getTrackTime() - corr_eleClusterTime, weight); + _reg_vtx_histos[region]->Fill1DHisto("pos_track_clus_dt_h", pos_trk->getTrackTime() - corr_posClusterTime, weight); + + + //TODO put this in the Vertex! + TVector3 vtxPosSvt; + vtxPosSvt.SetX(vtx->getX()); + vtxPosSvt.SetY(vtx->getY()); + vtxPosSvt.SetZ(vtx->getZ()); + vtxPosSvt.RotateY(-0.0305); + + //Just for the selected vertex + if (makeFlatTuple_){ + if(!isData_){ + _reg_tuples[region]->setVariableValue("ap_true_vtx_z", apZ); + _reg_tuples[region]->setVariableValue("ap_true_vtx_mass", apMass); + _reg_tuples[region]->setVariableValue("ap_true_vtx_energy", apEnergy); + _reg_tuples[region]->setVariableValue("vd_true_vtx_z", vdZ); + _reg_tuples[region]->setVariableValue("vd_true_vtx_mass", vdMass); + _reg_tuples[region]->setVariableValue("vd_true_vtx_energy", vdEnergy); + _reg_tuples[region]->setVariableValue("hitCode", float(L1L2hitCode)); + _reg_tuples[region]->setVariableValue("L1hitCode", float(L1hitCode)); + _reg_tuples[region]->setVariableValue("L2hitCode", float(L2hitCode)); + } + + _reg_tuples[region]->setVariableValue("unc_vtx_mass", vtx->getInvMass()); + _reg_tuples[region]->setVariableValue("unc_vtx_z" , vtxPosSvt.Z()); + _reg_tuples[region]->setVariableValue("unc_vtx_chi2", vtx->getChi2()); + _reg_tuples[region]->setVariableValue("unc_vtx_psum", p_ele.P()+p_pos.P()); + _reg_tuples[region]->setVariableValue("unc_vtx_px", vtx->getP().X()); + _reg_tuples[region]->setVariableValue("unc_vtx_py", vtx->getP().Y()); + _reg_tuples[region]->setVariableValue("unc_vtx_pz", vtx->getP().Z()); + _reg_tuples[region]->setVariableValue("unc_vtx_x", vtx->getX()); + _reg_tuples[region]->setVariableValue("unc_vtx_y", vtx->getY()); + _reg_tuples[region]->setVariableValue("unc_vtx_proj_x", vtx_proj_x); + _reg_tuples[region]->setVariableValue("unc_vtx_proj_y", vtx_proj_y); + _reg_tuples[region]->setVariableValue("unc_vtx_proj_x_sig", vtx_proj_x_sig); + _reg_tuples[region]->setVariableValue("unc_vtx_proj_y_sig", vtx_proj_y_sig); + _reg_tuples[region]->setVariableValue("unc_vtx_proj_sig", vtx_proj_sig); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_pos_clust_dt", corr_eleClusterTime - corr_posClusterTime); + + _reg_tuples[region]->setVariableValue("unc_vtx_cxx", cxx); + _reg_tuples[region]->setVariableValue("unc_vtx_cyy", cyy); + _reg_tuples[region]->setVariableValue("unc_vtx_czz", czz); + _reg_tuples[region]->setVariableValue("unc_vtx_cyx", cyx); + _reg_tuples[region]->setVariableValue("unc_vtx_czy", czy); + _reg_tuples[region]->setVariableValue("unc_vtx_czx", czx); + _reg_tuples[region]->setVariableValue("unc_vtx_deltaZ", deltaZ); + + //track vars + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_p", ele_trk->getP()); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_t", ele_trk->getTrackTime()); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_d0", ele_trk->getD0()); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_phi0", ele_trk->getPhi()); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_omega", ele_trk->getOmega()); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_tanLambda", ele_trk->getTanLambda()); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_z0", ele_trk->getZ0()); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_chi2ndf", ele_trk->getChi2Ndf()); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_clust_dt", ele_trk->getTrackTime() - corr_eleClusterTime); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_z0Err",ele_trk->getZ0Err()); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_d0Err", ele_trk->getD0Err()); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_tanLambdaErr", ele_trk->getTanLambdaErr()); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_PhiErr", ele_trk->getPhiErr()); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_OmegaErr", ele_trk->getOmegaErr()); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_nhits",ele2dHits); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_lastlayer",ele_lastlayer); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_si0",ele_Si0); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_si1",ele_Si1); + + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_p", pos_trk->getP()); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_t", pos_trk->getTrackTime()); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_d0", pos_trk->getD0()); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_phi0", pos_trk->getPhi()); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_omega", pos_trk->getOmega()); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_tanLambda", pos_trk->getTanLambda()); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_z0", pos_trk->getZ0()); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_chi2ndf", pos_trk->getChi2Ndf()); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_clust_dt", pos_trk->getTrackTime() - corr_posClusterTime); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_z0Err",pos_trk->getZ0Err()); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_d0Err", pos_trk->getD0Err()); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_tanLambdaErr", pos_trk->getTanLambdaErr()); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_PhiErr", pos_trk->getPhiErr()); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_OmegaErr", pos_trk->getOmegaErr()); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_nhits",pos2dHits); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_lastlayer",pos_lastlayer); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_si0",pos_Si0); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_si1",pos_Si1); + + //clust vars + _reg_tuples[region]->setVariableValue("unc_vtx_ele_clust_E", eleClus.getEnergy()); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_clust_x", eleClus.getPosition().at(0)); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_clust_corr_t",corr_eleClusterTime); + + _reg_tuples[region]->setVariableValue("unc_vtx_pos_clust_E", posClus.getEnergy()); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_clust_x", posClus.getPosition().at(0)); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_clust_corr_t",corr_posClusterTime); + _reg_tuples[region]->setVariableValue("run_number", evth_->getRunNumber()); + + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_ecal_x", ele_trk->getPositionAtEcal().at(0)); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_ecal_y", ele_trk->getPositionAtEcal().at(1)); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_z", ele_trk->getPosition().at(2)); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_ecal_x", pos_trk->getPositionAtEcal().at(0)); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_ecal_y", pos_trk->getPositionAtEcal().at(1)); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_z", pos_trk->getPosition().at(2)); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_px", ele_trk->getMomentum().at(0)); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_py", ele_trk->getMomentum().at(1)); + _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_pz", ele_trk->getMomentum().at(2)); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_px", pos_trk->getMomentum().at(0)); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_py", pos_trk->getMomentum().at(1)); + _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_pz", pos_trk->getMomentum().at(2)); + + _reg_tuples[region]->fill(); + } + } + + }// regions + + return true; +} + +void NewVertexAnaProcessor::finalize() { + + //TODO clean this up a little. + outF_->cd(); + _vtx_histos->saveHistos(outF_,_vtx_histos->getName()); + outF_->cd(_vtx_histos->getName().c_str()); + vtxSelector->getCutFlowHisto()->Write(); + + outF_->cd(); + if(!isData_) + _mc_vtx_histos->saveHistos(outF_, _mc_vtx_histos->getName()); + //delete histos; + //histos = nullptr; + + + for (reg_it it = _reg_vtx_histos.begin(); it!=_reg_vtx_histos.end(); ++it) { + std::string dirName = anaName_+"_"+it->first; + (it->second)->saveHistos(outF_,dirName); + outF_->cd(dirName.c_str()); + _reg_vtx_selectors[it->first]->getCutFlowHisto()->Write(); + //Save tuples + if (makeFlatTuple_) + _reg_tuples[it->first]->writeTree(); + + } + + if(!isData_){ + for (reg_mc_it it = _reg_mc_vtx_histos.begin(); it!=_reg_mc_vtx_histos.end(); ++it) { + std::string dirName = anaName_+"_mc_"+it->first; + (it->second)->saveHistos(outF_,dirName); + outF_->cd(dirName.c_str()); + } + } + + outF_->Close(); + +} + +DECLARE_PROCESSOR(NewVertexAnaProcessor); diff --git a/processors/src/SimpZBiOptimizationProcessor.cxx b/processors/src/SimpZBiOptimizationProcessor.cxx index a4b546ee3..d9b269baf 100644 --- a/processors/src/SimpZBiOptimizationProcessor.cxx +++ b/processors/src/SimpZBiOptimizationProcessor.cxx @@ -61,9 +61,6 @@ void SimpZBiOptimizationProcessor::configure(const ParameterSet& parameters) { new_variables_ = parameters.getVString("add_new_variables", new_variables_); new_variable_params_ = parameters.getVDouble("new_variable_params", new_variable_params_); - //Dev - testSpecialCut_ = parameters.getInteger("testSpecialCut", testSpecialCut_); - } catch (std::runtime_error& error) { @@ -74,6 +71,7 @@ void SimpZBiOptimizationProcessor::configure(const ParameterSet& parameters) { //USE JSON FILE TO LOAD IN NEW VARIABLES AND VARIABLE CONFIGURATIONS void SimpZBiOptimizationProcessor::addNewVariables(SimpAnaTTree* MTT, std::string variable, double param){ + std::cout << "[SimpZBiOptimizationProcessor]::addNewVariable " << variable << " with param " << param << std::endl; if(variable == "unc_vtx_ele_zalpha") MTT->addVariable_unc_vtx_ele_zalpha(param); if(variable == "unc_vtx_pos_zalpha") @@ -117,9 +115,12 @@ void SimpZBiOptimizationProcessor::addNewVariables(SimpAnaTTree* MTT, std::strin if(variable == "unc_vtx_abs_delta_z0tanlambda") MTT->addVariable_unc_vtx_abs_delta_z0tanlambda(); - else - std::cout << "[SimpZBiOptimization]::ERROR::NEW VARIABLE " << variable << " IS NOT DEFINED IN SimpAnaTTree.cxx" - <addVariable_unc_vtx_proj_significance(); + + //else + // std::cout << "[SimpZBiOptimization]::ERROR::NEW VARIABLE " << variable << " IS NOT DEFINED IN SimpAnaTTree.cxx" + // < histos, SimpAnaTTree* MTT){ @@ -131,62 +132,98 @@ void SimpZBiOptimizationProcessor::fillEventHistograms(std::shared_ptrFill1DHisto(var+"_h", MTT->getValue(var)); } - //always fill + //Impact Parameter histos->Fill2DHisto("z0_v_recon_z_hh",MTT->getValue("unc_vtx_z"),MTT->getValue("unc_vtx_ele_track_z0")); histos->Fill2DHisto("z0_v_recon_z_hh",MTT->getValue("unc_vtx_z"),MTT->getValue("unc_vtx_pos_track_z0")); - histos->Fill2DHisto("ele_track_z0_v_pos_track_z0_hh",MTT->getValue("unc_vtx_pos_track_z0"), - MTT->getValue("unc_vtx_ele_track_z0")); + + //Inv mass histos->Fill2DHisto("vtx_InvM_vtx_z_hh",MTT->getValue("unc_vtx_mass")*1000.0, MTT->getValue("unc_vtx_z")); + //track z0 + histos->Fill2DHisto("ele_track_z0_v_pos_track_z0_hh", MTT->getValue("unc_vtx_ele_track_z0"), + MTT->getValue("unc_vtx_pos_track_z0")); - //Investigating new variables + //Zalpha if(MTT->variableExists("unc_vtx_ele_zalpha")){ - histos->Fill2DHisto("z0_v_recon_zalpha_hh",MTT->getValue("unc_vtx_ele_zalpha"), + histos->Fill2DHisto("z0_v_unc_vtx_zalpha_hh",MTT->getValue("unc_vtx_ele_zalpha"), MTT->getValue("unc_vtx_ele_track_z0")); } if(MTT->variableExists("unc_vtx_pos_zalpha")){ - histos->Fill2DHisto("z0_v_recon_zalpha_hh",MTT->getValue("unc_vtx_pos_zalpha"), + histos->Fill2DHisto("z0_v_unc_vtx_zalpha_hh",MTT->getValue("unc_vtx_pos_zalpha"), MTT->getValue("unc_vtx_pos_track_z0")); } - - //isolation cut - if(MTT->variableExists("unc_vtx_ele_iso_z0err")){ - histos->Fill2DHisto("recon_z_v_iso_z0err_hh", MTT->getValue("unc_vtx_ele_iso_z0err"),MTT->getValue("unc_vtx_z")); - histos->Fill2DHisto("recon_z_v_iso_z0err_hh", MTT->getValue("unc_vtx_pos_iso_z0err"),MTT->getValue("unc_vtx_z")); + if(MTT->variableExists("unc_vtx_zalpha_max")){ + histos->Fill2DHisto("recon_z_v_unc_vtx_zalpha_max_hh",MTT->getValue("unc_vtx_zalpha_max"), + MTT->getValue("unc_vtx_z")); } - if(MTT->variableExists("unc_vtx_ele_z0_z0err")){ - histos->Fill2DHisto("recon_z_v_z0_z0err_hh", MTT->getValue("unc_vtx_ele_z0_z0err"),MTT->getValue("unc_vtx_z")); - histos->Fill2DHisto("recon_z_v_z0_z0err_hh", MTT->getValue("unc_vtx_pos_z0_z0err"),MTT->getValue("unc_vtx_z")); + + //v0 projection + if(MTT->variableExists("unc_vtx_proj_x")){ + histos->Fill2DHisto("unc_vtx_proj_x_v_unc_vtx_proj_y_hh", + MTT->getValue("unc_vtx_proj_x"), MTT->getValue("unc_vtx_proj_y")); } - if(MTT->variableExists("unc_vtx_ele_z0_z0err")){ - histos->Fill2DHisto("recon_z_v_z0_z0err_hh", MTT->getValue("unc_vtx_ele_z0_z0err"),MTT->getValue("unc_vtx_z")); - histos->Fill2DHisto("recon_z_v_z0_z0err_hh", MTT->getValue("unc_vtx_pos_z0_z0err"),MTT->getValue("unc_vtx_z")); + if(MTT->variableExists("unc_vtx_proj_x_sig")){ + histos->Fill2DHisto("unc_vtx_proj_x_y_significance_hh", + MTT->getValue("unc_vtx_proj_x_sig"), MTT->getValue("unc_vtx_proj_y_sig")); } - if(MTT->variableExists("unc_vtx_ele_isolation_cut")){ - histos->Fill2DHisto("recon_z_v_isolation_cut_hh", MTT->getValue("unc_vtx_ele_isolation_cut"),MTT->getValue("unc_vtx_z")); - histos->Fill2DHisto("recon_z_v_isolation_cut_hh", MTT->getValue("unc_vtx_pos_isolation_cut"),MTT->getValue("unc_vtx_z")); + if(MTT->variableExists("unc_vtx_proj_sig")){ + histos->Fill2DHisto("recon_z_v_proj_sig_hh", + MTT->getValue("unc_vtx_proj_sig"), MTT->getValue("unc_vtx_z")); } + //Vertex Errors if(MTT->variableExists("unc_vtx_cxx")){ histos->Fill2DHisto("recon_z_v_cxx_hh", MTT->getValue("unc_vtx_cxx"),MTT->getValue("unc_vtx_z")); histos->Fill2DHisto("recon_z_v_cyy_hh", MTT->getValue("unc_vtx_cyy"),MTT->getValue("unc_vtx_z")); histos->Fill2DHisto("recon_z_v_czz_hh", MTT->getValue("unc_vtx_czz"),MTT->getValue("unc_vtx_z")); + histos->Fill2DHisto("recon_z_v_czx_hh", MTT->getValue("unc_vtx_czx"),MTT->getValue("unc_vtx_z")); + histos->Fill2DHisto("recon_z_v_czy_hh", MTT->getValue("unc_vtx_czy"),MTT->getValue("unc_vtx_z")); + histos->Fill2DHisto("recon_z_v_cyx_hh", MTT->getValue("unc_vtx_cyx"),MTT->getValue("unc_vtx_z")); } - //Z0TanLambda - histos->Fill2DHisto("recon_z_v_z0tanlambda_hh", MTT->getValue("unc_vtx_ele_track_z0")/MTT->getValue("unc_vtx_ele_track_tanLambda"),MTT->getValue("unc_vtx_z")); - histos->Fill2DHisto("recon_z_v_z0tanlambda_hh", MTT->getValue("unc_vtx_pos_track_z0")/MTT->getValue("unc_vtx_pos_track_tanLambda"),MTT->getValue("unc_vtx_z")); - if(MTT->variableExists("unc_vtx_ele_z0tanlambda_right")){ - histos->Fill2DHisto("recon_z_v_z0tanlambda_right_hh", MTT->getValue("unc_vtx_ele_z0tanlambda_right"),MTT->getValue("unc_vtx_z")); - histos->Fill2DHisto("recon_z_v_z0tanlambda_right_hh", MTT->getValue("unc_vtx_pos_z0tanlambda_right"),MTT->getValue("unc_vtx_z")); - } - if(MTT->variableExists("unc_vtx_ele_z0tanlambda_left")){ - histos->Fill2DHisto("recon_z_v_z0tanlambda_left_hh", MTT->getValue("unc_vtx_ele_z0tanlambda_left"),MTT->getValue("unc_vtx_z")); - histos->Fill2DHisto("recon_z_v_z0tanlambda_left_hh", MTT->getValue("unc_vtx_pos_z0tanlambda_left"),MTT->getValue("unc_vtx_z")); - } - - if(MTT->variableExists("unc_vtx_abs_delta_z0tanlambda")){ - histos->Fill2DHisto("recon_z_v_abs_delta_z0tanlambda_hh", MTT->getValue("unc_vtx_abs_delta_z0tanlambda"),MTT->getValue("unc_vtx_z")); + histos->Fill2DHisto("recon_z_v_z0tanlambda_hh", + MTT->getValue("unc_vtx_ele_track_z0")/MTT->getValue("unc_vtx_ele_track_tanLambda"),MTT->getValue("unc_vtx_z")); + histos->Fill2DHisto("recon_z_v_z0tanlambda_hh", + MTT->getValue("unc_vtx_pos_track_z0")/MTT->getValue("unc_vtx_pos_track_tanLambda"),MTT->getValue("unc_vtx_z")); + if(MTT->variableExists("unc_vtx_deltaZ")){ + histos->Fill2DHisto("recon_z_v_unc_vtx_deltaZ_hh", MTT->getValue("unc_vtx_deltaZ"),MTT->getValue("unc_vtx_z")); } + //z0 error + histos->Fill2DHisto("recon_z_v_Z0err_hh", + MTT->getValue("unc_vtx_ele_track_z0Err"),MTT->getValue("unc_vtx_z")); + histos->Fill2DHisto("recon_z_v_Z0err_hh", + MTT->getValue("unc_vtx_pos_track_z0Err"),MTT->getValue("unc_vtx_z")); + //track time + histos->Fill2DHisto("recon_z_v_track_t_hh", MTT->getValue("unc_vtx_ele_track_t"),MTT->getValue("unc_vtx_z")); + histos->Fill2DHisto("recon_z_v_track_t_hh", MTT->getValue("unc_vtx_pos_track_t"),MTT->getValue("unc_vtx_z")); + + //Track parameters + histos->Fill2DHisto("recon_z_v_ele_track_d0_hh",MTT->getValue("unc_vtx_ele_track_d0"), MTT->getValue("unc_vtx_z")); + histos->Fill2DHisto("recon_z_v_pos_track_d0_hh",MTT->getValue("unc_vtx_pos_track_d0"), MTT->getValue("unc_vtx_z")); + histos->Fill2DHisto("recon_z_v_ele_track_phi0_hh",MTT->getValue("unc_vtx_ele_track_phi0"), MTT->getValue("unc_vtx_z")); + histos->Fill2DHisto("recon_z_v_pos_track_phi0_hh",MTT->getValue("unc_vtx_pos_track_phi0"), MTT->getValue("unc_vtx_z")); + histos->Fill2DHisto("recon_z_v_ele_track_px_hh",MTT->getValue("unc_vtx_ele_track_px"), MTT->getValue("unc_vtx_z")); + histos->Fill2DHisto("recon_z_v_pos_track_px_hh",MTT->getValue("unc_vtx_pos_track_px"), MTT->getValue("unc_vtx_z")); + histos->Fill2DHisto("recon_z_v_ele_track_py_hh",MTT->getValue("unc_vtx_ele_track_py"), MTT->getValue("unc_vtx_z")); + histos->Fill2DHisto("recon_z_v_pos_track_py_hh",MTT->getValue("unc_vtx_pos_track_py"), MTT->getValue("unc_vtx_z")); + histos->Fill2DHisto("recon_z_v_ele_track_pz_hh",MTT->getValue("unc_vtx_ele_track_pz"), MTT->getValue("unc_vtx_z")); + histos->Fill2DHisto("recon_z_v_pos_track_pz_hh",MTT->getValue("unc_vtx_pos_track_pz"), MTT->getValue("unc_vtx_z")); + histos->Fill2DHisto("recon_z_v_ele_track_nhits_hh",MTT->getValue("unc_vtx_ele_track_nhits"), MTT->getValue("unc_vtx_z")); + histos->Fill2DHisto("recon_z_v_pos_track_nhits_hh",MTT->getValue("unc_vtx_pos_track_nhits"), MTT->getValue("unc_vtx_z")); + + //track params vs params + histos->Fill2DHisto("ele_tanlambda_vs_phi0_hh",MTT->getValue("unc_vtx_ele_track_phi0"), + MTT->getValue("unc_vtx_ele_track_tanLambda")); + histos->Fill2DHisto("pos_tanlambda_vs_phi0_hh",MTT->getValue("unc_vtx_pos_track_phi0"), + MTT->getValue("unc_vtx_pos_track_tanLambda")); + histos->Fill2DHisto("ele_cluster_energy_v_track_p_hh",MTT->getValue("unc_vtx_ele_track_p"), + MTT->getValue("unc_vtx_ele_clust_E")); + histos->Fill2DHisto("pos_cluster_energy_v_track_p_hh",MTT->getValue("unc_vtx_pos_track_p"), + MTT->getValue("unc_vtx_pos_clust_E")); + histos->Fill2DHisto("ele_z0_vs_tanlambda_hh",MTT->getValue("unc_vtx_ele_track_tanLambda"), + MTT->getValue("unc_vtx_ele_track_z0")); + histos->Fill2DHisto("pos_z0_vs_tanlambda_hh",MTT->getValue("unc_vtx_pos_track_tanLambda"), + MTT->getValue("unc_vtx_pos_track_z0")); + } double SimpZBiOptimizationProcessor::countControlRegionBackgroundRate(std::string inFilename, std::string tree_name, @@ -262,12 +299,16 @@ void SimpZBiOptimizationProcessor::initialize(std::string inFilename, std::strin int param_idx = std::distance(new_variables_.begin(), it); std::cout << "[SimpZBiOptimization]::Attempting to add new variable " << *it << " with parameter " << new_variable_params_.at(param_idx) << std::endl; + //signalMTT_->addVariable(*it, new_variable_params_.at(param_idx)); + //bkgMTT_->addVariable(*it, new_variable_params_.at(param_idx)); addNewVariables(signalMTT_, *it, new_variable_params_.at(param_idx)); addNewVariables(bkgMTT_, *it, new_variable_params_.at(param_idx)); } //Finalize Initialization of New Mutable Tuples std::cout << "[SimpZBiOptimization]::Finalizing Initialization of New Mutable Tuples" << std::endl; + signalMTT_->shiftVariable("unc_vtx_ele_track_z0", -0.07); + signalMTT_->shiftVariable("unc_vtx_pos_track_z0", -0.07); signalMTT_->Fill(); bkgMTT_->Fill(); @@ -365,7 +406,7 @@ bool SimpZBiOptimizationProcessor::process(){ std::cout << "max iteration: " << max_iteration_ << std::endl; //Iteratively cut n% of the signal distribution for a given Test Cut variable - for(int iteration = 0; iteration < max_iteration_; iteration ++){ + for(int iteration = 0; iteration < max_iteration_+1; iteration ++){ double cutSignal = (double)iteration*step_size_*100.0; cutSignal = round(cutSignal); if(debug_) std::cout << "## ITERATION " << iteration << " ##" << std::endl; @@ -404,6 +445,12 @@ bool SimpZBiOptimizationProcessor::process(){ //Fill Signal variable distributions fillEventHistograms(signalHistos_, signalMTT_); } + if(iteration == max_iteration_){ + //Write iteration histos + signalHistos_->writeHistos(outFile_,"signal_pct_sig_cut_"+std::to_string(cutSignal)); + bkgHistos_->writeHistos(outFile_,"background_pct_sig_cut_"+std::to_string(cutSignal)); + break; + } //Initialize Signal Integrals if(iteration == 0){ @@ -515,8 +562,9 @@ bool SimpZBiOptimizationProcessor::process(){ //Build Background Model, used to estimate nbkg in Signal Region if(debug_) std::cout << "Build Background Model" << std::endl; - TF1* bkg_model = (TF1*)testCutHistos_->fitExponentialTail("background_zVtx_"+cutname, 100.0*background_sf_); - //TF1* bkg_model = (TF1*)testCutHistos_->fitExponentialPlusConst("background_zVtx_"+cutname, 1000.0); + double start_fit = 500.0*background_sf_; + TF1* bkg_model = (TF1*)testCutHistos_->fitExponentialTail("background_zVtx_"+cutname, start_fit); + //TF1* bkg_model = (TF1*)testCutHistos_->fitExponentialPlusConst("background_zVtx_"+cutname, start_fit); if(debug_) std::cout << "END Build Background Model" << std::endl; //Get signal unc_vtx_z vs true_vtx_z @@ -566,18 +614,17 @@ bool SimpZBiOptimizationProcessor::process(){ //Find maximum position of Zcut --> ZBi calculation requires non-zero background //Start the Zcut position at the target - double max_zcut = -4.0; - double zcut_step = 0.5; + double zcut_step = 0.1; TH1F* bkg_zVtx_h = (TH1F*)testCutHistos_->get1dHisto("testCutHistos_background_zVtx_"+cutname+"_h"); - //double endIntegral = bkg_zVtx_h->GetBinCenter(bkg_zVtx_h->FindLastBinAbove(0.0)); - double endIntegral = 100.0; + double max_zcut = bkg_model->GetXmin(); + double endIntegral = bkg_zVtx_h->GetBinLowEdge(bkg_zVtx_h->FindLastBinAbove(0.0)) + bkg_zVtx_h->GetBinWidth(1); + //double endIntegral = 100.0 double testIntegral = bkg_model->Integral(max_zcut, endIntegral); - if(debug_) std::cout << "Background between " << max_zcut << "and end of histo is " << - testIntegral << std::endl; + if(debug_) std::cout << "Background between " << max_zcut << "and end of histo is " << testIntegral << std::endl; while(testIntegral > min_ztail_events_){ max_zcut = max_zcut+zcut_step; testIntegral = bkg_model->Integral(max_zcut, endIntegral); - if(testIntegral == 0.0){ + if(testIntegral < min_ztail_events_){ max_zcut = max_zcut-zcut_step; testIntegral = bkg_model->Integral(max_zcut, endIntegral); break; @@ -586,7 +633,7 @@ bool SimpZBiOptimizationProcessor::process(){ if(debug_) std::cout << "Maximum Zcut: " << max_zcut << " gives " << testIntegral << " background events" << std::endl; //If configuration does not specify scanning Zcut values, use single Zcut position at maximum position. - double min_zcut = 10.0; + double min_zcut = bkg_model->GetXmin(); if(!scan_zcut_) min_zcut = max_zcut; std::cout << "Minimum Zcut position: " << min_zcut << std::endl; @@ -617,8 +664,8 @@ bool SimpZBiOptimizationProcessor::process(){ if(debug_) std::cout << "Scanning zcut position" << std::endl; for(double zcut = min_zcut; zcut < (max_zcut+zcut_step); zcut = zcut+zcut_step){ double Nbkg = bkg_model->Integral(zcut,endIntegral); - std::cout << "zcut: " << zcut << std::endl; - std::cout << "Nbkg: " << Nbkg << std::endl; + //std::cout << "zcut: " << zcut << std::endl; + //std::cout << "Nbkg: " << Nbkg << std::endl; //Get the Signal truth vertex z distribution beyond the reconstructed vertex Zcut TH1F* true_vtx_z_h = (TH1F*)vtx_z_hh->ProjectionY((std::to_string(zcut)+"_"+cutname+"_"+"true_vtx_z_projy").c_str(),vtx_z_hh->GetXaxis()->FindBin(zcut)+1,vtx_z_hh->GetXaxis()->GetNbins(),""); @@ -632,7 +679,7 @@ bool SimpZBiOptimizationProcessor::process(){ } //Get Signal Selection Efficiency, as a function of truth vertex Z, F(z) - if(debug_) std::cout << "Get Signal Selection Efficiency" << std::endl; + //if(debug_) std::cout << "Get Signal Selection Efficiency" << std::endl; TEfficiency* effCalc_h = new TEfficiency(*signalSelZ_h, *signalSimZ_h_); //if(zcut == min_zcut){ // outFile_->cd(("testCuts_pct_sig_cut_"+std::to_string(cutSignal)).c_str()); @@ -644,11 +691,11 @@ bool SimpZBiOptimizationProcessor::process(){ //Calculate expected signal for Neutral Dark Vector "rho" double nSigRho = simpEqs_->expectedSignalCalculation(signal_mass_, - eps, true, false, E_Vd_, effCalc_h, dNdm_, radFrac_, radAcc_, -4.3, zcut); + eps, true, E_Vd_, effCalc_h, dNdm_, radFrac_, radAcc_, -4.3, zcut); //Calculate expected signal for Neutral Dark Vector "rho" double nSigPhi = simpEqs_->expectedSignalCalculation(signal_mass_, - eps, false, true, E_Vd_, effCalc_h, dNdm_, radFrac_, radAcc_, -4.3, zcut); + eps, false, E_Vd_, effCalc_h, dNdm_, radFrac_, radAcc_, -4.3, zcut); /* //Calculate expected signal for Neutral Dark Vector "rho" @@ -662,14 +709,14 @@ bool SimpZBiOptimizationProcessor::process(){ */ double Nsig = nSigRho + nSigPhi; - if(debug_){ - std::cout << "nSigRho: " << nSigRho << std::endl; - std::cout << "nSigPhi: " << nSigPhi << std::endl; - std::cout << "Nsig: " << Nsig << std::endl; - } + //if(debug_){ + // std::cout << "nSigRho: " << nSigRho << std::endl; + // std::cout << "nSigPhi: " << nSigPhi << std::endl; + // std::cout << "Nsig: " << Nsig << std::endl; + //} Nsig = Nsig*signal_sf_; - if(debug_) std::cout << "Nsig after scale factor: " << Nsig << std::endl; + //if(debug_) std::cout << "Nsig after scale factor: " << Nsig << std::endl; //CLEAR POINTERS @@ -687,6 +734,7 @@ bool SimpZBiOptimizationProcessor::process(){ std::cout << "ZBi before rounding: " << ZBi << std::endl; ZBi = round(ZBi); + /* std::cout << "[SimpZBiOptimization]::Iteration Results:" << std::endl; std::cout << "Zcut = " << zcut << std::endl; std::cout << "Nsig = " << Nsig << std::endl; @@ -694,6 +742,7 @@ bool SimpZBiOptimizationProcessor::process(){ std::cout << "n_on: " << n_on << std::endl; std::cout << "n_off: " << n_off << std::endl; std::cout << "ZBi: " << ZBi << std::endl; + */ //Update Test Cut with best scan values if(ZBi > best_scan_zbi){ @@ -825,15 +874,6 @@ bool SimpZBiOptimizationProcessor::failPersistentCuts(SimpAnaTTree* MTT){ } } - //Testing 2016 impact parameter cut 04/10/2023 - //THIS SHOULD EVENTUALLY BE REMOVED - if(testSpecialCut_){ - //if(!MTT->impactParameterCut2016Canonical(signal_mass_)) - // failCuts = true; - if(!MTT->testImpactParameterCut()) - failCuts = true; - } - return failCuts; } diff --git a/processors/src/TrackingAnaProcessor.cxx b/processors/src/TrackingAnaProcessor.cxx index 105de11fa..d641a706b 100644 --- a/processors/src/TrackingAnaProcessor.cxx +++ b/processors/src/TrackingAnaProcessor.cxx @@ -1,6 +1,7 @@ #include "TrackingAnaProcessor.h" #include #include "utilities.h" +#include "AnaHelpers.h" TrackingAnaProcessor::TrackingAnaProcessor(const std::string& name, Process& process) : Processor(name, process) { @@ -20,16 +21,21 @@ void TrackingAnaProcessor::configure(const ParameterSet& parameters) { doTruth_ = (bool) parameters.getInteger("doTruth",doTruth_); truthHistCfgFilename_ = parameters.getString("truthHistCfg",truthHistCfgFilename_); selectionCfg_ = parameters.getString("selectionjson",selectionCfg_); + isData_ = parameters.getInteger("isData",isData_); + ecalCollName_ = parameters.getString("ecalCollName",ecalCollName_); + regionSelections_ = parameters.getVString("regionDefinitions",regionSelections_); } catch (std::runtime_error& error) { std::cout << error.what() << std::endl; } + if (!isData_) + time_offset_ = 5.; } void TrackingAnaProcessor::initialize(TTree* tree) { - + //Init histos trkHistos_ = new TrackHistos(trkCollName_); trkHistos_->loadHistoConfig(histCfgFilename_); @@ -45,7 +51,7 @@ void TrackingAnaProcessor::initialize(TTree* tree) { } if (doTruth_) { - truthHistos_ = new TrackHistos(trkCollName_+"_truthComparison"); + truthHistos_ = new TrackHistos(trkCollName_+"_truthComparison"); truthHistos_->loadHistoConfig(histCfgFilename_); truthHistos_->DefineHistos(); truthHistos_->loadHistoConfig(truthHistCfgFilename_); @@ -55,13 +61,85 @@ void TrackingAnaProcessor::initialize(TTree* tree) { //tree->SetBranchAddress(truthCollName_.c_str(),&truth_tracks_,&btruth_tracks_); } + // Setup track selections plots + for (unsigned int i_reg = 0; + i_reg < regionSelections_.size(); + i_reg++) { + std::string regname = AnaHelpers::getFileName(regionSelections_[i_reg],false); + std::cout<< "Setting up region "<< regname<(regname, regionSelections_[i_reg]); + reg_selectors_[regname]->setDebug(false); + reg_selectors_[regname]->LoadSelection(); + + reg_histos_[regname] = std::make_shared(regname); + reg_histos_[regname]->loadHistoConfig(histCfgFilename_); + reg_histos_[regname]->doTrackComparisonPlots(false); + reg_histos_[regname]->DefineTrkHitHistos(); + + regions_.push_back(regname); + + } + + + + + //Get event header information for trigger + tree->SetBranchAddress("EventHeader", &evth_ , &bevth_); + + //Get cluster information for FEEs + if (!ecalCollName_.empty()) + tree->SetBranchAddress(ecalCollName_.c_str(),&ecal_, &becal_); + } bool TrackingAnaProcessor::process(IEvent* ievent) { - + double weight = 1.; // Loop over all the LCIO Tracks and add them to the HPS event. int n_sel_tracks = 0; + + + //Trigger requirements - Singles 0 and 1. + //TODO use cutFlow + if (isData_ && (!evth_->isSingle0Trigger() && !evth_->isSingle1Trigger())) + return true; //true is correct? + + //Ask for 1 cluster p > 1.2 GeV with time [40,70] + //TODO Use Cutflow + + double minTime = 40; + double maxTime = 70; + + if (!isData_) { + minTime = 30; + maxTime = 50; + } + + if (ecal_->size() <= 2) + return true; + + bool foundFeeCluster = false; + + for (unsigned int iclu = 0; iclu < ecal_->size(); iclu++) { + if (ecal_->at(iclu)->getEnergy() > 1.5) + foundFeeCluster = true; + break; + } + + if (!foundFeeCluster) + return true; + + bool clusterInTime = true; + + for (unsigned int iclu = 0; iclu < ecal_->size(); iclu++) { + if (ecal_->at(iclu)->getTime() < 40 || ecal_->at(iclu)->getTime() > 70) + clusterInTime = false; + } + + if (!clusterInTime) + return true; + for (int itrack = 0; itrack < tracks_->size(); ++itrack) { if (trkSelector_) trkSelector_->getCutFlowHisto()->Fill(0.,weight); @@ -69,16 +147,41 @@ bool TrackingAnaProcessor::process(IEvent* ievent) { // Get a track Track* track = tracks_->at(itrack); int n2dhits_onTrack = !track->isKalmanTrack() ? track->getTrackerHitCount() * 2 : track->getTrackerHitCount(); + + TVector3 trk_mom; + trk_mom.SetX(track->getMomentum()[0]); + trk_mom.SetY(track->getMomentum()[1]); + trk_mom.SetZ(track->getMomentum()[2]); + //Track Selection if (trkSelector_ && !trkSelector_->passCutGt("n_hits_gt",n2dhits_onTrack,weight)) continue; + + if (trkSelector_ && !trkSelector_->passCutLt("chi2ndf_lt",track->getChi2Ndf(),weight)) + continue; - if (trkSelector_ && !trkSelector_->passCutGt("pt_gt",fabs(track->getPt()),weight)) - continue; + if (trkSelector_ && !trkSelector_->passCutGt("p_gt",trk_mom.Mag(),weight)) + continue; + if (trkSelector_ && !trkSelector_->passCutLt("p_lt",trk_mom.Mag(),weight)) + continue; + + if (trkSelector_ && !trkSelector_->passCutLt("trk_ecal_lt",track->getPositionAtEcal()[0],weight)) + continue; + + + if (trkSelector_ && !trkSelector_->passCutGt("trk_ecal_gt",track->getPositionAtEcal()[0],weight)) + continue; + + if (trkSelector_ && !trkSelector_->passCutGt("trk_time_gt",track->getTrackTime()-time_offset_,weight)) + continue; + + if (trkSelector_ && !trkSelector_->passCutLt("trk_time_lt",track->getTrackTime()-time_offset_,weight)) + continue; + Track* truth_track = nullptr; - + //Get the truth track if (doTruth_) { truth_track = (Track*) track->getTruthLink().GetObject(); @@ -104,11 +207,33 @@ bool TrackingAnaProcessor::process(IEvent* ievent) { } n_sel_tracks++; + + + + //Fill histograms for FEE smearing analysis + trkHistos_->Fill2DHisto("xypos_at_ecal_hh", + track->getPositionAtEcal()[0], + track->getPositionAtEcal()[1]); + + trkHistos_->Fill3DHisto("p_vs_TanLambda_Phi_hhh", + track->getPhi(), + track->getTanLambda(), + track->getP()); + + trkHistos_->Fill2DHisto("p_vs_nHits_hh", + track->getTrackerHitCount(), + track->getP()); + + trkHistos_->Fill3DHisto("p_vs_TanLambda_nHits_hhh", + track->getTanLambda(), + track->getTrackerHitCount(), + track->getP()); + + }//Loop on tracks trkHistos_->Fill1DHisto("n_tracks_h",n_sel_tracks); - return true; } @@ -125,6 +250,14 @@ void TrackingAnaProcessor::finalize() { delete truthHistos_; truthHistos_ = nullptr; } + + for (reg_it it = reg_histos_.begin(); it!=reg_histos_.end(); ++it) { + std::string dirName = it->first; + (it->second)->saveHistos(outF_,dirName); + outF_->cd(dirName.c_str()); + reg_selectors_[it->first]->getCutFlowHisto()->Write(); + } + //trkHistos_->Clear(); } diff --git a/processors/src/TrackingProcessor.cxx b/processors/src/TrackingProcessor.cxx index 859f49f66..146712391 100644 --- a/processors/src/TrackingProcessor.cxx +++ b/processors/src/TrackingProcessor.cxx @@ -26,13 +26,13 @@ void TrackingProcessor::configure(const ParameterSet& parameters) { truthTracksCollRoot_ = parameters.getString("truthTrackCollRoot",truthTracksCollRoot_); bfield_ = parameters.getDouble("bfield",bfield_); trackStateLocation_ = parameters.getString("trackStateLocation",trackStateLocation_); - + //Residual plotting is done in this processor for the moment. doResiduals_ = parameters.getInteger("doResiduals",doResiduals_); trackResDataLcio_ = parameters.getString("trackResDataLcio",trackResDataLcio_); resCfgFilename_ = parameters.getString("resPlots",resCfgFilename_); resoutname_ = parameters.getString("resoutname",resoutname_); - + } catch (std::runtime_error& error) { @@ -66,8 +66,9 @@ void TrackingProcessor::initialize(TTree* tree) { } bool TrackingProcessor::process(IEvent* ievent) { - - //Clean up + + + //Clean up if (tracks_.size() > 0 ) { for (std::vector::iterator it = tracks_.begin(); it != tracks_.end(); ++it) { delete *it; @@ -116,7 +117,7 @@ bool TrackingProcessor::process(IEvent* ievent) { // Heap an LCRelation navigator which will allow faster access rawTracker_hit_fits_nav = new UTIL::LCRelationNavigator(raw_svt_hit_fits); } - + EVENT::LCCollection* tracks{nullptr}; try { @@ -128,7 +129,7 @@ bool TrackingProcessor::process(IEvent* ievent) { std::cout << e.what() << std::endl; return false; } - + //Initialize map of shared hits std::map > SharedHits; diff --git a/processors/src/VertexAnaProcessor.cxx b/processors/src/VertexAnaProcessor.cxx index 5237426b1..913e41220 100644 --- a/processors/src/VertexAnaProcessor.cxx +++ b/processors/src/VertexAnaProcessor.cxx @@ -43,6 +43,9 @@ void VertexAnaProcessor::configure(const ParameterSet& parameters) { //region definitions regionSelections_ = parameters.getVString("regionDefinitions",regionSelections_); + //v0 projection fits + v0ProjectionFitsCfg_ = parameters.getString("v0ProjectionFitsCfg", v0ProjectionFitsCfg_); + //beamspot positions beamPosCfg_ = parameters.getString("beamPosCfg", beamPosCfg_); //track time bias corrections @@ -68,13 +71,20 @@ void VertexAnaProcessor::initialize(TTree* tree) { _vtx_histos->loadHistoConfig(histoCfg_); _vtx_histos->DefineHistos(); - if(!isData_){ + if(!isData_ && mc_reg_on_){ _mc_vtx_histos = std::make_shared(anaName_+"_mc_"+"vtxSelection"); _mc_vtx_histos->loadHistoConfig(mcHistoCfg_); _mc_vtx_histos->DefineHistos(); _mc_vtx_histos->Define2DHistos(); } + //Load Run Dependent V0 target projection fits from json + if(!v0ProjectionFitsCfg_.empty()){ + std::ifstream v0proj_file(v0ProjectionFitsCfg_); + v0proj_file >> v0proj_fits_; + v0proj_file.close(); + } + //Run Dependent Corrections //Beam Position if(!beamPosCfg_.empty()){ @@ -102,7 +112,8 @@ void VertexAnaProcessor::initialize(TTree* tree) { _reg_vtx_histos[regname]->DefineHistos(); - if(!isData_){ + bool mc_reg_on_ = false; + if(!isData_ && mc_reg_on_){ _reg_mc_vtx_histos[regname] = std::make_shared(anaName_+"_mc_"+regname); _reg_mc_vtx_histos[regname]->loadHistoConfig(mcHistoCfg_); _reg_mc_vtx_histos[regname]->DefineHistos(); @@ -130,6 +141,11 @@ void VertexAnaProcessor::initialize(TTree* tree) { _reg_tuples[regname]->addVariable("unc_vtx_cyx"); _reg_tuples[regname]->addVariable("unc_vtx_czy"); _reg_tuples[regname]->addVariable("unc_vtx_czx"); + _reg_tuples[regname]->addVariable("unc_vtx_proj_x"); + _reg_tuples[regname]->addVariable("unc_vtx_proj_y"); + _reg_tuples[regname]->addVariable("unc_vtx_proj_x_sig"); + _reg_tuples[regname]->addVariable("unc_vtx_proj_y_sig"); + _reg_tuples[regname]->addVariable("unc_vtx_proj_sig"); //track vars _reg_tuples[regname]->addVariable("unc_vtx_ele_track_p"); @@ -218,7 +234,7 @@ void VertexAnaProcessor::initialize(TTree* tree) { tree_->SetBranchAddress("EventHeader", &evth_ , &bevth_); if (brMap_.find(tsColl_.c_str()) != brMap_.end()) tree_->SetBranchAddress(tsColl_.c_str(), &ts_ , &bts_); tree_->SetBranchAddress(vtxColl_.c_str(), &vtxs_ , &bvtxs_); - tree_->SetBranchAddress(hitColl_.c_str(), &hits_ , &bhits_); + if (brMap_.find(hitColl_.c_str()) != brMap_.end()) tree_->SetBranchAddress(hitColl_.c_str(), &hits_ , &bhits_); tree_->SetBranchAddress(ecalColl_.c_str(), &ecal_ , &becal_); if(!isData_ && !mcColl_.empty()) tree_->SetBranchAddress(mcColl_.c_str() , &mcParts_, &bmcParts_); //If track collection name is empty take the tracks from the particles. TODO:: change this @@ -293,7 +309,7 @@ bool VertexAnaProcessor::process(IEvent* ievent) { } } - if (!isData_) _mc_vtx_histos->FillMCParticles(mcParts_, analysis_); + if (!isData_ && mc_reg_on_) _mc_vtx_histos->FillMCParticles(mcParts_, analysis_); } //Store processed number of events std::vector selected_vtxs; @@ -913,7 +929,7 @@ bool VertexAnaProcessor::process(IEvent* ievent) { { //Fill MC plots after all selections - if (!isData_) _reg_mc_vtx_histos[region]->FillMCParticles(mcParts_, analysis_); + if (!isData_ && mc_reg_on_) _reg_mc_vtx_histos[region]->FillMCParticles(mcParts_, analysis_); //Build map of hits and the associated MC part ids for later TRefArray ele_trk_hits = ele_trk_gbl->getSvtHits(); @@ -1076,7 +1092,7 @@ bool VertexAnaProcessor::process(IEvent* ievent) { int L2hitCode = 0; // hit code '1111' means truth in L2_ele_ax, L2_ele_ster, L2_pos_ax, L2_pos_ster if(!isData_){ //Get hit codes. Only sure this works for 2016 KF as is. - utils::get2016KFMCTruthHitCodes(ele_trk_gbl, pos_trk_gbl, hits_, L1L2hitCode, L1hitCode, L2hitCode); + utils::get2016KFMCTruthHitCodes(ele_trk_gbl, pos_trk_gbl, L1L2hitCode, L1hitCode, L2hitCode); //L1L2 truth hit selection if (!_reg_vtx_selectors[region]->passCutLt("hitCode_lt",((double)L1L2hitCode)-0.5, weight)) continue; if (!_reg_vtx_selectors[region]->passCutGt("hitCode_gt",((double)L1L2hitCode)+0.5, weight)) continue; @@ -1200,6 +1216,22 @@ bool VertexAnaProcessor::process(IEvent* ievent) { double ele_trk_z0err = ele_trk_gbl->getZ0Err(); double pos_trk_z0 = pos_trk_gbl->getZ0(); double pos_trk_z0err = pos_trk_gbl->getZ0Err(); + + //Project vertex to target + double vtx_proj_x = -999.9; + double vtx_proj_y = -999.9; + double vtx_proj_x_sig = -999.9; + double vtx_proj_y_sig = -999.9; + double vtx_proj_sig = -999.9; + if(!v0ProjectionFitsCfg_.empty()) + vtx_proj_sig = utils::v0_projection_to_target_significance(v0proj_fits_, evth_->getRunNumber(), + vtx_proj_x, vtx_proj_y, vtx_proj_x_sig, vtx_proj_y_sig, vtx->getX(), vtx->getY(), + reconz, vtx->getP().X(), vtx->getP().Y(), vtx->getP().Z()); + + _reg_vtx_histos[region]->Fill2DHisto("unc_vtx_x_v_unc_vtx_y_hh", vtx->getX(), vtx->getY()); + _reg_vtx_histos[region]->Fill2DHisto("unc_vtx_proj_x_v_unc_vtx_proj_y_hh", vtx_proj_x, vtx_proj_y); + _reg_vtx_histos[region]->Fill2DHisto("unc_vtx_proj_x_y_significance_hh", vtx_proj_x_sig, vtx_proj_y_sig); + _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_vtx_proj_significance_hh", vtx_proj_sig, reconz); _reg_vtx_histos[region]->Fill2DHisto("vtx_track_reconz_v_Z0err_hh", ele_trk_z0err, reconz); _reg_vtx_histos[region]->Fill2DHisto("vtx_track_reconz_v_Z0err_hh", pos_trk_z0err, reconz); _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_z0_hh", ele_trk_z0, reconz); @@ -1288,6 +1320,11 @@ bool VertexAnaProcessor::process(IEvent* ievent) { _reg_tuples[region]->setVariableValue("unc_vtx_cyx", cyx); _reg_tuples[region]->setVariableValue("unc_vtx_czy", czy); _reg_tuples[region]->setVariableValue("unc_vtx_czx", czx); + _reg_tuples[region]->setVariableValue("unc_vtx_proj_x", vtx_proj_x); + _reg_tuples[region]->setVariableValue("unc_vtx_proj_y", vtx_proj_y); + _reg_tuples[region]->setVariableValue("unc_vtx_proj_x_sig", vtx_proj_x_sig); + _reg_tuples[region]->setVariableValue("unc_vtx_proj_y_sig", vtx_proj_y_sig); + _reg_tuples[region]->setVariableValue("unc_vtx_proj_sig", vtx_proj_sig); //track vars _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_p", ele_trk_gbl->getP()); @@ -1363,7 +1400,7 @@ void VertexAnaProcessor::finalize() { vtxSelector->getCutFlowHisto()->Write(); outF_->cd(); - if(!isData_) + if(!isData_ && mc_reg_on_) _mc_vtx_histos->saveHistos(outF_, _mc_vtx_histos->getName()); //delete histos; //histos = nullptr; @@ -1380,7 +1417,7 @@ void VertexAnaProcessor::finalize() { } - if(!isData_){ + if(!isData_ && mc_reg_on_){ for (reg_mc_it it = _reg_mc_vtx_histos.begin(); it!=_reg_mc_vtx_histos.end(); ++it) { std::string dirName = anaName_+"_mc_"+it->first; (it->second)->saveHistos(outF_,dirName); diff --git a/processors/src/VertexProcessor.cxx b/processors/src/VertexProcessor.cxx index 248991015..329046124 100644 --- a/processors/src/VertexProcessor.cxx +++ b/processors/src/VertexProcessor.cxx @@ -25,7 +25,11 @@ void VertexProcessor::configure(const ParameterSet& parameters) { kinkRelCollLcio_ = parameters.getString("kinkRelCollLcio", kinkRelCollLcio_); trkRelCollLcio_ = parameters.getString("trkRelCollLcio", trkRelCollLcio_); trackStateLocation_= parameters.getString("trackStateLocation", trackStateLocation_); - + hitFitsCollLcio_ = parameters.getString("hitFitsCollLcio", hitFitsCollLcio_); + trkhitCollRoot_ = parameters.getString("trkhitCollRoot",trkhitCollRoot_); + rawhitCollRoot_ = parameters.getString("rawhitCollRoot",rawhitCollRoot_); + mcPartRelLcio_ = parameters.getString("mcPartRelLcio", mcPartRelLcio_); + bfield_ = parameters.getDouble("bfield",bfield_); } catch (std::runtime_error& error) { @@ -37,17 +41,38 @@ void VertexProcessor::initialize(TTree* tree) { // Add branches to tree tree->Branch(vtxCollRoot_.c_str(), &vtxs_); tree->Branch(partCollRoot_.c_str(), &parts_); + if (!trkhitCollRoot_.empty()) + tree->Branch(trkhitCollRoot_.c_str(), &hits_); + + if (!rawhitCollRoot_.empty()) + tree->Branch(rawhitCollRoot_.c_str(), &rawhits_); } bool VertexProcessor::process(IEvent* ievent) { if (debug_ > 0) std::cout << "VertexProcessor: Clear output vector" << std::endl; + + if (hits_.size() > 0) { + for (std::vector::iterator it = hits_.begin(); it != hits_.end(); ++it) { + delete *it; + } + hits_.clear(); + } + + if (rawhits_.size() > 0) { + for (std::vector::iterator it = rawhits_.begin(); it != rawhits_.end(); ++it) { + delete *it; + } + rawhits_.clear(); + } + for(int i = 0; i < vtxs_.size(); i++) delete vtxs_.at(i); vtxs_.clear(); for(int i = 0; i < parts_.size(); i++) delete parts_.at(i); parts_.clear(); Event* event = static_cast (ievent); + UTIL::LCRelationNavigator* mcPartRel_nav; // Get the collection of vertices from the LCIO event. If no such collection // exist, a DataNotAvailableException is thrown @@ -66,7 +91,7 @@ bool VertexProcessor::process(IEvent* ievent) { // Get the collection of LCRelations between GBL tracks and kink data and track data variables. EVENT::LCCollection* gbl_kink_data{nullptr}; EVENT::LCCollection* track_data{nullptr}; - + try { if (!kinkRelCollLcio_.empty()) @@ -81,10 +106,35 @@ bool VertexProcessor::process(IEvent* ievent) { std::cout<<"Failed retrieving " << kinkRelCollLcio_ <getLCEvent()->getCollectionNames(); + auto it = std::find (evColls->begin(), evColls->end(), hitFitsCollLcio_.c_str()); + bool hasFits = true; + if(it == evColls->end()) hasFits = false; + if(hasFits) + { + raw_svt_hit_fits = event->getLCCollection(hitFitsCollLcio_.c_str()); + } + //Check to see if MC Particles are in the file + it = std::find (evColls->begin(), evColls->end(), mcPartRelLcio_.c_str()); + bool hasMCParts = true; + EVENT::LCCollection* mcPartRel; + if(it == evColls->end()) hasMCParts = false; + if(hasMCParts) + { + mcPartRel = event->getLCCollection(mcPartRelLcio_.c_str()); + // Heap an LCRelation navigator which will allow faster access + mcPartRel_nav = new UTIL::LCRelationNavigator(mcPartRel); + + } + + if (debug_ > 0) std::cout << "VertexProcessor: Converting Verteces" << std::endl; for (int ivtx = 0 ; ivtx < lc_vtxs->getNumberOfElements(); ++ivtx) { @@ -99,9 +149,60 @@ bool VertexProcessor::process(IEvent* ievent) { std::vector lc_parts = lc_vtx->getAssociatedParticle()->getParticles(); for(auto lc_part : lc_parts) { - if (debug_ > 0) std::cout << "VertexProcessor: Build particle" << std::endl; - Particle * part = utils::buildParticle(lc_part,trackStateLocation_, gbl_kink_data, track_data); - if (debug_ > 0) std::cout << "VertexProcessor: Add particle" << std::endl; + if (debug_ > 0) std::cout << "VertexProcessor: Build particle" << std::endl; + Particle * part = utils::buildParticle(lc_part,trackStateLocation_, gbl_kink_data, track_data); + //============================================= + if (lc_part->getTracks().size()>0){ + EVENT::Track* lc_track = static_cast(lc_part->getTracks()[0]); + Track* track = utils::buildTrack(lc_track,trackStateLocation_,gbl_kink_data,track_data); + int nHits = 0; + if (bfield_ > 0.0) track->setMomentum(bfield_); + if (track->isKalmanTrack()) hitType = 1; //SiClusters + EVENT::TrackerHitVec lc_tracker_hits = lc_track->getTrackerHits(); + for (auto lc_tracker_hit : lc_tracker_hits) { + TrackerHit* tracker_hit = utils::buildTrackerHit(static_cast(lc_tracker_hit),rotateHits,hitType); + std::vector rawSvthitsOn3d; + utils::addRawInfoTo3dHit(tracker_hit,static_cast(lc_tracker_hit), + raw_svt_hit_fits,&rawSvthitsOn3d,hitType); + + int hitLayer = tracker_hit->getLayer(); + nHits++; + EVENT::LCObjectVec rawHits = lc_tracker_hit->getRawHits(); + for(int irawhit = 0; irawhit < rawHits.size(); ++irawhit){ + IMPL::TrackerHitImpl* rawhit = static_cast(rawHits.at(irawhit)); + if(debug_ > 0) + std::cout << "rawhit on track has lcio id: " << rawhit->id() << std::endl; + + if (hasMCParts) + { + // Get the list of fit params associated with the raw tracker hit + EVENT::LCObjectVec lc_simtrackerhits = mcPartRel_nav->getRelatedToObjects(rawhit); + + //Loop over SimTrackerHits to get MCParticles + for(int isimhit = 0; isimhit < lc_simtrackerhits.size(); isimhit++){ + IMPL::SimTrackerHitImpl* lc_simhit = static_cast(lc_simtrackerhits.at(isimhit)); + IMPL::MCParticleImpl* lc_mcp = static_cast(lc_simhit->getMCParticle()); + if(lc_mcp == nullptr) + std::cout << "mcp is null" << std::endl; + track->addMcpHit(hitLayer, lc_mcp->id()); + if(debug_ > 0) { + std::cout << "simtrackerhit lcio id: " << lc_simhit->id() << std::endl; + std::cout << "mcp lcio id: " << lc_mcp->id() << std::endl; + } + } + } + } + + track->addHitLayer(hitLayer); + hits_.push_back(tracker_hit); + rawSvthitsOn3d.clear(); + // loop on j>i tracks + } + track->setTrackerHitCount(nHits); + part->setTrack(track); + } + //============================================= + if (debug_ > 0) std::cout << "VertexProcessor: Add particle" << std::endl; parts_.push_back(part); vtx->addParticle(part); } @@ -110,6 +211,7 @@ bool VertexProcessor::process(IEvent* ievent) { vtxs_.push_back(vtx); } + if(hasMCParts) delete mcPartRel_nav; if (debug_ > 0) std::cout << "VertexProcessor: End process" << std::endl; return true; } diff --git a/processors/src/utilities.cxx b/processors/src/utilities.cxx index 08d0f1e56..4ba738858 100644 --- a/processors/src/utilities.cxx +++ b/processors/src/utilities.cxx @@ -609,6 +609,8 @@ double utils::getKalmanTrackL1Isolations(Track* track, std::vector* double closest_dt = 999999.9; double isohit_y = 999999.9; + if ( siClusters == nullptr ) + return isohit_dy; for(int j = 0; j < siClusters->size(); j++){ TrackerHit* althit = siClusters->at(j); int althit_id = althit->getID(); @@ -670,34 +672,21 @@ double utils::getKalmanTrackL1Isolations(Track* track, std::vector* return L1_stereo_iso; } -void utils::get2016KFMCTruthHitCodes(Track* ele_trk, Track* pos_trk, std::vector* hits, int& L1L2hitCode, int& L1hitCode, int& L2hitCode){ - //Build map of hits and the associated MC part ids for later - TRefArray ele_trk_hits = ele_trk->getSvtHits(); - TRefArray pos_trk_hits = pos_trk->getSvtHits(); - std::map > trueHitIDs; - for(int i = 0; i < hits->size(); i++) - { - TrackerHit* hit = hits->at(i); - trueHitIDs[hit->getID()] = hit->getMCPartIDs(); - } +void utils::get2016KFMCTruthHitCodes(Track* ele_trk, Track* pos_trk, int& L1L2hitCode, int& L1hitCode, int& L2hitCode){ //Count the number of hits per part on the ele track std::map nHits4part_ele; - for(int i = 0; i < ele_trk_hits.GetEntries(); i++) + for(int i =0; i < ele_trk->getMcpHits().size(); i++) { - TrackerHit* eleHit = (TrackerHit*)ele_trk_hits.At(i); - for(int idI = 0; idI < trueHitIDs[eleHit->getID()].size(); idI++ ) + int partID = ele_trk->getMcpHits().at(i).second; + if ( nHits4part_ele.find(partID) == nHits4part_ele.end() ) { - int partID = trueHitIDs[eleHit->getID()].at(idI); - if ( nHits4part_ele.find(partID) == nHits4part_ele.end() ) - { - // not found - nHits4part_ele[partID] = 1; - } - else - { - // found - nHits4part_ele[partID]++; - } + // not found + nHits4part_ele[partID] = 1; + } + else + { + // found + nHits4part_ele[partID]++; } } @@ -715,22 +704,18 @@ void utils::get2016KFMCTruthHitCodes(Track* ele_trk, Track* pos_trk, std::vector //Count the number of hits per part on the pos track std::map nHits4part_pos; - for(int i = 0; i < pos_trk_hits.GetEntries(); i++) + for(int i =0; i < pos_trk->getMcpHits().size(); i++) { - TrackerHit* posHit = (TrackerHit*)pos_trk_hits.At(i); - for(int idI = 0; idI < trueHitIDs[posHit->getID()].size(); idI++ ) + int partID = pos_trk->getMcpHits().at(i).second; + if ( nHits4part_pos.find(partID) == nHits4part_pos.end() ) { - int partID = trueHitIDs[posHit->getID()].at(idI); - if ( nHits4part_pos.find(partID) == nHits4part_pos.end() ) - { - // not found - nHits4part_pos[partID] = 1; - } - else - { - // found - nHits4part_pos[partID]++; - } + // not found + nHits4part_pos[partID] = 1; + } + else + { + // found + nHits4part_pos[partID]++; } } @@ -755,13 +740,17 @@ void utils::get2016KFMCTruthHitCodes(Track* ele_trk, Track* pos_trk, std::vector bool pos_trueStereoL1 = false; bool pos_trueAxialL2 = false; bool pos_trueStereoL2 = false; - if (ele_trk->isKalmanTrack()){ - for(int i = 0; i < ele_trk_hits.GetEntries(); i++) - { - TrackerHit* hit = (TrackerHit*)ele_trk_hits.At(i); - int trackhit_layer = hit->getLayer(); - int trackhit_volume = hit->getVolume(); + if(ele_trk->isKalmanTrack()){ + for(int i = 0; i < ele_trk->getMcpHits().size(); i++){ + int mcpid = ele_trk->getMcpHits().at(i).second; + int layer = ele_trk->getMcpHits().at(i).first; + int volume = -1; + if(ele_trk->isTopTrack() == 1) + volume = 0; + else + volume = 1; + bool isAxial = false; bool isStereo = false; bool isL1 = false; @@ -769,34 +758,28 @@ void utils::get2016KFMCTruthHitCodes(Track* ele_trk, Track* pos_trk, std::vector bool isGood = false; //L1 and L2 only - if(trackhit_layer < 2) + if(layer < 2) isL1 = true; - else if(trackhit_layer > 1 && trackhit_layer < 4) + else if(layer > 1 && layer < 4) isL2 = true; else continue; - if(trackhit_volume == 1){ - if(trackhit_layer%2 == 1) + if(volume == 1){ + if(layer%2 == 1) isAxial = true; else isStereo = true; } - if(trackhit_volume == 0){ - if(trackhit_layer%2 == 0) + if(volume == 0){ + if(layer%2 == 0) isAxial = true; else isStereo = true; } - //Check if truth hit is from truth track MCP - for(int idI = 0; idI < trueHitIDs[hit->getID()].size(); idI++ ) - { - int partID = trueHitIDs[hit->getID()].at(idI); - //hit is from best track MCP - if(partID == maxID_ele){ - isGood = true; - } - } + + if(mcpid == maxID_ele) + isGood = true; if(isGood){ if(isAxial){ @@ -815,14 +798,16 @@ void utils::get2016KFMCTruthHitCodes(Track* ele_trk, Track* pos_trk, std::vector } } - //Positron - if (pos_trk->isKalmanTrack()){ + if(pos_trk->isKalmanTrack()){ + for(int i = 0; i < pos_trk->getMcpHits().size(); i++){ + int mcpid = pos_trk->getMcpHits().at(i).second; + int layer = pos_trk->getMcpHits().at(i).first; + int volume = -1; + if(pos_trk->isTopTrack() == 1) + volume = 0; + else + volume = 1; - for(int i = 0; i < pos_trk_hits.GetEntries(); i++) - { - TrackerHit* hit = (TrackerHit*)pos_trk_hits.At(i); - int trackhit_layer = hit->getLayer(); - int trackhit_volume = hit->getVolume(); bool isAxial = false; bool isStereo = false; bool isL1 = false; @@ -830,34 +815,28 @@ void utils::get2016KFMCTruthHitCodes(Track* ele_trk, Track* pos_trk, std::vector bool isGood = false; //L1 and L2 only - if(trackhit_layer < 2) + if(layer < 2) isL1 = true; - else if(trackhit_layer > 1 && trackhit_layer < 4) + else if(layer > 1 && layer < 4) isL2 = true; else continue; - if(trackhit_volume == 1){ - if(trackhit_layer%2 == 1) + if(volume == 1){ + if(layer%2 == 1) isAxial = true; else isStereo = true; } - if(trackhit_volume == 0){ - if(trackhit_layer%2 == 0) + if(volume == 0){ + if(layer%2 == 0) isAxial = true; else isStereo = true; } - //Check if truth hit is from truth track MCP - for(int idI = 0; idI < trueHitIDs[hit->getID()].size(); idI++ ) - { - int partID = trueHitIDs[hit->getID()].at(idI); - //hit is from best track MCP - if(partID == maxID_pos){ - isGood = true; - } - } + + if(mcpid == maxID_pos) + isGood = true; if(isGood){ if(isAxial){ @@ -875,6 +854,8 @@ void utils::get2016KFMCTruthHitCodes(Track* ele_trk, Track* pos_trk, std::vector } } } + + //Require both Axial and Stereo truth hits to be 'Good' hit if(ele_trueAxialL1 && ele_trueStereoL1) L1L2hitCode = L1L2hitCode | (0x1 << 3); if(pos_trueAxialL1 && pos_trueStereoL1) L1L2hitCode = L1L2hitCode | (0x1 << 2); @@ -894,3 +875,47 @@ void utils::get2016KFMCTruthHitCodes(Track* ele_trk, Track* pos_trk, std::vector if(pos_trueAxialL2) L2hitCode = L2hitCode | (0x1 << 1); if(pos_trueStereoL2) L2hitCode = L2hitCode | (0x1 << 0); } + +double utils::v0_projection_to_target_significance(json v0proj_fits, int run, double &vtx_proj_x, double &vtx_proj_y, + double &vtx_proj_x_signif, double &vtx_proj_y_signif, double vtx_x, double vtx_y, double vtx_z, + double vtx_px, double vtx_py, double vtx_pz){ + //V0 Projection fit parameters are calculated externally by projecting vertices to the target z position, + //and then fitting the 2D distribution vtx_x vs vtx_y with a rotated 2D Gaussian. + //The fit parameters are defined along the rotated coordinate system. + //Therefore, the vertex position must be rotated into this coordinate system before calculating significance. + //The rotation angle corresponding to the fit is provided in the json file containing the rotated fit values. + + //Read v0 projection fits from json file + int closest_run; + for(auto entry : v0proj_fits.items()){ + int check_run = std::stoi(entry.key()); + if(check_run > run) + break; + else{ + closest_run = check_run; + } + } + double target_pos = v0proj_fits[std::to_string(closest_run)]["target_position"]; + double rot_mean_x = v0proj_fits[std::to_string(closest_run)]["rotated_mean_x"]; + double rot_mean_y = v0proj_fits[std::to_string(closest_run)]["rotated_mean_y"]; + double rot_sigma_x = v0proj_fits[std::to_string(closest_run)]["rotated_sigma_x"]; + double rot_sigma_y = v0proj_fits[std::to_string(closest_run)]["rotated_sigma_y"]; + double rotation_angle = (double)v0proj_fits[std::to_string(closest_run)]["rotation_angle_mrad"]/1000.0; + + //project vertex to target position + vtx_proj_x = vtx_x - ((vtx_z - target_pos)*(vtx_px/vtx_pz)); + vtx_proj_y = vtx_y - ((vtx_z - target_pos)*(vtx_py/vtx_pz)); + + //Rotate projected vertex by angle corresponding to run number + double rot_vtx_proj_x = vtx_proj_x*std::cos(rotation_angle) - vtx_proj_y*std::sin(rotation_angle); + double rot_vtx_proj_y = vtx_proj_x*std::sin(rotation_angle) + vtx_proj_y*std::cos(rotation_angle); + + //Calculate significance + vtx_proj_x_signif = (rot_vtx_proj_x - rot_mean_x)/rot_sigma_x; + vtx_proj_y_signif = (rot_vtx_proj_y - rot_mean_y)/rot_sigma_y; + + // + double significance = std::sqrt( vtx_proj_x_signif*vtx_proj_x_signif + vtx_proj_y_signif*vtx_proj_y_signif ); + + return significance; +} diff --git a/scripts/run_jobPool.py b/scripts/run_jobPool.py index bd499dae9..9b561c831 100644 --- a/scripts/run_jobPool.py +++ b/scripts/run_jobPool.py @@ -107,7 +107,7 @@ def launchTestsArgs(options, infilename, fileN, jobN): parser.add_option("-c", "--config", type="string", dest="configFile", help="Configuration file to use", metavar="configFile", default="dst.py") parser.add_option("-f", "--fileList", type="string", dest="fileList", - help="List of files to run on.", metavar="fileList", default="fileList.txt") + help="List of files to run on.", metavar="fileList", default="") parser.add_option("-m", "--fileMod", type="string", dest="fileMod", help="Modifier for output file names.", metavar="fileMod", default="_hpstr") parser.add_option("-r", "--fileExt", dest="fileExt", @@ -129,14 +129,21 @@ def launchTestsArgs(options, infilename, fileN, jobN): cfgList = [] fnList = range(1, 10001) - - listfiles = glob.glob(options.inDir+"/*"+options.fileExt) - if (len(listfiles) == 0): - print("Try */*"+options.fileExt) - listfiles = glob.glob(options.inDir+"/*/*"+options.fileExt) - - print(options.inDir) - + listfiles=[] + + if (options.fileList==""): + + listfiles = glob.glob(options.inDir+"/*"+options.fileExt) + if (len(listfiles) == 0): + print("Try */*"+options.fileExt) + listfiles = glob.glob(options.inDir+"/*/*"+options.fileExt) + + print(options.inDir) + else: + with open(options.fileList,'r') as inputFileList: + for line in inputFileList: + listfiles.append(line.strip()) + fnList = range(1, len(listfiles)+1) #create folder if doesn't exists