import ROOT, glob, math, sys, os, argparse
#from mpl_toolkits.mplot3d import Axes3D
#import matplotlib.pyplot as plt

def drawCdcSectorBorders():
    ll=[]
    lc=[]
    for iB in range(0,16):
        iX=-20.25*math.cos(iB*22.5/180*math.pi)
        iY=20.25*math.sin(iB*22.5/180*math.pi)

        oX=-80*math.cos(8./180*math.pi+iB*22.5/180*math.pi)
        oY=80*math.sin(8./180*math.pi+iB*22.5/180*math.pi)

        iXc=-20.25*math.cos(iB*22.5/180*math.pi+11.25/180*math.pi)
        iYc=20.25*math.sin(iB*22.5/180*math.pi+11.25/180*math.pi)
        oXc=-80*math.cos(8./180*math.pi+iB*22.5/180*math.pi+11.25/180*math.pi)
        oYc=80*math.sin(8./180*math.pi+iB*22.5/180*math.pi+11.25/180*math.pi)
###        print iX,iY,oX,oY
        ll.append(ROOT.TLine(iX,iY,oX,oY))
        ll[iB].SetLineColor(ROOT.kRed)
        ll[iB].Draw()
        lc.append(ROOT.TLine(iXc,iYc,oXc,oYc))
        lc[iB].SetLineColor(ROOT.kBlack)
        lc[iB].Draw()
    return ll,lc

parser=argparse.ArgumentParser(description="A simple 2D and 3D event display using pyplot and Root-graphs")
parser.add_argument("-p","--path",help="path to the reco files of interest", default="/nfs/mds/data/tpc/fopi/2011/sdorheim/cdcHitPersistence/")
parser.add_argument("-f","--file-name",help="name of reco file")
parser.add_argument("-pFopi","--path-fopi",help="path to the reco files of interest", default="/nfs/mds/data/tpc/fopi/2011/fopi_reco_sverre/rev170/")
parser.add_argument("-fFopi","--file-name-fopi",help="name of reco file", default="")
parser.add_argument("-cdcTB","--cdc-track-branch",help="name of cdc genfit track branch", default="CdcTrackPostFit")
parser.add_argument("-cdcPB","--cdc-point-branch", help="name of cdc reco point branch", default="CdcSpacepoint")
parser.add_argument("-tpcTB","--tpc-track-branch",help="name of tpc genfit track branch", default="TpcTrackPostFit")
parser.add_argument("-tpcPB","--tpc-point-branch", help="name of tpc reco point branch", default="TpcSPHit")
parser.add_argument("-fopiTB","--fopi-track-branch",help="name of Fopi genfit track branch", default="FopiTrackPostFit")
parser.add_argument("-tupleB","--matching-tuple-branch",help="name of matching tuple branch", default="FopiTrackTuples")
parser.add_argument("-alFile","--alignmentFile",help="name of optional alignmentfile ", default="")
args=parser.parse_args()
print args
#getting list of files


ROOT.gROOT.ProcessLine(".x $VMCWORKDIR/gconfig/rootlogon.C")
alignment=False
#ROOT.gROOT.SetBatch()
if len(args.alignmentFile)!=0:
    ROOT.TpcAlignmentManager.init(args.alignmentFile)
    alMan=ROOT.TpcAlignmentManager.getInstance()
    alignment=True
#Getting absolute path of script
#path_=os.path.dirname(sys.argv[0])
#path_=os.path.abspath(path)
#load workaround to return by reference used in Genfit
ROOT.gROOT.ProcessLine(".L $VMCWORKDIR/macro/tpc/FOPI/cand_workaround.C+")
ROOT.gROOT.ProcessLine(".L $VMCWORKDIR/macro/tpc/FOPI/trackExtrapWrapper.C+")
full_filename=args.path+"/"+args.file_name
gfile=ROOT.TFile(full_filename,"read")
fopi_full_filename=args.path_fopi+"/"+args.file_name_fopi
gFopiFile=ROOT.TFile(fopi_full_filename,"read")
if gfile.IsZombie() :
    print "zooooombie"
    exit()
dummy=ROOT.TVector3()
tree=gfile.Get("cbmsim")
if tree == None:
    print "no tree"
    exit()
