diff --git a/python/Fedra2sndsw.py b/python/Fedra2sndsw.py new file mode 100644 index 0000000000..a15201c6a6 --- /dev/null +++ b/python/Fedra2sndsw.py @@ -0,0 +1,118 @@ +#!/usr/bin/env python +'''Conversion tools from FEDRA to SNDSW''' + +from xml.sax.handler import feature_external_ges +import ROOT as r +import fedrarootlogon +import numpy as np + +#getting class instance +emureader = r.EmulsionDet() + +def importgeofile(geofile): + '''importing geometry file, can be replaced with sndsw config function if needed''' + r.gGeoManager.Import(geofile) + + +def convertseg(seg,brickID,refplate = 60): + '''convert position and angles from an EdbSegP into the global system + Seg: EdbSegP instance; + brickID: number of the brick the track was reconstructed in (ex. 31); + returns a 6-values array (x,y,z,Vx,Vy,Vz) + refplate: which is our reference plate + ''' + detID = int(brickID * 1e3 + refplate) + + + localarr = np.array([seg.X(), seg.Y(), seg.Z()]) + globalarr = np.zeros(3) + + emureader.GetPosition(detID,localarr,globalarr) + + anglocalarr = np.array([seg.TX(), seg.TY(), 1.]) + angglobalarr = np.zeros(3) + + emureader.GetAngles(detID,anglocalarr,angglobalarr) + + + return np.concatenate([globalarr,angglobalarr]) + +def converttrack(track,brickID, refplate=60, fittedsegments = False): + '''convert all segments from a volumetrack: + track: EdbTrackP instance; + brickID: number of the brick the track was reconstructed in (ex. 31); + fittedsegments: used fitted segments instead of true segments (default False) + refplate: which is our reference plate + ''' + + nseg = track.N() + globaltracklist = [] + for iseg in range(nseg): + myseg = track.GetSegment(iseg) + if fittedsegments: #replace true segment with fittedsegment + myseg = track.GetSegmentF(iseg) + globalsegarr = convertseg(myseg,brickID,refplate) + globaltracklist.append(globalsegarr) + + globaltrackarr = np.array(globaltracklist) + return globaltrackarr + +def convertvertex(vertex, brickID, refplate=60, fittedsegments=False): + '''Converting EdbVertex position into SNDSW system + vertex: EdbVertex instance; + brickID: number of the brick the vertex was reconstructed in (ex. 31) + fittedsegments: used fitted segments instead of true segments (default False) + refplate: which is our reference plate + ''' + detID = int(brickID * 1e3 + refplate) + #vertex position + localvarr = np.array([vertex.VX(),vertex.VY(),vertex.VZ()]) + globalvarr = np.zeros(3) + emureader.GetPosition(detID,localvarr,globalvarr) + #track info + ntracks = vertex.N() + globaltracklist = [] + for itrack in range(ntracks): + track = vertex.GetTrack(itrack) + globaltrackarr = converttrack(track,brickID,refplate,fittedsegments) + globaltracklist.append(globaltrackarr) + + globaltracksarr = np.array(globaltracklist) + return globalvarr, globaltracksarr + +def convertmcpoint(emupoint, MCEventID, EmulsionID = 0, weight = 70): + ''' Converting EmuPoint from MCEventID into EdbSegP''' + + detID = emupoint.GetDetectorID() + momentum = r.TMath.Sqrt(pow(emupoint.GetPx(),2) + pow(emupoint.GetPy(),2) + pow(emupoint.GetPz(),2)) + #first, convert positions + globalpos = np.array([emupoint.GetX(),emupoint.GetY(),emupoint.GetZ()]) + localpos = np.zeros(3) + + emureader.GetLocalPosition(detID, globalpos, localpos) + + xem = localpos[0] + yem = localpos[1] + + #second, convert angles + globalang = np.array([emupoint.GetPx(),emupoint.GetPy(),emupoint.GetPz()]) + localang = np.zeros(3) + + emureader.GetLocalAngles(detID, globalang, localang) + #angles in TX, TY format + tx = localang[0]/localang[2] #px/pz + ty = localang[1]/localang[2] #py/pz + #MCTruth info + pdgcode = emupoint.PdgCode() + trackID = emupoint.GetTrackID() + + #ready to write the segment + emuseg = r.EdbSegP() + emuseg.Set(EmulsionID, xem, yem, tx, ty, weight, 1) + emuseg.SetMC(MCEventID, emupoint.GetTrackID()) + emuseg.SetP(momentum) + emuseg.SetVid(pdgcode,0) + + emuseg.SetZ(localpos[2]) #0 is the center of the volume (center of plastic base in emulsion film) + + return emuseg diff --git a/python/fedra2scipy_utils.py b/python/fedra2scipy_utils.py new file mode 100644 index 0000000000..4992dc0e47 --- /dev/null +++ b/python/fedra2scipy_utils.py @@ -0,0 +1,385 @@ +import numpy as np +import pandas as pd +import fedrarootlogon +import ROOT as r + +def builddataframe(brick, path = "..", cutstring = "1", major = 0, minor = 0, newzprojection = None, charmsim = False): + """build pandas dataframe starting from couples and scanset + brick = Number of brick as in b0000* + path = input path to the folder containing theb b0000* folder + cutsring = eventual selection to couples + newzprojection = list of projection to a new z reference system + """ + nplate =0 + print("Reading ScanSet at path ",path) + + #reading scanset + sproc = r.EdbScanProc() + sproc.eProcDirClient=path + id = r.EdbID(brick,nplate,major,minor) + ss = sproc.ReadScanSet(id) + ss.Brick().SetID(brick) + + #preparing patterns + npl = ss.eIDS.GetEntries() + + cut = r.TCut(cutstring) + + #intial empty arrays + IDall = np.zeros(0,dtype=int) + PIDall = np.zeros(0,dtype=int) + + xall = np.zeros(0,dtype=np.float32) + yall = np.zeros(0,dtype=np.float32) + zall = np.zeros(0,dtype=np.float32) + TXall = np.zeros(0,dtype=np.float32) + TYall = np.zeros(0,dtype=np.float32) + + MCEvtall = np.zeros(0,dtype=int) + MCTrackall = np.zeros(0,dtype=int) + Pall = np.zeros(0,dtype=np.float32) + PdgCodeall = np.zeros(0,dtype=int) + MCMotherIDall = np.zeros(0,dtype=int) + + print ("Cut on couples ") + cut.Print() + + print("Try to open folders at path ",path+"/b0000"+str(brick)) + for i in range(npl): + idplate = ss.GetID(i) + + nplate = idplate.ePlate + plate = ss.GetPlate(idplate.ePlate) + #read pattern information + p = r.EdbPattern() + + ect = r.EdbCouplesTree() + if (nplate) <10: + ect.InitCouplesTree("couples",path+"/b{:06d}/p00{}/{}.{}.{}.{}.cp.root".format(brick,nplate,brick,nplate,major,minor),"READ") + elif (nplate) < 100: + ect.InitCouplesTree("couples",path+"/b{:06d}/p0{}/{}.{}.{}.{}.cp.root".format(brick,nplate,brick,nplate,major,minor),"READ") + else: + ect.InitCouplesTree("couples",path+"/b{:06d}/p{}/{}.{}.{}.{}.cp.root".format(brick,nplate,brick,nplate,major,minor),"READ") + + #addingcut + ect.eCut = cut + cutlist = ect.InitCutList() + + nsegcut = cutlist.GetN() + nseg = ect.eTree.GetEntries() + + IDarray_plate = np.zeros(nsegcut,dtype=int) + PIDarray_plate = np.zeros(nsegcut,dtype=int) + + xarray_plate = np.zeros(nsegcut,dtype=np.float32) + yarray_plate = np.zeros(nsegcut,dtype=np.float32) + zarray_plate = np.zeros(nsegcut,dtype=np.float32) + TXarray_plate = np.zeros(nsegcut,dtype=np.float32) + TYarray_plate = np.zeros(nsegcut,dtype=np.float32) + + MCEvtarray_plate = np.zeros(nsegcut,dtype=int) + MCTrackarray_plate = np.zeros(nsegcut,dtype=int) + Parray_plate = np.zeros(nsegcut,dtype=np.float32) + PdgCodearray_plate = np.zeros(nsegcut,dtype=int) + MCMotherIDarray_plate = np.zeros(nsegcut,dtype=int) + + print ("loop on {} segments over {} for plate {}".format(nsegcut, nseg,nplate)) + for ientry in range(nsegcut): + iseg = cutlist.GetEntry(ientry) + ect.GetEntry(iseg) + + seg=ect.eS + #//setting z and affine transformation + seg.SetZ(plate.Z()) + seg.SetPID(i) + seg.Transform(plate.GetAffineXY()) + + #transform angles too! + afftxty = plate.GetAffineTXTY() + tx = afftxty.A11()*seg.TX() + afftxty.A12()*seg.TY() + afftxty.B1() + ty = afftxty.A21()*seg.TX() + afftxty.A22()*seg.TY() + afftxty.B2() + seg.SetTX(tx) + seg.SetTY(ty) + + if(newzprojection is not None): + seg.PropagateTo(newzprojection[i]) + + IDarray_plate[ientry] = seg.ID() + PIDarray_plate[ientry] = seg.PID() + + xarray_plate[ientry] = seg.X() + yarray_plate[ientry] = seg.Y() + zarray_plate[ientry] = seg.Z() + TXarray_plate[ientry] = seg.TX() + TYarray_plate[ientry] = seg.TY() + + MCEvtarray_plate[ientry] = seg.MCEvt() + MCTrackarray_plate[ientry] = seg.MCTrack() + Parray_plate[ientry] = seg.P() + PdgCodearray_plate[ientry] = seg.Vid(0) + MCMotherIDarray_plate[ientry] = seg.Aid(0) + + #end of loop, storing them in global arrays + IDall = np.concatenate((IDall,IDarray_plate),axis=0) + PIDall = np.concatenate((PIDall,PIDarray_plate),axis=0) + + xall = np.concatenate((xall,xarray_plate),axis=0) + yall = np.concatenate((yall,yarray_plate),axis=0) + zall = np.concatenate((zall,zarray_plate),axis=0) + TXall = np.concatenate((TXall,TXarray_plate),axis=0) + TYall = np.concatenate((TYall,TYarray_plate),axis=0) + MCEvtall = np.concatenate((MCEvtall,MCEvtarray_plate),axis=0) + MCTrackall = np.concatenate((MCTrackall,MCTrackarray_plate),axis=0) + Pall = np.concatenate((Pall,Parray_plate),axis=0) + PdgCodeall = np.concatenate((PdgCodeall,PdgCodearray_plate),axis=0) + MCMotherIDall = np.concatenate((MCMotherIDall,MCMotherIDarray_plate),axis=0) + + data = {'ID':IDall,'PID':PIDall,'x':xall,'y':yall,'z':zall,'TX':TXall,'TY':TYall, + 'MCEvent':MCEvtall,'MCTrack':MCTrackall,'MCMotherID':MCMotherIDall,'P':Pall,'PdgCode':PdgCodeall} + df = pd.DataFrame(data, columns = ['ID','PID','x','y','z','TX','TY','MCEvent','MCTrack','MCMotherID','P','PdgCode'] ) + + return df + + +def addtrueMCinfo(df,simfile, simfilebkg = 0): + '''getting additional true MC info from source file, If simfilebkg is not 0, signal/bkg merge is taken into account to assign event numbers + + ''' + import pandas as pd + import numpy as np + import ROOT as r + + simtree = simfile.Get("cbmsim") + simtreebkg = r.TTree() + if (simfilebkg): + simtreebkg = simfilebkg.Get("cbmsim") #adding merged background + + #position differences from FairShip2Fedra: initialized to 0 + xoffset = 0. + yoffset = 0. + zoffset = 0. + + evID_multiplier = 1e+3 + #computing zoffset: in our couples, most downstream plate has always z=0 + simtree.GetEntry(0) + emulsionhits = simtree.EmulsionDetPoint + ihit = 0 + while (zoffset >= 0.): + hit = emulsionhits[ihit] + if (int(str(hit.GetDetectorID())[-2:])==60): #last two digits equal to 60 + zoffset = 0. - hit.GetZ() + ihit = ihit + 1 + + #virtual TM parameters (only ship-charm simulations) + spilldy = 1 + targetmoverspeed = 2.6 + + print("ZOffset between FairShip and FEDRA",zoffset) + + df = df.sort_values(["MCEvent","MCTrack","PID"],ascending=[True,True,False]) #sorting by MCEventID allows to access each event only once + df.reset_index() + + currentevent=-1 + + nsegments = len(df) + isegment = 0 + print("dataframe prepared, starting loop over {} segments".format(nsegments)) + + #preparing arrays with new columns + arr_MotherID = np.zeros(nsegments, dtype=int) + arr_MotherPDG = np.zeros(nsegments, dtype=int) + arr_ProcID = np.zeros(nsegments, dtype=int) + + arr_startX = np.zeros(nsegments,dtype=np.float32) + arr_startY = np.zeros(nsegments,dtype=np.float32) + arr_startZ = np.zeros(nsegments,dtype=np.float32) + arr_startT = np.zeros(nsegments,dtype=np.float32) + + arr_startPx = np.zeros(nsegments,dtype=np.float32) + arr_startPy = np.zeros(nsegments,dtype=np.float32) + arr_startPz = np.zeros(nsegments,dtype=np.float32) + + for (MCEvent, MCTrack) in zip(df['MCEvent'], df['MCTrack']): + + if (MCEvent != currentevent): + currentevent = MCEvent + eventtracks = r.TClonesArray() #first initialization to empty, after filled + + if (simfilebkg and MCEvent >= evID_multiplier): #need to use bkg tree + #decoding event number + inmuonpart = int(currentevent / evID_multiplier) - 1 + simtreebkg.GetEntry(inmuonpart) #conversion formula was (imuon +1) *evID_multiplier + eventtracks = simtreebkg.MCTrack + else: #using signal tree + simtree.GetEntry(currentevent) + eventtracks = simtree.MCTrack + + if(currentevent%10000 == 0): + print("Arrived at event ",currentevent) + + if(MCTrack >= 0 ): + #adding values + mytrack = eventtracks[MCTrack] + arr_MotherID[isegment] = mytrack.GetMotherId() + arr_ProcID[isegment] = mytrack.GetProcID() + + arr_startX[isegment] = (mytrack.GetStartX() + xoffset) * 1e+4 + 473000 #we need also to convert cm to micron + arr_startY[isegment] = (mytrack.GetStartY() + yoffset) * 1e+4 - 158000 + arr_startZ[isegment] = (mytrack.GetStartZ() + zoffset) * 1e+4 + arr_startT[isegment] = mytrack.GetStartT() + + arr_startPx[isegment] = mytrack.GetPx() + arr_startPy[isegment] = mytrack.GetPy() + arr_startPz[isegment] = mytrack.GetPz() + #getting mother track info + if (arr_MotherID[isegment] >= 0): + mothertrack = eventtracks[int(arr_MotherID[isegment])] + arr_MotherPDG[isegment] = mothertrack.GetPdgCode() + if(MCTrack < 0 ): + #missing values, putting defaults + arr_MotherID[isegment] = -2 + arr_ProcID[isegment] = -2 + + arr_startX[isegment] = -1.1 + arr_startY[isegment] = -1.1 + arr_startZ[isegment] = -1.1 + arr_startT[isegment] = -1.1 + + arr_startPx[isegment] = -1.1 + arr_startPy[isegment] = -1.1 + arr_startPz[isegment] = -1.1 + + arr_MotherPDG[isegment] = -2 + + isegment = isegment + 1 + + #adding the new columns to the dataframe + df["MotherID"] = arr_MotherID + df["MotherPDG"] = arr_MotherPDG + df["ProcID"] = arr_ProcID + + df["StartX"] = arr_startX + df["StartY"] = arr_startY + df["StartZ"] = arr_startZ + df["StartTime"] = arr_startT + + df["startPx"] = arr_startPx + df["startPy"] = arr_startPy + df["startPz"] = arr_startPz + + return df + +def addtrackindex(df, trackfilename): + ''' adding track index to dataframe, if tracking was performed''' + trackfile = r.TFile.Open(trackfilename) + tracktree = trackfile.Get("tracks") + + ntracks = tracktree.GetEntries() + #initial empty arrays, to be filled with segments from all tracks + IDall = np.zeros(0,dtype=int) + PIDall = np.zeros(0,dtype=int) + TrackIDall = np.zeros(0,dtype=int) + + print("start loop on {} tracks".format(tracktree.GetEntries())) + for track in tracktree: + nseg = track.nseg + segments = track.s + + #initial arrays, length given by number of segments of this tracks + IDarr = np.zeros(nseg,dtype=int) + PIDarr = np.zeros(nseg,dtype=int) + TrackIDarr = np.zeros(nseg,dtype=int) + + #start loop on segments + for iseg, seg in enumerate(segments): + IDarr[iseg] = seg.ID() + PIDarr[iseg] = seg.PID() + TrackIDarr[iseg] = seg.Track() + #concatenate, adding segments for this track to the global arrays + IDall = np.concatenate((IDall,IDarr),axis=0) + PIDall = np.concatenate((PIDall,PIDarr),axis=0) + TrackIDall = np.concatenate((TrackIDall,TrackIDarr),axis=0) + + labels = ["ID","PID","FEDRATrackID"] + dftracks = pd.DataFrame({"ID":IDall,"PID":PIDall,"FEDRATrackID":TrackIDall},columns = labels) + print("Track dataframe ready: merging it with all couples dataframe: not tracked segments will be labelled as NA") + #Now I need to merge them, however I want to keep all the segments, not only the ones which have been tracked. Luckily, there are many ways to do a merge (default is inner) + dfwithtracks = df.merge(dftracks,how = 'left', on=["PID","ID"]) + return dfwithtracks + +def retrievetrackinfo(df, trackfilename): + '''from trackindex, make a dataframe with fitted position and angles of the tracks''' + trackfile = r.TFile.Open(trackfilename) + tracktree = trackfile.Get("tracks") + + dftracked = df.query("TrackID>=0") #we access the subset with tracks + #we keep only the first segment of eacht track + dftracked = dftracked.sort_values("PID",ascending=False) + dftracked = dftracked.groupby("TrackID").first() + dftracked = dftracked.reset_index() + #how many tracks we have? What are their TrackIDs? + ntracks = len(dftracked) + trackID = dftracked["TrackID"].to_numpy(dtype=int) + #preparing arrays to store informatino + Xarr = np.zeros(ntracks,dtype=np.float32) + Yarr = np.zeros(ntracks,dtype=np.float32) + TXarr = np.zeros(ntracks,dtype=np.float32) + TYarr = np.zeros(ntracks,dtype=np.float32) + #loop over tracks + for trid in trackID: + #retrieving trackinfo + tracktree.GetEntry(trid) + track = tracktree.t + + Xarr[trid]=track.X() + Yarr[trid]=track.Y() + TXarr[trid]=track.TX() + TYarr[trid]=track.TY() + #end of loop, preparing outputdataframe and returning it + trackdf = pd.DataFrame({"TrackID":trackID,"X":Xarr,"Y":Yarr,"TX":TXarr,"TY":TYarr},columns = ["TrackID","X","Y","TX","TY"]) + return trackdf + + + +def addvertexindex(df,vertexfilename): + '''adding vertex index to dataframe. Requires Track Index''' + vertexfile = r.TFile.Open(vertexfilename) + vertextree = vertexfile.Get("vtx") + + nvertices = vertextree.GetEntries() + #initial empty arrays, to be filled with segments from all tracks + #TrackIDall = np.zeros(0,dtype=int) + VertexS = {} + VertexE = {} + + #start loop on vertices + for vtx in vertextree: + + vID = vtx.vID + ntracks = vtx.n + trackIDs = vtx.TrackID + trackedges = vtx.incoming + #start loop on tracks + for itrack in range(ntracks): + trackID = trackIDs[itrack] + #initial values, placeholder for vertex indexes + if trackID not in VertexS: + VertexS[trackID] = -1 + if trackID not in VertexE: + VertexE[trackID] = -1 + #filling tracks + if (trackedges[itrack] == 0): + VertexE[trackID] = vID + else: + VertexS[trackID] = vID + + #getting array with all tracks + labels = ["FEDRATrackID","VertexS","VertexE"] + TrackIDarr = list(VertexS.keys()) + VertexSarr = list(VertexS.values()) + VertexEarr = list(VertexE.values()) + dfvertices = pd.DataFrame({"FEDRATrackID":TrackIDarr, "VertexS":VertexSarr, "VertexE":VertexEarr},columns = labels, dtype = int) + + #dfwithvertices = df.merge(dfvertices,how='left', on = ["TrackID"]) + + return dfvertices diff --git a/shipLHC/EmulsionDet.cxx b/shipLHC/EmulsionDet.cxx index f280cfb4d9..4929e85d27 100644 --- a/shipLHC/EmulsionDet.cxx +++ b/shipLHC/EmulsionDet.cxx @@ -124,8 +124,8 @@ void EmulsionDet::DecodeBrickID(Int_t detID, Int_t &NWall, Int_t &NRow, Int_t &N Int_t NTransverse = (detID - NWall*1E4)/1E3; switch (NTransverse){ case (1): - NColumn = 1; - NRow = 2; + NColumn = 2; + NRow = 1; break; case (2): @@ -139,8 +139,8 @@ void EmulsionDet::DecodeBrickID(Int_t detID, Int_t &NWall, Int_t &NRow, Int_t &N break; case (4): - NColumn = 2; - NRow = 1; + NColumn = 1; + NRow = 2; break; } @@ -451,4 +451,69 @@ EmulsionDetPoint* EmulsionDet::AddHit(Int_t trackID,Int_t detID, return new(clref[size]) EmulsionDetPoint(trackID,detID, pos, mom, time, length, eLoss, pdgCode,Lpos,Lmom); } + +void EmulsionDet::GetPosition(Int_t id, const Double_t* localpos, Double_t* globalpos){ + + TGeoBBox *emubox = (TGeoBBox*) gGeoManager->GetVolume("Emulsion")->GetShape(); + + Double_t EmulsionDX = emubox->GetDX(); + Double_t EmulsionDY = emubox->GetDY(); + //from emulsion micron to sndsw cm and from corner to center (our scanning takes 0,0 in a corner) + Double_t centerpos[3]; + centerpos[0] = localpos[0]*1e-4 - EmulsionDX; + centerpos[1] = localpos[1]*1e-4 - EmulsionDY; + centerpos[2] = localpos[2]*1e-4; + + //get full path + TString pathtoplate = PathBrickID(id); + TGeoNavigator* nav = gGeoManager->GetCurrentNavigator(); + //going there + nav->cd(pathtoplate.Data()); + //returning position to master + nav->LocalToMaster(centerpos, globalpos); +} + +void EmulsionDet::GetAngles(Int_t id, const Double_t* localang, Double_t* globalang){ + + //get full path + TString pathtoplate = PathBrickID(id); + TGeoNavigator* nav = gGeoManager->GetCurrentNavigator(); + //going there + nav->cd(pathtoplate.Data()); + //returning position to master + nav->LocalToMasterVect(localang, globalang); +} + +void EmulsionDet::GetLocalPosition(Int_t id, const Double_t* globalpos, Double_t* localpos){ + Double_t centerpos[3]; + //get full path + TString pathtoplate = PathBrickID(id); + TGeoNavigator* nav = gGeoManager->GetCurrentNavigator(); + //going there + nav->cd(pathtoplate.Data()); + //returning position to local + nav->MasterToLocal(globalpos, centerpos); + + TGeoBBox *emubox = (TGeoBBox*) gGeoManager->GetVolume("Emulsion")->GetShape(); + + Double_t EmulsionDX = emubox->GetDX(); + Double_t EmulsionDY = emubox->GetDY(); + + localpos[0] = (centerpos[0] + EmulsionDX)*1e+4; + localpos[1] = (centerpos[1] + EmulsionDY)*1e+4; + localpos[2] = centerpos[2]*1e+4; + + //from sndsw cm to emulsion micron and from center to corner (our scanning takes 0,0 in a corner) +} + +void EmulsionDet::GetLocalAngles(Int_t id, const Double_t* globalang, Double_t* localang){ + //get full path + TString pathtoplate = PathBrickID(id); + TGeoNavigator* nav = gGeoManager->GetCurrentNavigator(); + //going there + nav->cd(pathtoplate.Data()); + //returning position to local + nav->MasterToLocalVect(globalang, localang); +} + ClassImp(EmulsionDet) diff --git a/shipLHC/EmulsionDet.h b/shipLHC/EmulsionDet.h index 706e8f8668..4a4a881039 100644 --- a/shipLHC/EmulsionDet.h +++ b/shipLHC/EmulsionDet.h @@ -73,6 +73,12 @@ class EmulsionDet : public FairDetector EmulsionDet(const EmulsionDet&); EmulsionDet& operator=(const EmulsionDet&); + + void GetPosition(Int_t id, const Double_t* localpos, Double_t* globalpos); // use localtomaster to convert positions + void GetAngles(Int_t id, const Double_t* localang, Double_t* globalang);//converting angles + + void GetLocalPosition(Int_t id, const Double_t* globalpos, Double_t* localpos); // use mastertolocal to convert positions + void GetLocalAngles(Int_t id, const Double_t* globalang, Double_t* localang);//converting angles ClassDef(EmulsionDet,5) diff --git a/shipLHC/scripts/2dEventDisplay.py b/shipLHC/scripts/2dEventDisplay.py index b33e384c70..c00cc43c94 100644 --- a/shipLHC/scripts/2dEventDisplay.py +++ b/shipLHC/scripts/2dEventDisplay.py @@ -873,3 +873,87 @@ def drawInfo(pad, k, run, event, timestamp): if timestamp_start: textInfo.DrawLatex(0, 0.2, 'Time (GMT): {}'.format(time_event)) pad.cd(k) +gemuzx = ROOT.TGraphErrors() +gemuzy = ROOT.TGraphErrors() + +def drawEmuTrack(itrack, brickID, trackfile): + import fedrarootlogon + import numpy as np + import Fedra2sndsw as EmuConv + #needed initializations for ReadTracksTree + dproc = ROOT.EdbDataProc() + ali = ROOT.EdbPVRec() + scancond = ROOT.EdbScanCond() + ali.SetScanCond(scancond) + + dproc.ReadTracksTree(ali,trackfile,"trid=={}".format(itrack)) + #in this case, the list of tracks is only one track + track = ali.eTracks[0] + #read track, applying conversion + globaltrackarr = EmuConv.converttrack(track,brickID) + for x,y,z,tx,ty,tz in globaltrackarr: + gemuzx.AddPoint(z,x) + gemuzy.AddPoint(z,y) + + gemuzx.Print() + h['simpleDisplay'].cd(1) + gemuzx.SetMarkerStyle(20) + gemuzx.SetMarkerSize(0.1) + gemuzx.Draw("sameP") + h['simpleDisplay'].cd(2) + gemuzy.SetMarkerStyle(20) + gemuzy.SetMarkerSize(0.1) + gemuzy.Draw("sameP") + gemuzy.Print() + + h['simpleDisplay'].Update() + h['simpleDisplay'].Draw() + +def drawEmuVertex(ivertex, brickID, vertexfile): + import fedrarootlogon + import numpy as np + import Fedra2sndsw as EmuConv + #needed initializations for ReadTracksTree + dproc = ROOT.EdbDataProc() + ali = ROOT.EdbPVRec() + scancond = ROOT.EdbScanCond() + ali.SetScanCond(scancond) + + #setting parameters as in vertexing.C + vertexrec = ROOT.EdbVertexRec() + vertexrec.SetPVRec(ali) + vertexrec.eDZmax=3000. + vertexrec.eProbMin=0.0001 + vertexrec.eImpMax=15. + vertexrec.eUseMom=False + vertexrec.eUseSegPar=True + vertexrec.eQualityMode=0 + + dproc.ReadVertexTree(vertexrec, vertexfile, "vID=={}".format(ivertex)) #22044 + #in this case, the list of tracks is only one vertex + vertex = ali.eVTX[0] + ntracks = vertex.N() + for itrack in range(ntracks): + track = vertex.GetTrack(itrack) + #read track, applying conversion + globaltrackarr = EmuConv.converttrack(track,brickID) + for x,y,z,tx,ty,tz in globaltrackarr: + gemuzx.AddPoint(z,x) + gemuzy.AddPoint(z,y) + + gemuzx.Print() + h['simpleDisplay'].cd(1) + gemuzx.SetMarkerStyle(20) + gemuzx.SetMarkerSize(0.1) + gemuzx.Draw("sameP") + h['simpleDisplay'].cd(2) + gemuzy.SetMarkerStyle(20) + gemuzy.SetMarkerSize(0.1) + gemuzy.Draw("sameP") + gemuzy.Print() + + h['simpleDisplay'].Update() + h['simpleDisplay'].Draw() + + + diff --git a/shipLHC/scripts/fromsndsw2FEDRA/FairShip2Fedra.rootrc b/shipLHC/scripts/fromsndsw2FEDRA/FairShip2Fedra.rootrc new file mode 100644 index 0000000000..a901921e5f --- /dev/null +++ b/shipLHC/scripts/fromsndsw2FEDRA/FairShip2Fedra.rootrc @@ -0,0 +1,10 @@ +FairShip2Fedra.nbricks: 20 +FairShip2Fedra.nplates: 60 +FairShip2Fedra.useefficiencymap: 0 +FairShip2Fedra.emuefficiency: 0.85 +FairShip2Fedra.dosmearing: 1 +FairShip2Fedra.maxtheta: 0.78 +FairShip2Fedra.minmomentum: 0.1 +FairShip2Fedra.ngrains: 70 +FairShip2Fedra.angres: 0.003 +FairShip2Fedra.useLocalBrickCoordinates: 1 diff --git a/shipLHC/scripts/fromsndsw2FEDRA/GetEntries.C b/shipLHC/scripts/fromsndsw2FEDRA/GetEntries.C new file mode 100644 index 0000000000..b422a3f200 --- /dev/null +++ b/shipLHC/scripts/fromsndsw2FEDRA/GetEntries.C @@ -0,0 +1,16 @@ +int GetEntries(int nbrick, int nplate){ + TFile *cpfile = TFile::Open(Form("p%03i/%i.%i.0.0.cp.root",nplate,nbrick,nplate)); + TTree *couples = (TTree*) cpfile->Get("couples"); + + return couples->GetEntries(); +} + + +int GetEntries(int nbrick){ + const int nplates = 60; + for (int iplate = 1; iplate <= nplates; iplate++){ + int ncouples = GetEntries(nbrick,iplate); + if (ncouples == 0) return 0; + } //end of loop, all ok + return 1; +} diff --git a/shipLHC/scripts/fromsndsw2FEDRA/SND2FEDRA_masterconv.sh b/shipLHC/scripts/fromsndsw2FEDRA/SND2FEDRA_masterconv.sh new file mode 100644 index 0000000000..90da5ae15f --- /dev/null +++ b/shipLHC/scripts/fromsndsw2FEDRA/SND2FEDRA_masterconv.sh @@ -0,0 +1,40 @@ +#!/bin/bash + +# Bash script wrapping all of "FairShip2Fedra" conversion tools +# you'll need to copy in this folder the sim file you want to +# convert in FEDRA +RED=$'\033[0;31m' #Red coloring +NC=$'\033[0m' #No coloring + +SIMFILE=$(ls sndLHC*.root) +GEOFILE=$(ls geofile*.root) +[ ! -d "$SIMFILE" ] && echo "${RED}++${NC} You are about to convert $SIMFILE ${RED}++${NC}" +echo "${RED}++${NC} Bricks and plates folders will be created in the following directory: $(pwd) ${RED}++${NC}" +cp $SNDSW_ROOT/shipLHC/scripts/fromsndsw2FEDRA/preparebricks.sh . +cp $SNDSW_ROOT/shipLHC/scripts/fromsndsw2FEDRA/doreco.sh . +cp $SNDSW_ROOT/shipLHC/scripts/fromsndsw2FEDRA/csvconversion.sh . +cp $SNDSW_ROOT/shipLHC/scripts/fromsndsw2FEDRA/addvertexinfo.sh . +cp $SNDSW_ROOT/shipLHC/scripts/fromsndsw2FEDRA/fromsndsw2FEDRA.C . +cp $SNDSW_ROOT/shipLHC/scripts/fromsndsw2FEDRA/FairShip2Fedra.rootrc . +cp $SNDSW_ROOT/shipLHC/scripts/fromsndsw2FEDRA/GetEntries.C . +for i in $(seq 1 5) + do + mkdir b0000${i}{1..4} + mkdir b0000${i}{1..4}/p0{01..60} + done +echo "++++++++++++++++" +echo "FairShip2Fedra parameters are the following:" +cat FairShip2Fedra.rootrc +ulimit -n 1500 #to create many files together +root -l -q fromsndsw2FEDRA.C\(\"$SIMFILE\",\"$GEOFILE\"\) +echo "${RED}++${NC} Proceeding to doreco.sh ${RED}++${NC}" +source doreco.sh +echo "${RED}++${NC} Done ${RED}++${NC}" +echo "${RED}++${NC} Now performing csvconversion ${RED}++${NC}" +source csvconversion.sh +echo "${RED}++${NC} Done ${RED}++${NC}" +echo "${RED}++${NC} Now adding vertex informations ${RED}++${NC}" +source addvertexinfo.sh +echo "${RED}++${NC} Merging all vertextree files ${RED}++${NC}" +hadd vertextree_allbricks.root b0000*/vertextree.root +echo "${RED}++${NC} SND2FEDRA conversion successfully ended ${RED}++${NC}" diff --git a/shipLHC/scripts/fromsndsw2FEDRA/add_brickID_number.py b/shipLHC/scripts/fromsndsw2FEDRA/add_brickID_number.py new file mode 100644 index 0000000000..5e73a0e810 --- /dev/null +++ b/shipLHC/scripts/fromsndsw2FEDRA/add_brickID_number.py @@ -0,0 +1,18 @@ +'''small script to add information to vertextree as a friend TNtuple''' + +import ROOT as r +import fedrarootlogon +import sys + +brickID = int(sys.argv[1]) + +vtxfile = r.TFile.Open("vertextree.root","UPDATE") +vtxtree = vtxfile.Get("vtx") + +outputtree = r.TNtuple("brickinfo","Information about the brick","vID:brickID") +#main loop, filling additional information +for vtx in vtxtree: + outputtree.Fill(vtx.vID, brickID) + +vtxfile.cd() +outputtree.Write("",r.TObject.kOverwrite) diff --git a/shipLHC/scripts/fromsndsw2FEDRA/addvertexinfo.sh b/shipLHC/scripts/fromsndsw2FEDRA/addvertexinfo.sh new file mode 100644 index 0000000000..e74557b8d5 --- /dev/null +++ b/shipLHC/scripts/fromsndsw2FEDRA/addvertexinfo.sh @@ -0,0 +1,11 @@ +#!/bin/bash +brickIDs=({1..5}{1..4}) +for ibrick in $(seq 0 19) + do + echo "processing brick ${brickIDs[ibrick]}" + cd b0000${brickIDs[ibrick]} + + python $SNDSW_ROOT/shipLHC/scripts/fromsndsw2FEDRA/add_brickID_number.py ${brickIDs[ibrick]} + + cd .. + done diff --git a/shipLHC/scripts/fromsndsw2FEDRA/csvconversion.py b/shipLHC/scripts/fromsndsw2FEDRA/csvconversion.py new file mode 100644 index 0000000000..75b10306ee --- /dev/null +++ b/shipLHC/scripts/fromsndsw2FEDRA/csvconversion.py @@ -0,0 +1,27 @@ +import fedra2scipy_utils +import sys +import pandas as pd +import ROOT as r + +def applyconversion(nbrick, simfile, simbkgfile): + '''convert couples ROOT files into a csv''' + + df = fedra2scipy_utils.builddataframe(nbrick) + df = fedra2scipy_utils.addtrueMCinfo(df,simfile, simbkgfile) + df = fedra2scipy_utils.addtrackindex(df,"linked_tracks.root") + dfvertices = fedra2scipy_utils.addvertexindex(df,"vertextree.root") + + return df, dfvertices + +#the two steps can now be done together, without an intermediate file +nbrick = int(sys.argv[1]) +simfile = r.TFile.Open(sys.argv[2]) #neutrino signal simulation +simbkgfile = r.TFile.Open(sys.argv[3]) #muon background simulation +#starting all conversion steps +df,dfvertices = applyconversion(nbrick, simfile, simbkgfile) + +#df = df.drop(columns = ["P","Flag"]) +df.to_csv('brick{}.csv'.format(nbrick),index=False) +dfvertices.to_csv('brick{}_vertices.csv'.format(nbrick),index=False) + + diff --git a/shipLHC/scripts/fromsndsw2FEDRA/csvconversion.sh b/shipLHC/scripts/fromsndsw2FEDRA/csvconversion.sh new file mode 100644 index 0000000000..f5f08db0c9 --- /dev/null +++ b/shipLHC/scripts/fromsndsw2FEDRA/csvconversion.sh @@ -0,0 +1,11 @@ +#!/bin/bash +#first argument signal file, second background file +brickIDs=({1..5}{1..4}) +for ibrick in $(seq 0 19) + do + cd b0000${brickIDs[ibrick]} + cp $SNDSW_ROOT/shipLHC/scripts/fromsndsw2FEDRA/csvconversion.py ./ + echo "now converting " b0000${brickIDs[ibrick]} + python csvconversion.py ${brickIDs[ibrick]} ../$1 ../$2 #first parameter signal simulation, second background simulation + cd .. + done diff --git a/shipLHC/scripts/fromsndsw2FEDRA/doreco.sh b/shipLHC/scripts/fromsndsw2FEDRA/doreco.sh new file mode 100644 index 0000000000..b0df33b299 --- /dev/null +++ b/shipLHC/scripts/fromsndsw2FEDRA/doreco.sh @@ -0,0 +1,25 @@ +#!/bin/bash +SCRIPTFOLDER=$SNDSW_ROOT/shipLHC/scripts/fromsndsw2FEDRA/ +echo "Using default parameters for tracking, as listed here:" +cp $SCRIPTFOLDER/track.rootrc ./ +cat ./track.rootrc +brickIDs=({1..5}{1..4}) +for ibrick in $(seq 0 19) + do + echo "processing brick ${brickIDs[ibrick]}" + cd b0000${brickIDs[ibrick]} + echo "makescanset -set=${brickIDs[ibrick]}.0.0.0 -from_plate=60 -to_plate=1 -dz=-1315 -suff=cp.root" > scanset.sh + source scanset.sh + cp ../track.rootrc ./ + ncouples=$(root -l -q ../GetEntries.C\(${brickIDs[ibrick]}\)) + if [[ "${ncouples: -1}" == "0" ]]; then + echo "${RED}No hits found in brick ${brickIDs[ibrick]}, skipping${NC}" + cd .. + continue + fi + emtra -set=${brickIDs[ibrick]}.0.0.0 -new -v=2 + ln -s b0000${brickIDs[ibrick]}.0.0.0.trk.root linked_tracks.root + cp $SCRIPTFOLDER/vertexing.C ./ + root -l -q vertexing.C\(${brickIDs[ibrick]}\) + cd .. + done diff --git a/shipLHC/scripts/fromsndsw2FEDRA/fromsndsw2FEDRA.C b/shipLHC/scripts/fromsndsw2FEDRA/fromsndsw2FEDRA.C new file mode 100644 index 0000000000..73c2a404eb --- /dev/null +++ b/shipLHC/scripts/fromsndsw2FEDRA/fromsndsw2FEDRA.C @@ -0,0 +1,241 @@ +//tool to load hit from emulsion to FEDRA (13 aprile 2018) //updated to work with the new TTreeReader structure (6 March 2019) +//to use it, go in a directory and create the folders b000001/p001 to b000001/p029 +//then launch it from the directory mother of b000001 + +void fromsndsw2FEDRA(TString simfilename, TString geofilename); +void smearing (Float_t &TX, Float_t &TY, const float angres); +bool efficiency(const float emuefficiency); +bool efficiency(const float tantheta, TH1D * emuefficiency); + +void set_default(TEnv &cenv){ //setting default parameters, if not presents from file + cenv.SetValue("FairShip2Fedra.nbricks",20);//5 walls in 2x2 configuration + cenv.SetValue("FairShip2Fedra.nplates",60); + //cenv.SetValue("FairShip2Fedra.nevents",10000); // number of events to be passed to FEDRA + cenv.SetValue("FairShip2Fedra.useefficiencymap",0); + cenv.SetValue("FairShip2Fedra.emuefficiency",0.85); //only if useefficiency map is set to false + cenv.SetValue("FairShip2Fedra.dosmearing",1); + cenv.SetValue("FairShip2Fedra.maxtheta",1.); //angular max of scanning + cenv.SetValue("FairShip2Fedra.minmomentum",0.1); //do not pass particles beyond this value, track ID would be -2 + cenv.SetValue("FairShip2Fedra.ngrains",70); // to set weight + cenv.SetValue("FairShip2Fedra.angres",0.003);//used for smearing, if dosmearing = true + cenv.SetValue("FairShip2Fedra.useLocalBrickCoordinates",1); //local brick coordinates + +} + +using namespace TMath; +TRandom3 *grandom = new TRandom3(); //creating every time a TRandom3 is a bad idea +TFile *file = NULL; +TH1D *heff = NULL ; //efficiency at different angles +EmulsionDet emureader; +//starting script +void fromsndsw2FEDRA(TString simfilename, TString geofilename){ + + //now need also geofile for transfomrations + gGeoManager->Import(geofilename); + + grandom->SetSeed(0); + + TEnv cenv("FairShip2Fedra"); + set_default(cenv); + cenv.ReadFile("FairShip2Fedra.rootrc" ,kEnvLocal); + + //getting options from file + //const Int_t nevents = cenv.GetValue("FairShip2Fedra.nevents",10000); + const int nplates = cenv.GetValue("FairShip2Fedra.nplates",60); + const int nbricks = cenv.GetValue("FairShip2Fedra.nbricks",20); // to set b00000%i number + + float angres = cenv.GetValue("FairShip2Fedra.angres",0.003); //Used cases: 3, 5milliradians. Constant value overwritten if useresfunction=true + float minmomentum = cenv.GetValue("FairShip2Fedra.minmomentum",0.1); + float maxtheta = cenv.GetValue("FairShip2Fedra.maxtheta",1.); + const float ngrains = cenv.GetValue("FairShip2Fedra.ngrains",70.) ; //the same number for all the couples, so they have the same weigth. + const float emuefficiency = cenv.GetValue("FairShip2Fedra.emuefficiency",0.85); // flat value + + const bool useefficiencymap = cenv.GetValue("FairShip2Fedra.useefficiencymap",0); //use the map instead of the constant value down + const bool dosmearing = cenv.GetValue("FairShip2Fedra.dosmearing",1); //gaussian smearing or not + const bool useLocalBrickCoordinates = cenv.GetValue("FairShip2Fedra.useLocalBrickCoordinates",1); //local brick coordinates + + cout<<"Starting conversion with efficiency "<Get("heff"); + } + + map brickindex; + //building map + brickindex[11] = 0; + brickindex[12] = 1; + brickindex[13] = 2; + brickindex[14] = 3; + brickindex[21] = 4; + brickindex[22] = 5; + brickindex[23] = 6; + brickindex[24] = 7; + brickindex[31] = 8; + brickindex[32] = 9; + brickindex[33] = 10; + brickindex[34] = 11; + brickindex[41] = 12; + brickindex[42] = 13; + brickindex[43] = 14; + brickindex[44] = 15; + brickindex[51] = 16; + brickindex[52] = 17; + brickindex[53] = 18; + brickindex[54] = 19; + + int ihit = 0; + // ***********************CREATING FEDRA TREES************************** + //gInterpreter->AddIncludePath("/afs/cern.ch/work/a/aiulian/public/fedra/include"); + EdbCouplesTree *ect[nbricks][nplates]; //2D array->which brick and which plate? + int brickIDs[20] = {11,12,13,14,21,22,23,24,31,32,33,34,41,42,43,44,51,52,53,54}; + + for (int ibrick = 0; ibrick < nbricks; ibrick++){ + for (int i = 1; i <= nplates; i++){ + ect[ibrick][i-1] = new EdbCouplesTree(); + ect[ibrick][i-1]->InitCouplesTree("couples",Form("b0000%02i/p0%02i/%i.%i.0.0.cp.root",brickIDs[ibrick],i,brickIDs[ibrick],i),"RECREATE"); //i learned padding with %0i + } + } + + TFile * inputfile = TFile::Open(simfilename.Data()); + if (!inputfile) return; + + cout<<"opening file "< tracks(reader,"MCTrack"); + TTreeReaderArray emulsionhits(reader,"EmulsionDetPoint"); + + int nevents = reader.GetEntries(); + cout<<"Start processing nevents: "<= 0) motherID = tracks[trackID].GetMotherId(); + else motherID = -2; //hope I do not see them + + int NWall = 0, NRow = 0, NColumn = 0, NPlate = 0; + + emureader.DecodeBrickID(detID, NWall, NRow, NColumn, NPlate); //getting info about brick + int nbrick = (detID - NPlate)*1e-3; + + int pdgcode = emupoint.PdgCode(); + + if ((TDatabasePDG::Instance()->GetParticle(pdgcode))!=NULL){ + charge = TDatabasePDG::Instance()->GetParticle(pdgcode)->Charge(); + mass = TDatabasePDG::Instance()->GetParticle(pdgcode)->Mass(); + } + else{ + charge = 0.; + mass = 0.; + } + float kinenergy = TMath::Sqrt(pow(mass,2)+pow(momentum,2)) - mass; + + //first, convert positions + float xem, yem, tx, ty; + if (useLocalBrickCoordinates){ //each brick has its own reference system + double globalpos[3] = {emupoint.GetX(),emupoint.GetY(),emupoint.GetZ()}; + double localpos[3] = {0,0,0}; + + emureader.GetLocalPosition(detID, globalpos, localpos); + + xem = localpos[0]; + yem = localpos[1]; + + //second, convert angles + double globalang[3] = {emupoint.GetPx(),emupoint.GetPy(),emupoint.GetPz()}; + double localang[3] = {0,0,0}; + + emureader.GetLocalAngles(detID, globalang, localang); + //angles in TX, TY format + tx = localang[0]/localang[2]; + ty = localang[1]/localang[2]; + } + else{ //global reference system as sndsw, but in FEDRA units + xem = emupoint.GetX(); + yem = emupoint.GetY(); + + xem = xem* 1E+4 + 473000; + yem = yem* 1E+4 - 158000; + + tx = emupoint.GetPx()/emupoint.GetPz(); + ty = emupoint.GetPy()/emupoint.GetPz(); + } + float tantheta = pow(pow(tx,2) + pow(ty,2),0.5); + + // *************EXCLUDE HITS FROM BEING SAVED******************* + if (tantheta > TMath::Tan(maxtheta)) savehit = false; //we scan from theta 0 to a maximum of 1 rad + if(charge == 0.) savehit = false; //we do not track neutral particles + if(momentum < minmomentum) savehit = false; //particles with too low kin energy + //saving the hits for a plate in the corresponding couples (only one layer saved, the other has ID + 10000) + + if (useefficiencymap){ //efficiency map with angle + if(!efficiency(tantheta, heff)) savehit = false; //inserting some holes due to emulsion inefficiency + } + //constant value of efficiency + else if(!efficiency(emuefficiency)) savehit = false; + if (dosmearing){ + smearing(tx,ty,angres); + } + if (savehit){ + int whichbrick = brickindex[nbrick]; //finding index of the array for the brick of our hit; + ect[whichbrick][NPlate-1]->eS->Set(ihit,xem,yem,tx,ty,1,Flag); + ect[whichbrick][NPlate-1]->eS->SetMC(inECCevent, trackID); //objects used to store MC true information + ect[whichbrick][NPlate-1]->eS->SetAid(motherID, 0); //forcing areaID member to store mother MC track information + ect[whichbrick][NPlate-1]->eS->SetP(momentum); + ect[whichbrick][NPlate-1]->eS->SetVid(pdgcode,0); //forcing viewID[0] member to store pdgcode information + ect[whichbrick][NPlate-1]->eS->SetW(ngrains); //need a high weight to do tracking + ect[whichbrick][NPlate-1]->Fill(); + ihit++; //hit entry, increasing as the tree is filled + } + }//end of loop on emulsion points + } //end of loop on tree + for (int ibrick = 0; ibrick < nbricks; ibrick++){ + for (int iplate = 0; iplate < nplates; iplate++){ + //ect[nbrick][iplate]->Write(); + ect[ibrick][iplate]->Close(); + } + } + cout<<"end of script, saving rootrc wih used parameters"<Gaus(0,angres); //angular resolution, adding a gaussian offset to TX and TY + float deltaTY = grandom->Gaus(0,angres); + //cout<Uniform(0,1); + if (prob < emuefficiency) return true; //efficiency larger than probability, we take the event + else return false; +} + +bool efficiency(const float tantheta, TH1D * emuefficiency){ //for now, just a constant, to be replaced with an efficiency map (probably with the angle) + + float prob = grandom->Uniform(0,1); + int ibin = emuefficiency->FindBin(tantheta); + const float efficiency = emuefficiency->GetBinContent(ibin); + if (prob < efficiency) return true; //efficiency larger than probability, we take the event + else return false; +} diff --git a/shipLHC/scripts/fromsndsw2FEDRA/track.rootrc b/shipLHC/scripts/fromsndsw2FEDRA/track.rootrc new file mode 100644 index 0000000000..5088945f95 --- /dev/null +++ b/shipLHC/scripts/fromsndsw2FEDRA/track.rootrc @@ -0,0 +1,30 @@ +fedra.readCPcut: 1 +fedra.track.TrZmap: 8000 0 400000 8000 0 400000 30 +fedra.track.npass: 1 +fedra.track.minPlate: -999 +fedra.track.maxPlate: 999 +fedra.track.refPlate: 999 +fedra.track.nsegmin: 2 +fedra.track.ngapmax: 4 +fedra.track.DZGapMax: 5000 +fedra.track.DRmax: 45 +fedra.track.DTmax: 0.07 +fedra.track.Sigma0: 3 3 0.005 0.005 +fedra.track.PulsRamp0: 15 20 +fedra.track.PulsRamp04: 15 20 +fedra.track.Degrad: 4 +fedra.track.probmin: 0.001 +fedra.track.momentum: 2 +fedra.track.mass: 0.14 +fedra.track.do_use_mcs: 0 +fedra.track.RadX0: 5810 +fedra.track.erase: 0 +fedra.track.do_realign: 0 +fedra.track.do_misalign: 0 +fedra.track.misalign_offset: 500 +fedra.track.do_local_corr: 1 +fedra.track.do_comb: 0 +fedra.track.NsegMin: 2 +emtra.outdir: .. +emtra.env: track.rootrc +emtra.EdbDebugLevel: 1 diff --git a/shipLHC/scripts/fromsndsw2FEDRA/track_experimental.rootrc b/shipLHC/scripts/fromsndsw2FEDRA/track_experimental.rootrc new file mode 100644 index 0000000000..db87473326 --- /dev/null +++ b/shipLHC/scripts/fromsndsw2FEDRA/track_experimental.rootrc @@ -0,0 +1,30 @@ +fedra.readCPcut: 1 +fedra.track.TrZmap: 8000 0 400000 8000 0 400000 30 +fedra.track.npass: 1 +fedra.track.minPlate: -999 +fedra.track.maxPlate: 999 +fedra.track.refPlate: 999 +fedra.track.nsegmin: 2 #unused +fedra.track.ngapmax: 4 #unused in the current framework +fedra.track.DZGapMax: 6000 +fedra.track.DRmax: 45 +fedra.track.DTmax: 0.04 +fedra.track.Sigma0: 4 4 0.003 0.003 +fedra.track.PulsRamp0: 15 20 +fedra.track.PulsRamp04: 15 20 +fedra.track.Degrad: 6 +fedra.track.probmin: 0.001 +fedra.track.momentum: 2 +fedra.track.mass: 0.14 +fedra.track.do_use_mcs: 0 +fedra.track.RadX0: 5810 +fedra.track.erase: 0 +fedra.track.do_realign: 0 +fedra.track.do_misalign: 0 +fedra.track.misalign_offset: 500 +fedra.track.do_local_corr: 1 +fedra.track.do_comb: 0 +fedra.track.NsegMin: 4 +emtra.outdir: .. +emtra.env: track.rootrc +emtra.EdbDebugLevel: 1 diff --git a/shipLHC/scripts/fromsndsw2FEDRA/vertexing.C b/shipLHC/scripts/fromsndsw2FEDRA/vertexing.C new file mode 100644 index 0000000000..ef706b72f7 --- /dev/null +++ b/shipLHC/scripts/fromsndsw2FEDRA/vertexing.C @@ -0,0 +1,262 @@ +//---------------------------------------------------------------------------- +// +// Usage: root -l vertexing.C\(nbrick\) +// with nbrick number of the brick which we are analyzing (11,12,13,14,21,22...51,52,53,54) +//---------------------------------------------------------------------------- + +R__EXTERN Int_t gEDBDEBUGLEVEL; + +EdbDataProc *dproc=0; +EdbPVRec *gAli=0; +EdbVertexRec *gEVR=0; +EdbDisplay *ds=0; + +int BRICKID = 0; +TH1D *hip = new TH1D("hip","Impact parameters",500,0,5000); + +void trvol( const char *def, const char *rcut = "nseg>1" ); +void init( const char *def, int iopt, const char *rcut="1" ); +void set_segments_dz(float dz=300.); +void do_propagation(); +void do_vertex(); +void td(); +void sd(); +void checkpatterns(); +void vd( int trmin, float amin); + +#include +namespace VERTEX_PAR +{ + float DZmax = 3000.; + float ProbMinV = 0.0001; // minimum acceptable probability for chi2-distance between tracks + float ImpMax = 15.; // maximal acceptable impact parameter [microns] (for preliminary check) + bool UseMom = false; // use or not track momentum for vertex calculations + bool UseSegPar = true; // use only the nearest measured segments for vertex fit (as Neuchatel) + int QualityMode= 0; // vertex quality estimation method (0:=Prob/(sigVX^2+sigVY^2); 1:= inverse average track-vertex distance) +} + +//--------------------------------------------------------------------- +void vertexing(int brID = 0, char* dset = 0) +{ + + gEDBDEBUGLEVEL=0; + //trvol(dset); + //trvol(dset,"nseg>1 && TMath::Abs(s.eY-50000)<5000 && TMath::Abs(s.eX-50000)<5000"); + BRICKID = brID; +// SELECTION WHICH TRACKS TO USE FOR VERTEXING + trvol(dset,"nseg>1"); + // trvol(dset,"nseg>1 && s.eX > 62500 && s.eY > 50000"); //fourth quarter + // reconstruct vertexes starting from linked_tracks.root +} + +//--------------------------------------------------------------------- +void trvol( const char *def, const char *rcut ) +{ + // this function read volume tracks and do the vertex reconstruction + // from linked_tracks.root + + init(def, 100 ,rcut); // read tracks (option 100) + gAli->FillCell(30,30,0.009,0.009); + do_vertex(); + //vd(VERTEX_CUTS::nmintracks,VERTEX_CUTS::angleaperture); // draw reconstructed vertex + //td(); //to draw tracks + //dproc->MakeTracksTree(gAli,"linked_tracks_p.root"); +} + +//--------------------------------------------------------------------- +void init( const char *def, int iopt, const char *rcut) +{ + if(!def) dproc = new EdbDataProc(); + else dproc = new EdbDataProc(def); + + dproc->InitVolume(iopt, rcut); + gAli = dproc->PVR(); + if (BRICKID > 0) checkpatterns(); + set_segments_dz(300.); +} + +//--------------------------------------------------------------------- +void set_segments_dz(float dz) +{ + int np = gAli->Npatterns(); + for(int i=0; iGetPattern(i); + int ns = p->N(); + for(int j=0; jGetSegment(j)->SetDZ(dz); + } +} + +void checkpatterns(){ + //check for patterns, if there is one missing add it + int np = gAli->Npatterns(); + + TFile *setfile = TFile::Open(Form("b%06d.0.0.0.set.root",BRICKID)); + EdbScanSet *set = (EdbScanSet*) setfile->Get("set"); + + for(int i=0; iGetPattern(i); + if(!p) { + cout<<"missing pattern "<eB.GetPlate(i)->Z(); //note, set->eB->GetPlate(i) uses the PID, set->GetPlate(i) uses the number of plate, be careful! + EdbPattern *pat = new EdbPattern( 0., 0.,zmissingPID); + pat->SetID(i); + pat->SetScanID(0); + gAli->AddPatternAt(pat,i); + }//end if + } //end for loop +} + +//--------------------------------------------------------------------- +void do_vertex() +{ + using namespace VERTEX_PAR; + + //gAli->FitTracks(4.,0.139 ); + gEVR = new EdbVertexRec(); + gEVR->eEdbTracks = gAli->eTracks; + gEVR->eVTX = gAli->eVTX; + gEVR->SetPVRec(gAli); + + gEVR->eDZmax=DZmax; + gEVR->eProbMin=ProbMinV; + gEVR->eImpMax=ImpMax; + gEVR->eUseMom=UseMom; + gEVR->eUseSegPar=UseSegPar; + gEVR->eQualityMode=QualityMode; + + printf("%d tracks vor vertexing\n", gEVR->eEdbTracks->GetEntries() ); + int nvtx = gEVR->FindVertex(); + printf("%d 2-track vertexes was found\n",nvtx); + + if(nvtx == 0) return; + int nadd = gEVR->ProbVertexN(); + + int nl = gEVR->LinkedVertexes(); //should avoid same track associated to multiple vertices + printf("%d vertices are linked\n",nl); + + EdbDataProc::MakeVertexTree(*(gEVR->eVTX),"vertextree.root"); +} + +//--------------------------------------------------------------------- +void td() +{ + // draw all tracks + + TObjArray *trarr=gAli->eTracks; + gStyle->SetPalette(1); + const char *dsname="display-t"; + ds = EdbDisplay::EdbDisplayExist(dsname); + if(!ds) ds=new EdbDisplay(dsname,-50000.,50000.,-50000.,50000.,-4000.,80000.); + ds->SetVerRec(gEVR); + ds->SetDrawTracks(4); + ds->SetArrTr( trarr ); + printf("%d tracks to display\n", trarr->GetEntries() ); + ds->Draw(); +} + +//--------------------------------------------------------------------- +void sd() +{ + // draw all tracks and segments (basetracks) + + TObjArray *trarr = gAli->eTracks; + TObjArray *sarr = new TObjArray(); + + EdbSegP *s=0; + for(int i=0; iNpatterns(); i++) { + EdbPattern *pat = gAli->GetPattern(i); + for(int j=0; jN(); j++) { + s = pat->GetSegment(j); + if(s->Track()<0) // exclude segments already attached to tracks + sarr->Add(s); + } + } + + printf("%d tracks to display\n", trarr->GetEntries() ); + printf("%d segments to display\n", sarr->GetEntries() ); + + gStyle->SetPalette(1); + const char *dsname="display-s"; + ds = EdbDisplay::EdbDisplayExist(dsname); + if(!ds) ds=new EdbDisplay(dsname,-50000.,50000.,-50000.,50000.,-4000.,80000.); + ds->SetVerRec(gEVR); + ds->SetDrawTracks(4); + ds->SetArrSegP( sarr ); + ds->SetArrTr( trarr ); + ds->Draw(); +} + +//--------------------------------------------------------------------- +void vd( int trmin, float amin) +{ + // draw vertexes with multiplicity>=trmin, and aperture >= amin + + TObjArray *varr = new TObjArray(); + TObjArray *tarr = new TObjArray(); + + EdbVertex *v=0; + EdbTrackP *t=0; + + int nv = gEVR->Nvtx(); + printf("nv=%d\n",nv); + if(nv<1) return; + + for(int i=0; ieVTX->At(i)); + if(v->Flag()<0) continue; + if( v->N()N()<3 && v->MaxAperture()Add(v); + for(int j=0; jN(); j++) tarr->Add( v->GetTrack(j) ); + } + + gStyle->SetPalette(1); + + const char *dsname="display-v"; + ds = EdbDisplay::EdbDisplayExist(dsname); + if(!ds) ds=new EdbDisplay(dsname,-100000.,100000.,-100000.,100000.,-40000., 0.); + ds->SetVerRec(gEVR); + ds->SetArrTr( tarr ); + printf("%d tracks to display\n", tarr->GetEntries() ); + ds->SetArrV( varr ); + printf("%d vertex to display\n", varr->GetEntries() ); + ds->SetDrawTracks(4); + ds->SetDrawVertex(1); + ds->Draw(); +} + +// Experimental function to study segment distance to vertex +void add_proton(EdbVertex* vertex, EdbVertexRec* gEVR){ +cout<<"start analyzing a vertex"<eEdbTracks->GetEntries(); //number of tracks +// EdbVTA *vta = 0; + EdbSegP *potsegment; + + bool associated_track; //track near the vertex + + for(int itr=0; itreEdbTracks->At(itr)); + + if (tr->N() > 26){ + + float iptrack2vertex = vertex->DistTrack(tr,0.); + if (iptrack2vertex < track_ipmin) track_ipmin = iptrack2vertex; + + if (iptrack2vertex < 200.) associated_track = true; //this is not the right track + + if (associated_track) cout<<"Let's see each segment for this track, track distance = "<N(); iseg++){ //loop on segments + + potsegment = (EdbSegP*)tr->GetSegment(iseg); + float ipseg2vertex = vertex->DistSeg(potsegment,0.); + if (associated_track) cout<Z()<Fill(track_ipmin); +}