test/fakeAnalysis/main.py

0001 #!/usr/bin/env python3
0002
0003 import ROOT
0004 from array import array
0005
0006 from collections import OrderedDict
0007 import pickle
0008
0009 from Validation.RecoTrack.python.plotting.ntuple import *
0010
0011 import analysis as an
0012 import graphics as gr
0013 #import graphics_vertical as gr # Use this to plot helixes better, this, however, supports only vertical z-axis and no filters except fake filter.
0014
0015
0016 def main():
0017     '''
0018     The main function, below are commented templates of function calls to perform analysis.
0019     '''
0020
0021     ntuple = TrackingNtuple("/afs/cern.ch/work/j/jmaanpaa/public/trackingNtuple_810pre6_v5.root")
0022     #ntuple = TrackingNtuple("trackingNtuple_pre9_1000ev.root") # Uncomment this for 1000 event data
0023     pl = gr.EventPlotter()
0024
0025     ###### USEFUL COMMANDS, UNCOMMENT TO USE ######
0026
0027     PlotFakes(ntuple) # Plots the fakes in 3D one by one
0028
0029     #Create_Histograms_MaxMatched_FakeTracks(ntuple)
0030     #Create_Histograms_MaxMatched_RealTracks(ntuple)
0031
0032     #Create_Histograms_Classification(ntuple) # CREATES tHE MAIN CLASSIFICATION HISTOGRAMS
0033     #Create_Histograms_IncludedDecayHitFrac(ntuple)
0034     #Create_Histograms_UnmatchedInfo(ntuple)
0035     #Create_Histograms_PixInfo(ntuple)
0036
0037     #Create_EndOfTrackingHistograms(ntuple)
0038     #Create_EndOfTrackingHistogramsReal(ntuple)
0039     #Create_EOT_DistanceHistogram(ntuple)
0040
0041     #Create_ResolutionHistograms(ntuple)
0042
0043     #print an.Calculate_AmountOfFakes(ntuple)
0044
0045
0046     ###### NOT THAT USEFUL COMMANDS, NOT TESTED AFTER UPDATES ######
0047
0048     #ntuple_iter = iter(ntuple)#__iter__()
0049     #event = ntuple_iter.next()
0050     #tracks = event.tracks()
0051
0052     #Create_PtHistograms(ntuple)
0053
0054     #an.Save_Normalisation_Coefficients(ntuple)
0055
0056     #Create_EndOfTrackingPointsPlot(ntuple)
0057     #an.EndOfTrackingDetectorInfo(Load_File("End_list_file_layermiss.dmp"),[0],True)
0058
0059     #Create_Histograms_SharedHitsFracs_RealTracks(ntuple)
0060     #Create_Histograms_SharedHitsFracs(ntuple)
0061
0062     #PlotEvent3DHits(event, "PSG")
0063     #PlotXY(event, [-10,10])
0064     #PlotZY(event,flag="P")#, [-10,10])
0065     #PlotPixelGluedTracks3D(tracks)
0066     #PlotTracks3D(FindFakes(event))
0067     #PlotTracks3D(FindUntrackedParticles(event))
0068     #EfficiencyRate(ntuple, 100)
0069     #for track in an.FindFakes(event):
0070 #   pl.DrawTrackTest(track)
0071
0072     #pl.PlotEvent3DHits(event)
0073
0074     '''
0075     track = iter(tracks).next()
0076     pl.Plot3DHits(track)
0077     pl.Draw()
0078
0079     for track in tracks:
0080     if track.hits():
0081         #pl.PlotTrack3D(track)
0082         pl.Plot3DHelixes([track])
0083         pl.Plot3DHits(track)
0084         #pl.Plot3DHelixes(an.FindFakes(event))
0085         pl.Draw()
0086     '''
0087     #pl.ParticleTest(event.trackingParticles(),True)
0088     #pl.PlotFakes(event,1,1,None)#[])
0089     #pl.PlotZY(event)
0090     #pl.Draw()
0091
0092
0093 def Save_To_File(items, name):
0094     f = open(name,'w')
0095     pickle.dump(items,f)
0096     f.close()
0097
0098 def Load_File(name):
0099     f = open(name, 'r')
0100     items = pickle.load(f)
0101     f.close()
0102     return items
0103
0104 def PlotFakes(ntuple):
0105     ''' Plots fakes with the PlotFakes command in the graphics module. '''
0106     pl = gr.EventPlotter()
0107     for event in ntuple:
0108         # Fill the filters with corresponding class indexes or detector layers to filter everything else from the plotting
0109     pl.PlotFakes(event, True, fake_filter = [], end_filter = [], last_filter = "", particle_end_filter = "")
0110
0111 def Create_PtHistograms(ntuple):
0112     ''' Creates pt histograms for fake and true tracks. '''
0113     pl = gr.EventPlotter()
0114
0115     fake_pts, true_pts = an.TrackPt(ntuple)
0116     fake_pts = [min(pt, 4.99) for pt in fake_pts]
0117     true_pts = [min(pt, 4.99) for pt in true_pts]
0118     pl.ClassHistogram(fake_pts, "Fake p_{t}", False, False)
0119     pl.ClassHistogram(true_pts, "True p_{t}", False, False)
0120
0121 def Create_ResolutionHistograms(ntuple):
0122     '''
0123     Performs resolution analysis and plots histograms from the results.
0124     Comment and uncomment lines depending on if you have saved the files.
0125     '''
0126     pl = gr.EventPlotter()
0127
0128     res_data_fakes = an.Resolution_Analysis_Fakes(ntuple, 100)
0129     #Save_To_File(res_data_fakes, "resolutions_fake_matches1000ev.dmp")
0130     #res_data_fakes = Load_File("resolutions_fake_matches1000ev.dmp")
0131
0132     #pl.ResolutionHistograms(res_data_fakes)
0133     #pl.HistogramSigmaPt(res_data_fakes)
0134
0135     res_data_trues = an.Resolution_Analysis_Trues(ntuple, 100)
0136     #Save_To_File(res_data_trues, "resolutions_trues1000ev.dmp")
0137     #res_data_trues = Load_File("resolutions_trues1000ev.dmp")
0138
0139     res_data_fakes_NLay = an.Resolution_Analysis_Fakes(ntuple, 100, real_criterion = ["NLay",6])
0140     #Save_To_File(res_data_fakes_NLay, "resolutions_fake_NLay_matches1000ev.dmp")
0141     #res_data_fakes_NLay = Load_File("resolutions_fake_NLay_matches1000ev.dmp")
0142
0143     #pl.ResolutionHistograms(res_data_fakes)
0144     #pl.HistogramSigmaPt(res_data_fakes)
0145
0146     #pl.ResolutionHistograms(res_data_trues)
0147     #pl.HistogramSigmaPt(res_data_trues)
0148
0149     #pl.ResolutionHistogramsCombine(res_data_trues, res_data_fakes)
0150     pl.ResolutionHistogramsNormalised(res_data_trues, res_data_fakes)
0151     pl.ResolutionHistogramsNormalisedCumulative(res_data_trues, res_data_fakes)
0152
0153     #pl.HistogramSigmaPtCombined(res_data_trues, res_data_fakes)
0154     #pl.HistogramSigmaPtDual(res_data_trues, res_data_fakes)
0155     pl.ResolutionsSigmaVsPt(res_data_trues, res_data_fakes)
0156
0157     pl.ResolutionHistogramsNormalised(res_data_trues, res_data_fakes_NLay)
0158     pl.ResolutionHistogramsNormalisedCumulative(res_data_trues, res_data_fakes_NLay)
0159
0160     pl.ResolutionsSigmaVsPt(res_data_trues, res_data_fakes_NLay)
0161
0162
0163
0164 def Create_EndOfTrackingHistograms(ntuple):
0165     '''
0166     Creates various histograms related to end of tracking
0167     '''
0168     pl = gr.EventPlotter()
0169
0170     ''' Comment/uncomment these lines depending on if you have saved the end of tracking file'''
0171     end_list = an.Get_EndOfTrackingPoints(ntuple, 100, [])
0172     #Save_To_File(end_list,"End_list_file_layermiss_invalid4_fix_4class.dmp")
0173
0174     #end_list = Load_File("End_list_file_layermiss_invalid4_fix_4class.dmp") #an.Get_EndOfTrackingPoints(ntuple, 10, [])
0175
0176     #end_list = an.Get_EndOfTrackingPointsReal(ntuple, 100)
0177     #Save_To_File(end_list,"End_list_file_real_invalid.dmp")
0178
0179     #end_list = Load_File("End_list_file_real_invalid.dmp")
0180     '''
0181     for end in end_list:
0182     if end.end_class == 0 and end.last_str == end.particle_end_str and end.last_str == "BPix3":
0183         print end.last
0184         print end.particle_end
0185         print end.fake_end
0186         print end.last_str, end.particle_end_str, end.fake_end_str
0187         print end.end_class
0188         print "*****"
0189     '''
0190     an.Analyse_EOT_ParticleDoubleHit(end_list)
0191     pl.EndOfTrackingHistogram(end_list, "end_class",[],[], False, False) #last false will adjust the bpix3 filter
0192     #pl.EndOfTrackingHistogram(end_list, "invalid_reason_last",[],[],False, True)
0193     pl.EndOfTrackingHistogram(end_list, "invalid_reason_fake",[0],[],False, False)
0194     pl.EndOfTrackingHistogram(end_list, "particle_pdgId",[0],[],False, False)
0195     #pl.EndOfTrackingHistogram(end_list, "invalid_reason_particle",[],[],False, False)
0196     #pl.EndOfTrackingHistogram(end_list, "particle_end_layer",[0],[],False, True)
0197
0198     pl.EndOfTrackingHistogram(end_list, "last_layer",[0],[],False, False)# 11,13,21
0199     pl.EndOfTrackingHistogram(end_list, "fake_end_layer",[0],[],False, False)
0200     pl.EndOfTrackingHistogram(end_list, "particle_end_layer",[0],[],False, False)
0201
0202     errors = an.Analyse_EOT_Error(end_list, True, "same", False)  # Same detectior id
0203     pl.ClassHistogram(errors,"Distance between EOT fake and particle hit", False, False)
0204     errors = an.Analyse_EOT_Error(end_list, True, "different", False)  # Different detectior id
0205     pl.ClassHistogram(errors,"Distance between EOT fake and particle hit", False, False)
0206     an.EndOfTrackingDetectorInfo(end_list,[0],True)
0207
0208
0209 def Create_EndOfTrackingHistogramsReal(ntuple):
0210     '''
0211     Create end of tracking histograms for true tracks.
0212     '''
0213     pl = gr.EventPlotter()
0214
0215     ''' Comment/uncomment these lines depending on if you have saved the end of tracking file'''
0216     end_list = an.Get_EndOfTrackingPointsReal(ntuple, 100)
0217     #Save_To_File(end_list,"End_list_file_real_fixed.dmp")
0218
0219     #end_list = Load_File("End_list_file_real_fixed.dmp") #an.Get_EndOfTrackingPoints(ntuple, 10, [])
0220
0221     '''
0222     for end in end_list:
0223     if end.end_class == 0 and end.last_str == end.particle_end_str and end.last_str == "BPix3":
0224         print end.last
0225         print end.particle_end
0226         print end.fake_end
0227         print end.last_str, end.particle_end_str, end.fake_end_str
0228         print end.end_class
0229         print "*****"
0230     '''
0231     pl.EndOfTrackingHistogram(end_list, "end_class",[],[])
0232     pl.EndOfTrackingHistogram(end_list, "last_layer",[0],[],False)
0233     pl.EndOfTrackingHistogram(end_list, "fake_end_layer",[0],[],False)
0234     pl.EndOfTrackingHistogram(end_list, "particle_end_layer",[0],[],False)
0235
0236 def Create_EOT_DistanceHistogram(ntuple):
0237     pl = gr.EventPlotter()
0238     end_list = an.Get_EndOfTrackingPoints(ntuple, 100, [])
0239     #Save_To_File(end_list,"End_list_file_layermiss.dmp")
0240     #end_list = Load_File("End_list_file_layermiss.dmp")
0241
0242     errors = an.Analyse_EOT_Error(end_list, True)
0243
0244     pl.ClassHistogram(errors,"Distance between EOT fake and particle hit", False, False)
0245
0246 def Create_EndOfTrackingPointsPlot(ntuple):
0247     # NOT TESTED AFTER UPDATE
0248     pl = gr.EventPlotter()
0249
0250     end_list = an.Get_EndOfTrackingPoints(ntuple, 100, [])
0251     last = [], fake_ends = [], particle_ends = [], end_classes = [], fake_classes = []
0252     for end in end_list:
0253     last.append(end.last)
0254     fake_ends.append(end.fake_end)
0255     particle_ends.append(end.particle_end)
0256     end_classes.append(end.end_class)
0257     fake_classes.append(end.fake_class)
0258
0259     '''
0260     pl.PlotEndOfTrackingPoints3D(last,[])#, fail)
0261     pl.PlotEndOfTrackingPoints3D(fake_ends,[],[],2)
0262     pl.PlotEndOfTrackingPoints3D(particle_ends,[],[],3)
0263     #pl.PlotDetectorRange()
0264     pl.Draw()
0265     '''
0266     pl.PlotEndOfTracking3D(last, fake_ends, particle_ends, end_classes, fake_classes)
0267     pl.PlotDetectorRange()
0268     pl.Draw()
0269
0270     pl.Reset()
0271     pl.PlotEndOfTracking3D(last, fake_ends, particle_ends, end_classes, fake_classes, [0],[])
0272     pl.PlotDetectorRange()
0273     pl.Draw()
0274
0275     pl.Reset()
0276     pl.PlotEndOfTracking3D(last, fake_ends, particle_ends, end_classes, fake_classes, [],[10])
0277     pl.PlotDetectorRange()
0278     pl.Draw()
0279
0280     pl.Reset()
0281     pl.PlotEndOfTracking3D(last, fake_ends, particle_ends, end_classes, fake_classes, [],[13])
0282     pl.PlotDetectorRange()
0283     pl.Draw()
0284
0285 def Create_Histograms_IncludedDecayHitFrac(ntuple):
0286     '''
0287     Creates plots related to the "multiple matches including a decay interaction" -class.
0288     '''
0289     pl = gr.EventPlotter()
0290
0291     daughter_frac, parent_frac, total_frac = an.Calculate_IndludedDecayHitFractions(ntuple, 100)
0292     pl.ClassHistogram(daughter_frac, "Daughter share of hits from fake", False, False)
0293     pl.ClassHistogram(parent_frac, "Parent share of hits from fake", False, False)
0294     pl.ClassHistogram(total_frac, "Total share of hits from fake", False, False)
0295
0296
0297 def Create_Histograms_UnmatchedInfo(ntuple):
0298     pl = gr.EventPlotter()
0299
0300     results = an.Calculate_UnmatchedSharedHits(ntuple, 100)
0301     pl.ClassHistogram(results, "Maximum number of shared hits in unmatched and unclassified classes", False, False)
0302
0303 def Create_Histograms_SharedHitsFracs(ntuple):
0304     pl = gr.EventPlotter()
0305
0306     hits_shared, hit_fractions = an.Calculate_SharedHits(ntuple, 100, [])
0307     pl.ClassHistogram(hits_shared, "Maximum number of shared hits in all classes", False, False)
0308     pl.ClassHistogram(hit_fractions, "Shared hit fractions in all classes", False, False)
0309
0310     hits_shared, hit_fractions = an.Calculate_SharedHits(ntuple, 100, [10,11,12,13])
0311     pl.ClassHistogram(hits_shared, "Maximum number of shared hits in single match classes", False, False)
0312     pl.ClassHistogram(hit_fractions, "Shared hit fractions in single match classes", False, False)
0313
0314     hits_shared, hit_fractions = an.Calculate_SharedHits(ntuple, 100, [20,21,22,23])
0315     pl.ClassHistogram(hits_shared, "Maximum number of shared hits in multiple match classes", False, False)
0316     pl.ClassHistogram(hit_fractions, "Shared hit fractions in multiple match classes", False, False)
0317
0318     hits_shared, hit_fractions = an.Calculate_SharedHits(ntuple, 100, [10])
0319     pl.ClassHistogram(hits_shared, "Maximum number of shared hits in other single match class", False, False)
0320     pl.ClassHistogram(hit_fractions, "Shared hit fractions in other single match class", False, False)
0321
0322     hits_shared, hit_fractions = an.Calculate_SharedHits(ntuple, 100, [13])
0323     pl.ClassHistogram(hits_shared, "Maximum number of shared hits in single decay match class", False, False)
0324     pl.ClassHistogram(hit_fractions, "Shared hit fractions in single decay match class", False, False)
0325
0326 def Create_Histograms_SharedHitsFracs_RealTracks(ntuple):
0327     pl = gr.EventPlotter()
0328
0329     hits_shared, hit_fractions = an.Calculate_SharedHits_RealTracks(ntuple, 100)
0330     pl.ClassHistogram(hits_shared, "Maximum number of shared hits in real reconstructed tracks", False, False)
0331     pl.ClassHistogram(hit_fractions, "Shared hit fractions in real reconstructed tracks", False, False)
0332
0333 def Create_Histograms_MaxMatched_FakeTracks(ntuple):
0334     pl = gr.EventPlotter()
0335     #results = an.Calculate_MaxMatchedHits(ntuple, 100)
0336     #pl.ClassHistogram(results, "Max Matched Hits", False, False)
0337
0338     results = an.Calculate_MaxMatchedConsecutiveHits(ntuple, 100)
0339     pl.ClassHistogram(results, "Matched consecutive hits of fake", False, False, True, 0, 20, 20, True)
0340
0341     results = an.Calculate_MaxMatchedConsecutiveHits(ntuple, 100, True)
0342     pl.ClassHistogram(results, "Consecutive hit fraction of fake", False, False, False, 0, 1.1, 11, True)
0343
0344 def Create_Histograms_MaxMatched_RealTracks(ntuple):
0345     pl = gr.EventPlotter()
0346     results = an.Calculate_MaxMatchedHits_RealTracks(ntuple, 100)
0347     pl.ClassHistogram(results, "Matched Hits of real reconstructed tracks", False, False)
0348
0349     results = an.Calculate_MaxMatchedConsecutiveHits_RealTracks(ntuple, 100)
0350     pl.ClassHistogram(results, "Matched consecutive hits of real reconstructed track", False, False)
0351
0352 def Create_Histograms_PixInfo(ntuple):
0353     pl = gr.EventPlotter()
0354
0355     nPix_data, nLay_data = an.Calculate_MatchPixelHits(ntuple, 100)
0356     pl.ClassHistogram(nPix_data, "Number of pixelhits in matches", False, False)
0357     pl.ClassHistogram(nLay_data, "Number of pixel layers in matches", False, False)
0358
0359 def Create_Histograms_Classification(ntuple):
0360     '''
0361     Creates histograms from classification of fakes.
0362     '''
0363     pl = gr.EventPlotter()
0364     '''
0365     results = an.ClassifyEventFakes(ntuple, 10, False)
0366     pl.ClassHistogram(results, "Classification of all fakes", an.classes, True)
0367
0368     results = an.Calculate_MaxMatchedConsecutiveHits(ntuple, 100)
0369     pl.ClassHistogram(results, "Max matched consecutive hits", False, False)
0370
0371     results = an.Calculate_MaxMatchedHits(ntuple, 100)
0372     pl.ClassHistogram(results, "Max Matched Hits", False, False)
0373
0374     results = an.Calculate_MaxMatchedConsecutiveHits(ntuple, 100, True)
0375     pl.ClassHistogram(results, "Hit fraction of fake", False, False) #"Max matched consecutive hits with respect to fake", False, False)
0376     '''
0377
0378     results = an.ClassifyEventFakes(ntuple, 100, False) # ORIGINAL MATCH CRITERION
0379     #results = an.ClassifyEventFakes(ntuple, 100, False, real_criterion = ["NLay", 6]) # SHARED LAYERS (NLay) MATCH CRITERION
0380
0381     #pl.ClassHistogram(results, "Classification of all fakes", an.classes, True)
0382
0383     all_classified = OrderedDict()
0384     all_classified["No matches"] = results[0]
0385     all_classified["One match"] = sum([results[i] for i in [10,11,12,13]])
0386     all_classified["Multiple matches"] = sum([results[i] for i in [20,21,22,23]])
0387     all_classified["Loose match"] = results[-1]
0388
0389     pl.ClassHistogram(all_classified, "Classification of all fakes", False, True)
0390
0391     single_classified = OrderedDict()
0392     single_classified["Decayed and reconstructed"] = results[11]
0393     single_classified["Only reconstructed"] = results[12]
0394     single_classified["Only decayed"] = results[13]
0395     single_classified["Neither"] = results[10]
0396     pl.ClassHistogram(single_classified, "Classification of fakes with single match", False, True)
0397
0398     multi_classified = OrderedDict()
0399     multi_classified["Includes a decay interaction"] = results[21]
0400     multi_classified["All matches reconstructed"] = results[22]
0401     multi_classified["Some matches reconstructed"] = results[23]
0402     multi_classified["Other"] = results[20]
0403     pl.ClassHistogram(multi_classified, "Classification of fakes with multiple matches", False, True)
0404
0405
0406
0407 if __name__ == "__main__":
0408     main()
0409