drawCdcPoints=True
if gFopiFile.IsZombie():
    print "Fopi file zoombie, no drawing of cdc points"
    drawCdcPoints=False
if drawCdcPoints:
    fopiTree=gFopiFile.Get("FOPITree")
    if fopiTree==None:
            print "Fopi tree not there, no drawing of cdc points"
            drawCdcPoints=False
if drawCdcPoints:
    tcaCdcHits=ROOT.TClonesArray("CdcHit")
    fopiTree.SetBranchAddress("CdcHitBr",ROOT.AddressOf(tcaCdcHits))
    major="FopiEventBr.GetSpillNum()"
    minor="FopiEventBr.GetEvtNum()"
    fopiTree.BuildIndex(major,minor)
    fopiIndex=fopiTree.GetTreeIndex()
                              
    

#dummy reco:
fRun=ROOT.FairRunAna()
fRun.SetInputFile(full_filename)
fMagField = ROOT.FOPIField(0.616); 
fMagField.SetTargetOffset(-61);
fRun.SetField(fMagField)
ROOT.GFFieldManager.getInstance().init(ROOT.PndFieldAdaptor(fRun.GetField()));
ROOT.GFMaterialEffects.getInstance().init(ROOT.GFTGeoMaterialInterface());
taroff=-61
targetpos=taroff+41
geoFile="/tmp/FopiGeom_s339"+str(targetpos)+".root"
fRun.SetGeomFile(geoFile)
rtdb = fRun.GetRuntimeDb();
parInput1=ROOT.FairParRootFileIo()
parInput1.open(full_filename.replace("\.reco","\.param"))
rtdb.setFirstInput(parInput1);
fRun.SetOutputFile("dummy.root");
fRun.Init()
#initialize TClonesArrays used
tcaCombTracks=ROOT.TClonesArray("GFTrack")
tcaCdcTracks=ROOT.TClonesArray("GFTrack")
tcaTpcTracks=ROOT.TClonesArray("GFTrack")
tcaTpcPoints=ROOT.TClonesArray("TpcSPHit")
tcaCdcPointsGF=ROOT.TClonesArray("CdcSpacepoint")
tcaRpcPoints=ROOT.TClonesArray("RpcPixHit")
tcaRpcPoints=ROOT.TClonesArray("RpcPixHit")
tcaBarPoints=ROOT.TClonesArray("BarPixHit")
tcaTuples=ROOT.TClonesArray("MatchingTuple")
tcaTpcEventId=ROOT.TClonesArray("TpcEventIdentifier")
tree.SetBranchAddress("TpcEventIdentifier",ROOT.AddressOf(tcaTpcEventId))
tree.SetBranchAddress(args.cdc_track_branch,ROOT.AddressOf(tcaCdcTracks))
tree.CdcTrack=tcaCdcTracks
tree.SetBranchAddress(args.tpc_track_branch,ROOT.AddressOf(tcaTpcTracks))
tree.TpcTrackPostFit=tcaTpcTracks
tree.SetBranchAddress(args.fopi_track_branch,ROOT.AddressOf(tcaCombTracks))
tree.FopiTrackPostFit=tcaCombTracks
#    tree.SetBranchAddress("TpcCluster",ROOT.AddressOf(tcaTpcClusters))
#    tree.TpcCluster=tcaTpcClusters
tree.SetBranchAddress(args.tpc_point_branch,ROOT.AddressOf(tcaTpcPoints))
tree.TpcSPHit=tcaTpcPoints
#    tree.SetBranchAddress("CdcHit",ROOT.AddressOf(tcaCdcPoints))
#    tree.CdcHit=tcaCdcPoints
tree.SetBranchAddress(args.cdc_point_branch,ROOT.AddressOf(tcaCdcPointsGF))
tree.CdcPseudoSPHit=tcaCdcPointsGF
tree.SetBranchAddress(args.matching_tuple_branch,ROOT.AddressOf(tcaTuples))
tree.SetBranchAddress("RpcPixHit", ROOT.AddressOf(tcaRpcPoints))#12
tree.SetBranchAddress("BarPixHit", ROOT.AddressOf(tcaBarPoints))#13


#simple colorMap for matplotlib                          
colors=[]
colors.append('r')
colors.append('b')
colors.append('g')
colors.append('c')
colors.append('m')
colors.append('y')
colors.append('k')
colors.append('w')
colors.append('#FF6600')#orange

totentries=tree.GetEntries()
totentriesFop=fopiTree.GetEntries()
print totentries,totentriesFop

#canvas for 2D event display
c=ROOT.TCanvas("All tracks 2D","All tracks 2D",0,0,1000,1000)
c1=ROOT.TCanvas("All tracks RZ","All tracks RZ",1000,0,1000,1000)
c2=ROOT.TCanvas("Matched tracks 2D","Matched tracks 2D",1920,000,1000,1000)
c3=ROOT.TCanvas("Matched tracks RZ","Matched tracks RZ",3000,000,1000,1000)
i=0
def dPhi(a,b):
    ret=a-b
    if ret > 180:
        ret-=360
    elif ret < -180:
        ret+= 180
    return ret
ndPhiPos=0
ndPhiMin=0
while i < totentries:
    tree.GetEntry(i)
    evId=tcaTpcEventId.At(0)
    fopiEntry=fopiIndex.GetEntryNumberWithIndex(int(evId.getSpill()),int(evId.getEventInSpill()))
    print "getting fopiEntry", fopiEntry, evId.getSpill(),evId.getEventInSpill(), drawCdcPoints
    fopiTree.GetEvent(fopiEntry)
      
#       if i%500==0:
    print i
    i=i+1
    if tcaCombTracks == None or tcaCdcTracks == None or tcaTpcTracks == None :
        continue
    
    nTpc=tcaTpcTracks.GetEntriesFast()
    nCdc=tcaCdcTracks.GetEntriesFast()
    nComb=tcaCombTracks.GetEntriesFast()
    
    #condition can be removed if needed
    if True:
        print "Event: "+str(i)+"\t nComb: "+str(nComb) +"\t nCdc: "+str(nCdc)+ "\t nTpc: " + str(nTpc)
        #Plot all TPC and CDC tracks, matched and unmatched in 2D
        c.cd()
        hist=ROOT.TH2D("All_tracks_2D","All tracks ;x (cm);y (cm)",1000,-105,105,1000,-105,105)
        histRZ=ROOT.TH2D("All_tracks_RZ","All tracks RZ;z (cm);r (cm)",1000,-105,105,1000,0,125)
        #plotting the TPC-tracks
        graphs = {} #dictionary for graphs with the points of the tpc tracks
        graphsCdc = {}
        graphsCdcTracks={}
        graphsRZ = {} #dictionary for graphs with the points of the tpc tracks
        graphsCdcRZ = {}
        graphsCdcTracksRZ={}
        trackNb=0
        hist.Draw()
        innerEdgeTpc1=ROOT.TEllipse(0,0,5,5)
        innerEdgeTpc1.SetFillStyle(0)
        outerEdgeTpc1=ROOT.TEllipse(0,0,15,15)
        outerEdgeTpc1.SetFillStyle(0)
        target1=ROOT.TBox(-2.25,-2.25,2.25,2.25)
        target1.SetFillColor(ROOT.kGray)
        target1.Draw()
        outerEdgeTpc1.Draw()
        innerEdgeTpc1.Draw()
        a,b= drawCdcSectorBorders()
        c1.cd()
        histRZ.Draw()
        target4=ROOT.TBox(-.5,0,.5,2.25)
        target4.SetFillColor(ROOT.kGray)
        target4.Draw()
        tpc4=ROOT.TBox(-40,5,30,15)
        tpc4.SetFillColor(ROOT.kGray)
        tpc4.Draw()

        hist.GetXaxis().SetRangeUser(-105,105)
        hist.GetYaxis().SetRangeUser(-105,105)
        tpcDetId=0
        for track in tcaTpcTracks:
            cand=track.getCand()
            graphs[trackNb]=ROOT.TGraph()
            graphs[trackNb].SetMarkerColor(trackNb+1)
            graphsRZ[trackNb]=ROOT.TGraph()
            graphsRZ[trackNb].SetMarkerColor(trackNb+1)
            
            for hitNb in range (cand.getNHits()):
                hitId=ROOT.workAroundGetHitId(cand,hitNb)
                tpcDetId=ROOT.workAroundGetDetId(cand,hitNb)
                hit=tcaTpcPoints.At(hitId)
                coord=hit.getRawHitCoord()
                coordVec=ROOT.TVector3(coord[0],coord[1],coord[2])
                if alignment:
                    coordVec=alMan.masterToLocal("orig",coordVec)
                    coordVec=alMan.localToMaster("tpc",coordVec)
                #                    print coord[0],coord[1]
                graphs[trackNb].SetPoint(hitNb,coordVec.x(),coordVec.y())
                graphsRZ[trackNb].SetPoint(hitNb,coordVec.z(),coordVec.Perp())                
                #print str(coord[0][0])+", "+str(coord[1][0])+", "+str(coord[2][0])
            print "Drawing tpc graph nb:"+str(trackNb)+" .GetN()  "+ str(graphs[trackNb].GetN())
            c.cd()
            graphs[trackNb].Draw("*")
            c1.cd()
            graphsRZ[trackNb].Draw("*")
            trackNb+=1
        #plotting the CDC-tracks
        trackNbCdc=0
        
        for track in tcaCdcTracks:
            graphsCdc[trackNbCdc]=ROOT.TGraph()
            graphsCdc[trackNbCdc].SetMarkerColor(trackNbCdc+1)
            graphsCdcTracks[trackNbCdc]=ROOT.TGraph()
            graphsCdcTracks[trackNbCdc].SetMarkerColor(trackNbCdc+1)
            graphsCdcTracks[trackNbCdc].SetLineColor(trackNbCdc+1)
            graphsCdcRZ[trackNbCdc]=ROOT.TGraph()
            graphsCdcRZ[trackNbCdc].SetMarkerColor(trackNbCdc+1)
            graphsCdcTracksRZ[trackNbCdc]=ROOT.TGraph()
            graphsCdcTracksRZ[trackNbCdc].SetMarkerColor(trackNbCdc+1)
            graphsCdcTracksRZ[trackNbCdc].SetLineColor(trackNbCdc+1)
            hitNb=0
            cand=track.getCand()
            hitCounter=-1
            for hitNb in range (cand.getNHits()):
                hitId=ROOT.workAroundGetHitId(cand,hitNb)
                cdcDetId=ROOT.workAroundGetDetId(cand,hitNb)

                if cdcDetId==22:
                    hit=tcaCdcHits.At(hitId)
                
                    if hit !=None:
                        hitCounter+=1
                        coord=hit.GetHitPos()
                        graphsCdc[trackNbCdc].SetPoint(hitCounter,coord.X(),coord.Y())
                        graphsCdcRZ[trackNbCdc].SetPoint(hitCounter,coord.Z(),coord.Perp())
                if cdcDetId==12 or cdcDetId==13:
                    hit=None
                    if cdcDetId==12:
                        hit=tcaRpcPoints.At(hitId)
                    elif cdcDetId==13:
                        hit=tcaBarPoints.At(hitId)
                    if hit !=None:
                        hitCounter+=1
                        coord=hit.getDetPlane(ROOT.RKTrackRep()).getO()
                        graphsCdc[trackNbCdc].SetPoint(hitCounter,coord.X(),coord.Y())
                        graphsCdcRZ[trackNbCdc].SetPoint(hitCounter,coord.Z(),coord.Perp())
                        
                
                
                        
            rep=track.getCardinalRep().Clone()
            repOut=track.getCardinalRep().Clone()

            if rep.getMom(rep.getReferencePlane()).Perp()>0.01 and rep.getStatusFlag()==0 and rep.getPos(rep.getReferencePlane()).Perp()<40 :
                try:
                    lastPos=ROOT.TVector3(99,99,99)
                    for iStep in range(0,201):
                        pos=rep.stepalongChange(-0.2)#,dummy,dummy)
#                        pos.Print()
                        if pos.Perp()<lastPos.Perp():
                            lastPos=pos
                        else:
                            pos=lastPos
                        
                        graphsCdcTracks[trackNbCdc].SetPoint(200-iStep,pos.X(),pos.Y())
                        graphsCdcTracksRZ[trackNbCdc].SetPoint(200-iStep,pos.Z(),pos.Perp())
                    for iStep in range(200,700):
                        pos=repOut.stepalongChange(0.3)#,dummy,dummy)
                        graphsCdcTracks[trackNbCdc].SetPoint(iStep,pos.X(),pos.Y())
                        graphsCdcTracksRZ[trackNbCdc].SetPoint(iStep,pos.Z(),pos.Perp())                    
                except Exception:
                        print "failed"
            else:
                print "not extrapolating track: statusFlag:",rep.getStatusFlag(),", pt: ", rep.getMom(rep.getReferencePlane()).Perp(), ", trackposR: ",  rep.getPos(rep.getReferencePlane()).Perp()
#            graphsCdcTracks[trackNbCdc].Print()
            del rep,repOut
            print "Drawing cdc graph nb:"+str(trackNbCdc)+" .GetN()  "+ str(graphsCdc[trackNbCdc].GetN()) + " and the cdc track with pt: "+ str(track.getCardinalRep().getMom(track.getCardinalRep().getReferencePlane()).Perp()) +", and theta: "+str(track.getCardinalRep().getMom(track.getCardinalRep().getReferencePlane()).Theta()*ROOT.TMath.RadToDeg())
            c.cd()
            graphsCdc[trackNbCdc].Draw("*")
            graphsCdc[trackNbCdc].SetMarkerStyle(24)
            graphsCdcTracks[trackNbCdc].Draw("C")
            graphsCdcTracks[trackNbCdc].SetLineWidth(2)
            graphsCdc[trackNbCdc].SetMarkerStyle(32)
            c1.cd()
            graphsCdcRZ[trackNbCdc].Draw("*")
            graphsCdcRZ[trackNbCdc].SetMarkerStyle(24)
            graphsCdcTracksRZ[trackNbCdc].Draw("C")
            graphsCdcTracksRZ[trackNbCdc].SetLineWidth(2)
            graphsCdcRZ[trackNbCdc].SetMarkerStyle(32)            
            trackNbCdc+=1
        c.Draw()
        c.Update()
        c1.Draw()
        c1.Update()
        #Plot all combined tracks, in 2D and 3D
        c2.cd()
        hist2=ROOT.TH2D("comb_tracks_2D","combined tracks 2d ;x (cm);y (cm)",1000,-90,90,1000,-90,90)
        hist2.Draw()
        hist2.GetXaxis().SetRangeUser(-15,15)
        hist2.GetYaxis().SetRangeUser(-15,15)
        innerEdgeTpc=ROOT.TEllipse(0,0,5,5)
        innerEdgeTpc.SetFillStyle(0)
        outerEdgeTpc=ROOT.TEllipse(0,0,15,15)
        outerEdgeTpc.SetFillStyle(0)
        outerEdgeTpc.Draw()
        innerEdgeTpc.Draw()
        target2=ROOT.TBox(-2.25,-2.25,2.25,2.25)
        target2.SetFillColor(ROOT.kGray)
        target2.Draw()
        hist3=ROOT.TH2D("comb_tracks_rz","combined tracks rz ;z (cm);r (cm)",1000,-90,90,1000,0,90)
        hist3.GetYaxis().SetRangeUser(0,60)
        c3.cd()
        hist3.Draw()
        target3=ROOT.TBox(-.5,0,.5,2.25)
        target3.SetFillColor(ROOT.kGray)
        target3.Draw()
        tpc3=ROOT.TBox(-40,5,30,15)
        tpc3.SetFillColor(ROOT.kGray)
        tpc3.Draw()

        graphsTot = {}
        graphsExtrap = {}
        graphsTot3d = {}
        graphsTotRZ = {}
        graphsExtrapRZ = {}
        
        trackNbTot=0
        #initializing 3D-matplotlib-figure
        
        #fig = plt.figure()
        #        ax = Axes3D(fig)
        cdcTpcMatches = {}
        for fopiTuple in tcaTuples:
            #arrays of points used by matplotlib
            if fopiTuple.isMatched():
                if fopiTuple.hasBranch("CdcTrack") and fopiTuple.hasBranch(args.tpc_track_branch):
                    if fopiTuple.getIdFromBranch("CdcTrack") in cdcTpcMatches.keys():
                        cdcTpcMatches[fopiTuple.getIdFromBranch("CdcTrack")].append(fopiTuple.getIdFromBranch(args.tpc_track_branch))
                    else:
                        cdcTpcMatches[fopiTuple.getIdFromBranch("CdcTrack")]=[fopiTuple.getIdFromBranch(args.tpc_track_branch)]
        iTuple=-1
        for fopiTuple in tcaTuples:
            iTuple+=1
            if not (fopiTuple.isMatched() and fopiTuple.hasBranch("CdcTrack") and fopiTuple.hasBranch(args.tpc_track_branch)):
                print "no match",fopiTuple.isMatched() , fopiTuple.hasBranch("CdcTrack") ,fopiTuple.hasBranch(args.tpc_track_branch)
                for string in  fopiTuple.getListOfBranchNames():
                    print string
                continue
            
            cdcId=fopiTuple.getIdFromBranch("CdcTrack")
            tpcId=fopiTuple.getIdFromBranch(args.tpc_track_branch)
            print tpcId,cdcId
            tpcTrack=tcaTpcTracks.At(tpcId)
            cdcTrack=tcaCdcTracks.At(cdcId)
            rep1=cdcTrack.getCardinalRep().Clone()
            repOut1=cdcTrack.getCardinalRep().Clone()
            combTrack=tcaCombTracks.At(iTuple)
            comRep=combTrack.getCardinalRep()
            # if not ROOT.hasDetId(cdcTrack.getCand(),12):
            #     print "No rpc hit ignoring"
            #     continue
            x=[]
            y=[]
            z=[]
            cand=tpcTrack.getCand()
            cdcCand=cdcTrack.getCand()
            graphsTot[trackNbTot]=ROOT.TGraph()
            graphsTot[trackNbTot].SetMarkerColor(trackNbTot+1)
            graphsExtrap[trackNbTot]=ROOT.TGraph()
            graphsExtrap[trackNbTot].SetLineColor(trackNbTot+1)
            graphsTotRZ[trackNbTot]=ROOT.TGraph()
            graphsTotRZ[trackNbTot].SetMarkerColor(trackNbTot+1)
            graphsExtrapRZ[trackNbTot]=ROOT.TGraph()
            graphsExtrapRZ[trackNbTot].SetLineColor(trackNbTot+1)
            nPoints=0
            if rep1.getMom(rep1.getReferencePlane()).Perp()<0.15 or rep1.getStatusFlag()!=0 or rep1.getPos(rep1.getReferencePlane()).Perp()>40 or cdcCand.getNHits()<45:
                continue
            for hitNb in reversed(range (cand.getNHits())):
                detId=ROOT.workAroundGetDetId(cand,hitNb)
                hitId=ROOT.workAroundGetHitId(cand,hitNb)
                hit=tcaTpcPoints.At(hitId)
                coord=hit.getRawHitCoord()
                coordVec=ROOT.TVector3(coord[0],coord[1],coordVec[2]);
                if alignment:
                    coordVec=alMan.masterToLocal("orig",coordVec)
                    if coordVec.Perp()<7 or coordVec>13:
                        continue
                    locPhi=coordVec.Phi()*ROOT.TMath.RadToDeg()
                    if locPhi>0 and locPhi<20:
                        continue
                    if locPhi>175 or locPhi<-165:
                        continue
                    coordVec=alMan.localToMaster("tpc",coordVec)

                refPoint=ROOT.extrapolateToCylinder(cdcTrack,coordVec)
                if refPoint.Perp() < 4:
                    continue
                dPhi_=dPhi(coordVec.Phi()*ROOT.TMath.RadToDeg(),refPoint.Phi()*ROOT.TMath.RadToDeg())
                
                if dPhi_<0:
                    ndPhiMin+=1
                else:
                    ndPhiPos+=1
                graphsTot[trackNbTot].SetPoint(nPoints,coordVec.x(),coordVec.y())
                graphsTotRZ[trackNbTot].SetPoint(nPoints,coordVec.z(),coordVec.Perp())
                x.append(coord[0])
                y.append(coord[1])
                z.append(coord[2])
                nPoints+=1
                    
            if drawCdcPoints:
                hitCounter=-1
                for hitNb in range(cdcCand.getNHits()):
                    detId=ROOT.workAroundGetDetId(cdcCand,hitNb)
                    if detId!=22:
                        print "detId",detId
                        continue

                    hitId=ROOT.workAroundGetHitId(cdcCand,hitNb)
                    hit=tcaCdcHits.At(hitId)
                    if hit!= None:
                        hitCounter+=1
                        coord=hit.GetHitPos()
                        graphsTot[trackNbTot].SetPoint(nPoints+hitCounter,coord.X(),coord.Y())
                        graphsTotRZ[trackNbTot].SetPoint(nPoints+hitCounter,coord.Z(),coord.Perp())
                        x.append(coord.X())
                        y.append(coord.Y())
                        z.append(coord.Z())
            #redo cdc extrap
   
            if rep1.getMom(rep1.getReferencePlane()).Perp()>0.01 and rep1.getStatusFlag()==0 and rep1.getPos(rep1.getReferencePlane()).Perp()<40:
                #and comRep.getCharge()>0:
                print "################################################################################ charge: ", str(comRep.getCharge())
                #rep1.Print()
#                rep1.getPos(rep1.getReferencePlane()).Print()
                try:
                    lastPos1=ROOT.TVector3(99,99,99)
                    lastLastPos=ROOT.TVector3(99,99,99)
                    for iStep in range(0,201):
                        pos=rep1.stepalongChange(-0.2)#,dummy,dummy)
                        if (pos.Perp()<lastPos1.Perp() or lastLastPos.Perp()<lastPos.Perp()):
                            lastLastPos=lastPos1
                            lastPos1=pos
                        else:
                            pos=lastPos1
#                        pos.Print()
                        graphsExtrap[trackNbTot].SetPoint(200-iStep,pos.X(),pos.Y())
                        graphsExtrapRZ[trackNbTot].SetPoint(200-iStep,pos.Z(),pos.Perp())
                    for iStep in range(200,600):
                        pos=repOut1.stepalongChange(0.2)#,dummy,dummy)
                        graphsExtrap[trackNbTot].SetPoint(iStep,pos.X(),pos.Y())
                        graphsExtrapRZ[trackNbTot].SetPoint(iStep,pos.Z(),pos.Perp())
                    
                except Exception, e:
                        print "failed",e

            else:
                print "not extrapolating track: statusFlag:",rep1.getStatusFlag(),", pt: ", rep1.getMom(rep1.getReferencePlane()).Perp(), ", trackposR: ",  rep1.getPos(rep1.getReferencePlane()).Perp()
            del rep1,repOut1
            #                     print "detId "+str(detId)
            #only plot tracks for wich one has color !!!
            #            if trackNbTot < len(colors):
            #                ax.scatter(x,y,z,c=colors[trackNbTot])
            c2.cd()
            graphsTot[trackNbTot].Draw("*")

            graphsExtrap[trackNbTot].Draw("*")
            graphsExtrap[trackNbTot].SetMarkerStyle(6)
            graphsExtrap[trackNbTot].SetMarkerColor(trackNbTot+1)
            c3.cd()
            graphsTotRZ[trackNbTot].Draw("*")
            graphsExtrapRZ[trackNbTot].Draw("*")
            graphsExtrapRZ[trackNbTot].SetMarkerStyle(6)
            graphsExtrapRZ[trackNbTot].SetMarkerColor(trackNbTot+1)
            trackNbTot+=1

        # ax.set_xlabel('X (cm)')
        # ax.set_ylabel('Y (cm)')
        # ax.set_zlabel('Z (cm)')
        # ax.set_xlim3d([-60,60])
        # ax.set_ylim3d([-60,60])

        c2.Draw()
        c2.Update()
        c3.Update()
        print "################################################################################  min: ",ndPhiMin,"pos:", ndPhiPos
        #plt.show() #show 3d plot
        if i%1==0:
            bah = raw_input("enter to continue")
            if 'q' in bah:
                exit()
            if '-' in bah:
                i=i-2
        del hist,hist2
            
input()