FEDRA emulsion software from the OPERA Collaboration
EdbEDAIO Class Reference

#include <EdbEDA.h>

Inheritance diagram for EdbEDAIO:

Public Types

enum  OutputFileMode { kBern , kOtherLabs }
 

Public Member Functions

 EdbEDAIO ()
 
EdbPatternGetPatternIPL (int pid)
 
int IdMuon (char *filename="../cs_info.txt")
 
void OpenTextEditor (char *filename)
 
EdbPatternReadCouples (int ipl, EdbPattern *pat=NULL)
 
EdbPatternReadCouplesPID (int pid, EdbPattern *pat=NULL)
 
void ReadFeedbackFile (char *filename=NULL)
 
EdbVertexReadFeedbackFile2008 (char *filename=NULL)
 
EdbPVRecReadFilteredText (char *filename=NULL)
 
void WriteFeedbackFile (char *filename=NULL)
 
void WriteFeedbackFile2008 (EdbVertex *v=NULL, char *filename=NULL)
 
void WriteFilteredText (char *filename=NULL)
 
virtual ~EdbEDAIO ()
 

Public Attributes

TCut eCutCP
 
int eOutputFileMode
 

Member Enumeration Documentation

◆ OutputFileMode

Enumerator
kBern 
kOtherLabs 
165 {
166 kBern,
168 };
@ kOtherLabs
Definition: EdbEDA.h:167
@ kBern
Definition: EdbEDA.h:166

Constructor & Destructor Documentation

◆ EdbEDAIO()

EdbEDAIO::EdbEDAIO ( )
inline
172:eOutputFileMode(kBern), eCutCP("eN1==1&&eN2==1&&eCHI2P<s.eW*0.12-1"){}
int eOutputFileMode
Definition: EdbEDA.h:164
TCut eCutCP
Definition: EdbEDA.h:170

◆ ~EdbEDAIO()

virtual EdbEDAIO::~EdbEDAIO ( )
inlinevirtual
173{}

Member Function Documentation

◆ GetPatternIPL()

EdbPattern * EdbEDAIO::GetPatternIPL ( int  pid)

return Pattern with Couples.
if Couples are not loaded yet, load it in PVRec for TS.

262 {
265
267 EdbPVRec *pvr = set->GetPVRec();
268
269 EdbPattern *pat=NULL;
270
271 for(int i=0;i<pvr->Npatterns();i++){
272 EdbPattern *p = pvr->GetPattern(i);
273 if(p==NULL) continue;
274 if(p->N()>0){
275 if(p->GetSegment(0)->Plate()==ipl) {
276 pat=p;
277 break;
278 }
279 }
280 }
281
282 if(pat==NULL) {
283 printf("No pattern in TS, pl=%d.\n", ipl);
284 return NULL;
285 }
286
287 // Check if couples data were already loaded.
288 if( IsIncludeCouples(pat)==0 ) {
289
290 if (set->GetDataSet()!=NULL) gEDA->ReadCouples(ipl,pat);
291 else if(set->GetScanSet()!=NULL) {
292 printf("ReadCouples SS ipl=%d\n", ipl);
293 EdbID id = *(EdbID *)set->GetScanSet()->eIDS.At(0);
294 id.ePlate=ipl;
295 id.Print();
296 float z = pat->Z();
297 gEDA->ScanProc()->ReadPatCP(*pat,id,eCutCP);
298 pat->SetZ(z);
299 pat->SetSegmentsZ();
300 printf("z = %lf\n", pat->Z());
301 pat->SetSegmentsPlate(ipl);
302 }
303 }
304 return pat;
305}
EdbEDA * gEDA
Definition: EdbEDA.C:3
EdbPattern * ReadCouples(int ipl, EdbPattern *pat=NULL)
Definition: EdbEDA.C:218
Definition: EdbEDATrackSet.h:178
EdbEDATrackSet * GetTrackSet(int i)
Definition: EdbEDA.h:617
EdbScanProc * ScanProc()
Definition: EdbEDA.h:578
Definition: EdbID.h:7
Definition: EdbPVRec.h:148
Definition: EdbPattern.h:273
Int_t Npatterns() const
Definition: EdbPattern.h:366
EdbPattern * GetPattern(int id) const
Definition: EdbPattern.cxx:1721
int ReadPatCP(EdbPattern &pat, int id[4], TCut c="1")
Definition: EdbScanProc.cxx:2359
TList eIDS
list of the identifiers to be processed
Definition: EdbScanSet.h:16
Float_t Z() const
Definition: EdbPattern.h:84
void SetZ(float z)
Definition: EdbPattern.h:41
void SetSegmentsZ()
Definition: EdbPattern.cxx:250
void SetSegmentsPlate(int plate)
Definition: EdbPattern.cxx:241
EdbScanSet * set
Definition: emtraceback.cpp:14
int IsIncludeCouples(EdbPattern *pat)
Definition: EdbEDAUtil.C:41
#define NULL
Definition: nidaqmx.h:84
p
Definition: testBGReduction_AllMethods.C:8

◆ IdMuon()

int EdbEDAIO::IdMuon ( char *  filename = "../cs_info.txt")

code by GL

1045 {
1047 char buf[256];
1048 char key[256];
1049 char m1[12], m2[10], m3[10],m4[10],m5[10],m6[10];
1050 static int val=-2;
1051 char mu[3];
1052
1053 if(val>=0) return val;
1054 if(val==-1) return val;
1055
1056 FILE *fp=fopen(filename,"r");
1057 if (fp==NULL) {
1058 printf("ERROR open file: %s \n", filename);
1059 val = -1;
1060 return val;
1061 }
1062
1063 // printf( "Read affine transformation from file: %s\n", fname );
1064 while( fgets(buf,256,fp)!=NULL ) {
1065
1066 for( int i=0; i<(int)strlen(buf); i++ )
1067 if( buf[i]=='#' ) {
1068 buf[i]='0'; // cut out comments starting from #
1069 break;
1070 }
1071
1072 if( sscanf(buf,"%s",key)!=1 ) continue;
1073 if ( !strcmp(key,"Prediction") )
1074 {
1075 sscanf(buf+strlen(key),"%s %s %s %s %s %s ",m1,m2,m3,m4,m5,m6);
1076 mu[0] = m2[3];
1077 mu[1] = m2[4];
1078 mu[2] = m2[5];
1079 val = atol(mu);
1080 // cout << m2[3] << " " << m2[4] << " " << m2[5]<< " " << val << endl;
1081
1082 }
1083 if ( !strcmp(key,"No") )
1084 {
1085 val=-1;
1086 }
1087
1088 }
1089 fclose(fp);
1090
1091 return val;
1092}
const char filename[256]
Definition: RecDispNU.C:83
fclose(pFile)

◆ OpenTextEditor()

void EdbEDAIO::OpenTextEditor ( char *  filename)
1094 {
1095 TEveBrowser* browser = gEve->GetBrowser();
1096 browser->StartEmbedding(1);
1097 new TGTextEditor(filename, gClient->GetRoot());
1098 browser->StopEmbedding();
1099
1100 unsigned int i,k=0;
1101 for(i=0;i<strlen(filename);i++){
1102 if(filename[i]=='\\'||filename[i]=='/') k=i+1;
1103 }
1104 browser->SetTabTitle(filename+k, 1);
1105}

◆ ReadCouples()

EdbPattern * EdbEDAIO::ReadCouples ( int  ipl,
EdbPattern pat = NULL 
)
218 {
219 printf("ReadCouples DS ipl=%d\n", ipl);
220 return ReadCouplesPID(gEDA->GetPID(ipl), pat0);
221}
int GetPID(int ipl)
Definition: EdbEDA.h:276
EdbPattern * ReadCouplesPID(int pid, EdbPattern *pat=NULL)
Definition: EdbEDA.C:225

◆ ReadCouplesPID()

EdbPattern * EdbEDAIO::ReadCouplesPID ( int  pid,
EdbPattern pat = NULL 
)
225 {
226
228 EdbDataPiece *piece = dset->GetPiece(pid);
229 EdbLayer *l = piece->GetLayer(0);
230 EdbPattern *pat = new EdbPattern( 0.,0., l->Z());
231 pat->SetPID(pid);
232 piece->InitCouplesTree("READ");
233 piece->GetCPData_new(pat);
234 piece->CloseCPData();
235 pat->SetSegmentsZ();
236 pat->SetID(pid);
237 pat->SetSegmentsPID();
238 pat->Transform( l->GetAffineXY() );
239 pat->TransformA( l->GetAffineTXTY() );
240 pat->TransformShr( l->Shr() );
241
242 if(pat0==0) return pat;
243
244 EdbScanCond cond;
245
246 // copy segments.
247 for(int i=0;i<pat->N();i++){
248 EdbSegP *s = pat->GetSegment(i);
249 if(s->Plate()==0) s->SetPlate(gEDA->GetIPL(pid));
250
251 // fill COV for vertexing
252 s->SetErrors();
253 cond.FillErrorsCov(s->TX(), s->TY(), s->COV());
254
255 pat0->AddSegment(*s);
256 }
257 delete pat;
258 return pat0;
259}
EdbDataProc * dset
Definition: RecDispEX.C:9
Edb raw data unit (scanned plate) associated with run file.
Definition: EdbDataSet.h:26
EdbLayer * GetLayer(int id)
Definition: EdbDataSet.h:87
int InitCouplesTree(const char *mode="READ")
Definition: EdbDataSet.cxx:1239
int GetCPData_new(EdbPattern *pat, EdbPattern *p1=0, EdbPattern *p2=0, TIndex2 *trseg=0)
Definition: EdbDataSet.cxx:792
void CloseCPData()
Definition: EdbDataSet.cxx:91
OPERA emulsion data set.
Definition: EdbDataSet.h:144
int GetIPL(int PID)
Definition: EdbEDA.h:252
EdbDataSet * GetDataSet()
Definition: EdbEDA.h:244
Definition: EdbLayer.h:39
float Shr() const
Definition: EdbLayer.h:89
EdbAffine2D * GetAffineTXTY()
Definition: EdbLayer.h:120
EdbAffine2D * GetAffineXY()
Definition: EdbLayer.h:119
float Z() const
Definition: EdbLayer.h:77
void SetSegmentsPID()
Definition: EdbPattern.cxx:1455
void SetID(int id)
Definition: EdbPattern.h:309
void SetPID(int pid)
Definition: EdbPattern.h:310
virtual void Transform(const EdbAffine2D *a)
Definition: EdbVirtual.cxx:155
Definition: EdbScanCond.h:10
void FillErrorsCov(float tx, float ty, TMatrixD &cov)
Definition: EdbScanCond.cxx:161
Definition: EdbSegP.h:21
void TransformA(const EdbAffine2D *affA)
Definition: EdbPattern.cxx:324
void TransformShr(const float shr)
Definition: EdbPattern.cxx:341
Int_t N() const
Definition: EdbPattern.h:86
EdbSegP * GetSegment(int i) const
Definition: EdbPattern.h:66
int pid[1000]
Definition: m2track.cpp:13
s
Definition: check_shower.C:55

◆ ReadFeedbackFile()

void EdbEDAIO::ReadFeedbackFile ( char *  filename = NULL)
757 {
758 // Read feedback file format (ver 2009 Oct).
759
760 gEDA->Reset();
761
762
763 char FlagPart[][10]={"", "mu", "charm", "e", "e-pair"};
764 char FlagCS[][10] ={"", "cs"};
765
766 // if filename is not given, open file browser.
767 if(filename==NULL){
768 TGFileInfo *fi=new TGFileInfo;
769 fi->fIniDir = StrDup(".");
770 const char *filetypes[] = { "Feedback File", "*.feedback", "All files","*",0,0};
771 fi->fFileTypes=filetypes;
772 new TGFileDialog(gClient->GetRoot(), gEve?gEve->GetMainWindow():0, kFDOpen, fi);
773 filename = fi->fFilename;
774 }
775
776 if(NULL==filename) return;
777
778 FILE *fp = fopen(filename,"rt");
779 if(fp==NULL) {
780 ErrorMessage(Form("File open error : %s . stop.\n", filename));
781 return;
782 }
783
784 EdbPVRec *pvr = new EdbPVRec;
785 EdbScanCond cond;
786 EdbTrackP *t=0;
787 EdbEDAText *tx = 0;
788
789 char buf[256];
790
791 char asetname[10];
792 EdbEDAAreaSet *aset=NULL;
793
794 while(fgets(buf,sizeof(buf),fp)){
795 TString line=buf;
796 printf("%s", line.Data());
797
798 // track comment = ""
799 if(line.Contains("track comment")){
800 int istart = line.Index("\"");
801 int iend = line.Index("\"",1,istart+1, TString::kExact);
802 if(istart!=kNPOS && iend!=kNPOS && istart+1!=iend){
803 char buf[256];
804 strncpy(buf, line.Data()+istart+1, iend-istart-1);
805 buf[iend-istart-1]='\0';
806 tx->AddText(buf);
807 }
808 }
809
810 if(line.Contains("Area information :")){
811 aset = new EdbEDAAreaSet;
812 aset->SetDraw(kTRUE);
813 sscanf( line.Data()+line.Index(":")+1, "%s", asetname);
814 printf("Area infor %s\n", asetname);
815 }
816 if(aset){
817 if(line.Contains(Form("area : %s", asetname))){
818 int ipl;
819 double xmin, xmax, ymin, ymax, z;
820 sscanf(line.Data() + line.Index(Form("area : %s", asetname)) + strlen(Form("area : %s", asetname)),
821 "%d %lf %lf %lf %lf %lf", &ipl, &xmin, &xmax, &ymin, &ymax, &z);
822 aset->AddArea(ipl, xmin, xmax, ymin, ymax, z);
823 }
824 }
825
826 // Remove comments
827 line.ReplaceAll("\t"," ");
828 int iposcmt = line.Index("//");
829 if(iposcmt!=kNPOS) line.Remove(iposcmt, line.Length()-iposcmt);
830
831 // Check number of columns.
832 TObjArray *tokens = line.Tokenize(" ");
833 int ntokens = tokens->GetEntries();
834 delete tokens;
835
836
837 if(ntokens==0) continue;
838
839 else if( ntokens == 10 ){
840 // Vertices
841 float x,y,z; int id,isprimary,istau, ischarm;
842 sscanf(line.Data(),"%d %f %f %f %d %d %d", &id, &x, &y, &z, &isprimary, &ischarm, &istau);
843 EdbVertex *v = new EdbVertex();
844 v->SetXYZ(x,y,z); v->SetID(id);
845 v->SetFlag(isprimary);
846 pvr->AddVertex(v);
847 printf("Vertex %d %f %f %f\n",v->ID(), v->X(), v->Y(), v->Z());
848 }
849
850 else if( ntokens == 27 ){
851 // Tracks
852 float x,y,z,ax,ay, ip_upvtx, ip_downvtx, p,pmin,pmax;
853 int id_track, id_upvtx, id_downvtx, manual, particle_id, scanback, darkness;
854 int OutOfBrick, LastPlate, nseg, iplRmax1, iplRmax2, result;
855 float RmaxT, RmaxL, rmst, rmsl;
856
857 int n,nn;
858 sscanf(line.Data(), "%d %d %d %f %f %f %f %f %f%n", &id_track, &id_upvtx, &id_downvtx, &x, &y, &z, &ax, &ay, &ip_upvtx, &nn);
859 sscanf(line.Data()+(n=nn), "%f %f %f %f %d %d %d %d %d%n", &ip_downvtx, &p,&pmin,&pmax, &manual, &particle_id, &scanback, &darkness, &OutOfBrick, &nn);
860 sscanf(line.Data()+(n+=nn), "%d %d %f %f %f %f %d %d %d", &LastPlate, &nseg, &RmaxT, &RmaxL, &rmst, &rmsl, &iplRmax1, &iplRmax2, &result);
861
862 // Create Track. "t" is defined out of main loop.
863 t = new EdbTrackP;
864 t->Set(id_track, x, y, ax, ay, 0, 0);
865 t->SetZ(z);
866 t->SetTrack(id_track);
867 pvr->AddTrack(t);
868
869 tx = new EdbEDAText(t->X(), t->Y(), t->Z(), "[");
870 if(0<=particle_id && particle_id<=4) { if(tx->N()>1)tx->AddText(", "); tx->AddText(FlagPart[particle_id]);}
871 if(scanback==1) { if(tx->N()>1)tx->AddText(", "); tx->AddText(FlagCS[scanback]);}
872 tx->AddText("]");
873 if(tx->N()==2) tx->SetText("");
874 tx->SetReference(tx);
875 gEDA->AddDrawObject(tx);
876
877 // fill COV for vertexing
878 t->SetErrors();
879 cond.FillErrorsCov(t->TX(), t->TY(), t->COV());
880
881 // associate to vertex.
882 for(int i=0; i<pvr->Nvtx(); i++){
883 EdbVertex *v = pvr->GetVertex(i);
884 if(id_upvtx==v->ID()||id_downvtx==v->ID()){
885 EdbVTA *vta = new EdbVTA(t, v);
886 vta->SetZpos( t->Z()>v->Z() ? 1 : 0);
887 vta->SetFlag(2);
888 vta->SetImp( id_upvtx==v->ID()? ip_upvtx : ip_downvtx);
889 v->AddVTA(vta);
890 }
891 }
892 }
893
894 else if( ntokens == 9 ){
895 // Segments
896 int ipl, type, irec, ph;
897 float x,y,z,ax,ay;
898 sscanf(line.Data(),"%d %f %f %f %f %f %d %d %d", &ipl, &x, &y, &z, &ax, &ay, &type, &irec, &ph);
899
900 EdbSegP *s = new EdbSegP(t->ID(),x,y,ax,ay,0,0);
901 s->SetZ(z);
902 s->SetPID(ipl);
903 s->SetPlate(ipl);
904 s->SetW(ph);
905 s->SetTrack(t->ID());
906 s->SetSide(type);
907
908 // fill COV for vertexing
909 s->SetErrors();
910 cond.FillErrorsCov(s->TX(), s->TY(), s->COV());
911
912 // Add segment in PVRec and Track, keeping consistency of pointer in them.
913 EdbSegP *s_in_pattern = pvr->AddSegment(*s);
914 t->AddSegment( s_in_pattern);
915
916 if(NULL!=tx){
917 EdbSegP seg(*s);
918 seg.PropagateTo(seg.Z()+1000);
919 tx->SetXYZ(seg.X(),seg.Y(),seg.Z());
920 }
921 }
922
923
924 }
925
926 for(int i=0;i<pvr->Npatterns(); i++) pvr->GetPattern(i)->SetPID(i);
927 for(int i=0;i<pvr->Npatterns(); i++) pvr->GetPattern(i)->SetSegmentsPID();
928 for(int i=0;i<pvr->Ntracks(); i++) pvr->GetTrack(i)->SetCounters();
929
930 printf("\nEdbPVRec summary. %d vertices, %d tracks.\n", pvr->Nvtx(), pvr->Ntracks());
931 pvr->Print();
932
933 fclose(fp);
934
935 if(pvr==NULL) return;
936 for(int i=0;i<pvr->Npatterns(); i++){
937 EdbPattern *pat = pvr->GetPattern(i);
938 pat->SetPID(gEDA->GetPID(pat->ID()));
939 }
940
941 gEDA->SetPVR(pvr);
942 gEDA->GetTrackSet("TS")->AddTracksPVR(pvr);
944
945 if(aset){
946 EdbEDATrackSet *set = gEDA->GetTrackSet(asetname);
947 if(set) set->SetAreaSet(aset);
948 }
949 gEDA->Redraw();
950 gEDA->UpdateScene();
951}
Definition: EdbEDASets.h:55
void AddArea(int ipl, double xmin, double xmax, double ymin, double ymax, double z)
Definition: EdbEDASets.h:110
void SetDraw(bool draw)
Definition: EdbEDASets.h:106
void SetPVR(EdbPVRec *pvr)
Definition: EdbEDA.h:224
void SetReference(void *ref)
Definition: EdbEDA.h:43
Definition: EdbEDA.h:50
void SetXYZ(float x, float y, float z)
Definition: EdbEDA.h:76
int N()
Definition: EdbEDA.h:79
void AddText(TString str)
Definition: EdbEDA.h:78
void SetText(TString str)
Definition: EdbEDA.h:77
void AddTracksPVR(EdbPVRec *pvr)
Definition: EdbEDATrackSet.h:312
void AddVertices(TObjArray *vertices)
Definition: EdbEDASets.h:189
EdbEDAVertexSet * GetVertexSet()
Definition: EdbEDA.h:660
void UpdateScene()
Definition: EdbEDA.h:681
void Redraw()
Definition: EdbEDA.h:680
void AddDrawObject(EdbEDAObject *o)
Definition: EdbEDA.h:582
void Reset()
Definition: EdbEDA.h:683
EdbVertex * GetVertex(Int_t &i)
Definition: EdbPVRec.h:256
TObjArray * eVTX
array of vertex
Definition: EdbPVRec.h:162
Int_t Ntracks() const
Definition: EdbPVRec.h:203
void AddVertex(EdbVertex *vtx)
Definition: EdbPVRec.h:251
Int_t Nvtx() const
Definition: EdbPVRec.h:255
EdbSegP * AddSegment(EdbSegP &s)
Definition: EdbPVRec.cxx:2925
EdbTrackP * GetTrack(int i) const
Definition: EdbPVRec.h:241
void AddTrack(EdbTrackP *track)
Definition: EdbPVRec.h:246
int ID() const
Definition: EdbPattern.h:319
void Print() const
Definition: EdbPattern.cxx:1693
Definition: EdbPattern.h:113
void SetCounters()
Definition: EdbPattern.h:159
Definition: EdbVertex.h:26
void SetImp(float imp)
Definition: EdbVertex.h:57
void SetZpos(int zpos)
Definition: EdbVertex.h:55
void SetFlag(int flag)
Definition: EdbVertex.h:56
Definition: EdbVertex.h:69
Int_t ID() const
Definition: EdbVertex.h:126
void AddVTA(EdbVTA *vta)
Definition: EdbVertex.cxx:355
Float_t X() const
Definition: EdbVertex.h:130
void SetID(int ID=0)
Definition: EdbVertex.h:156
void SetFlag(int flag=0)
Definition: EdbVertex.h:158
void SetXYZ(float x, float y, float z)
Definition: EdbVertex.h:157
Float_t Z() const
Definition: EdbVertex.h:132
Float_t Y() const
Definition: EdbVertex.h:131
TTree * t
Definition: check_shower.C:4
void ErrorMessage(char *title, char *message)
Definition: EdbEDAUtil.C:479
UInt_t id
Definition: tlg2couples.C:117
Int_t type
Definition: testBGReduction_By_ANN.C:15

◆ ReadFeedbackFile2008()

EdbVertex * EdbEDAIO::ReadFeedbackFile2008 ( char *  filename = NULL)
954 {
955 using namespace std;
956
957 // if filename is not given, open file browser.
958 if(filename==NULL){
959 TGFileInfo *fi=new TGFileInfo;
960 fi->fIniDir = StrDup(".");
961 const char *filetypes[] = { "Feedback File", "*.feedback", "All files","*",0,0};
962 fi->fFileTypes=filetypes;
963 new TGFileDialog(gClient->GetRoot(), gEve->GetMainWindow(), kFDOpen, fi);
964 filename = fi->fFilename;
965 }
966
967 if(NULL==filename) return NULL;
968
969 // check the file.
970 ifstream fin(filename);
971 if(!fin) {
972 printf("can not open file : %s\n", filename);
973 return NULL;
974 } printf("read feedback %s\n", filename);
975
976
978
980 set->ClearTracks();
981
982 // read file line by line.
983 char buf[512];
984 while( fin && fin.getline(buf,512) ){
985 string line(buf);
986
987 // ignore comment
988 if(line.find("//")!=line.npos) line[line.find("//")]='\0';
989
990 // check the number of columnline.
991 istringstream ss(line);
992 string s2;
993 int ncolumn;
994 for(ncolumn=0;ss>>s2;ncolumn++){
995 //cout <<s2 << endl;
996 }
997
998 // ncolumn == 7 --> vertex info. ncolumn == 17 --> track info
999 if(ncolumn==7){
1000 double x,y,z;
1001 int id,isprimary,istau, ischarm;
1002 sscanf(line.c_str(),"%d %lf %lf %lf %d %d %d", &id, &x, &y, &z, &isprimary, &ischarm, &istau);
1003 EdbVertex *v = new EdbVertex();
1004 v->SetXYZ(x,y,z);
1005 v->SetID(id);
1006 gEDA->AddVertex(v);
1007 }
1008 else if (ncolumn==17){
1009 double x,y,z,ax,ay, ip_upvtx, ip_downvtx;
1010 double p,pmin,pmax;
1011 int id_track, id_upvtx, id_downvtx;
1012 int manual, particle_id, scanback, darkness;
1013
1014 istringstream ss(line);
1015
1016 ss >> id_track >> id_upvtx >> id_downvtx >> x >> y >> z >> ax >> ay >> ip_upvtx >> ip_downvtx
1017 >> p >> pmin >> pmax >> manual >> particle_id >> scanback >> darkness;
1018
1019 EdbSegP s(id_track, x, y, ax, ay);
1020 s.SetZ(z);
1021
1022 EdbTrackP *t = set->FindTrack(&s);
1023
1024 if(t) {
1025 printf("%s track %d\n", line.c_str(), t->ID());
1026 set->SetTrack(t);
1027 for(int i=0; i<gEDA->GetVertexSet()->N(); i++){
1029 if(id_upvtx!=v->ID()) continue;
1030 EdbVTA *vta = new EdbVTA(t, v);
1031 vta->SetZpos(1);
1032 vta->SetFlag(2);
1033 vta->SetImp(EdbEDAUtil::CalcIP(&s,v));
1034 v->AddVTA(vta);
1035 }
1036 }
1037 else printf("track not found. id = %d\n", id_track);
1038 }
1039 }
1040 gEDA->Redraw();
1041 if( gEDA->GetVertexSet()->N()) return gEDA->GetVertexSet()->GetVertex(0);
1042 else return NULL;
1043}
int N()
Definition: EdbEDASets.h:215
EdbVertex * GetVertex(int i)
Definition: EdbEDASets.h:216
void AddVertex(EdbVertex *v)
Definition: EdbEDA.h:661
void ClearVertices()
Definition: EdbEDA.h:663
ss
Definition: energy.C:62
EdbSegP * s2
Definition: tlg2couples.C:30
double CalcIP(EdbSegP *s, double x, double y, double z)
Definition: EdbEDAUtil.C:85
Definition: AlignmentCint.cxx:51

◆ ReadFilteredText()

EdbPVRec * EdbEDAIO::ReadFilteredText ( char *  filename = NULL)
1108 {
1109 int i,j;
1110
1111 if(filename==NULL||strlen(filename)<1){ // File dialog
1112 TGFileInfo *fi=new TGFileInfo;
1113 fi->fIniDir = StrDup(".");
1114 const char *filetypes[] = { "Filtered text", "*.txt", "All files","*",0,0};
1115 fi->fFileTypes = filetypes;
1116 new TGFileDialog(gClient->GetRoot(), gEve->GetMainWindow(), kFDOpen, fi);
1117 if(fi->fFilename!=NULL) {
1118 filename = fi->fFilename;
1119 }
1120 }
1121 if(filename==NULL) return NULL;
1122
1123 // check the file.
1124 ifstream fin(filename);
1125 if(!fin) {
1126 printf("can not open file : %s\n", filename);
1127 return NULL;
1128 }
1129 else printf("Read Filtered Text File : %s\n", filename);
1130
1131
1132 TObjArray *vtxarr = new TObjArray;
1133 TObjArray *trkarr = new TObjArray;
1134 TObjArray *segarr = new TObjArray;
1135
1136
1137 EdbPVRec *PVR = new EdbPVRec;
1138
1140
1141 for(i=0;i<dset->N();i++){
1142 EdbDataPiece *piece = dset->GetPiece(i);
1143 EdbLayer* l = piece->GetLayer(0);
1144 EdbPattern *pat = new EdbPattern( 0.0, 0.0, l->Z());
1145 pat->SetPID(i);
1146 PVR->AddPattern(pat);
1147 }
1148
1149 char buf[512];
1150 while( fin && fin.getline(buf,512) ){
1151 string line(buf);
1152
1153 std::istringstream ssline(line);
1154 std::string column1;
1155 ssline >> column1;
1156
1157 std::cout <<line <<std::endl;
1158 if(column1.find("//")!=column1.npos) continue;
1159
1160 if(column1=="vertex"){
1161 int id;
1162 float x, y, z;
1163 ssline >> id >> x >> y >> z;
1164 EdbVertex *v = new EdbVertex;
1165 v->SetXYZ(x,y,z);
1166 v->SetID(id);
1167 vtxarr->Add(v);
1168 PVR->AddVertex(v);
1169 }
1170
1171 if(column1=="track"){
1172 int id, flag, ipl1, ipl2, nseg, ivtx1, ivtx2;
1173 double tx,ty, ip_1ry, ip1, ip2;
1174 ssline >> id >> flag >> tx >> ty >> ipl1 >> ipl2 >> nseg >> ip_1ry >> ivtx1 >> ip1 >> ivtx2 >> ip2;
1175 EdbTrackP *t = new EdbTrackP;
1176 t->SetID(id);
1177 t->SetFlag(flag);
1178 for(i=0;i<vtxarr->GetEntries();i++){
1179 EdbVertex *v = (EdbVertex *) vtxarr->At(i);
1180 if(v->ID()==ivtx1||v->ID()==ivtx2){
1181 EdbVTA *vta = new EdbVTA(t, v);
1182 vta->SetFlag(2);
1183 v->AddVTA(vta);
1184 }
1185 }
1186 trkarr->Add(t);
1187 PVR->AddTrack(t);
1188 }
1189
1190 if(column1=="segment"){
1191 int itrk, ipl, iseg, sflag, tflag;
1192 float x, y, z, ax, ay, chi2, w;
1193 ssline >> itrk >> tflag >> ipl >> iseg >> x >> y >> z >> ax >> ay >> w >> chi2 >> sflag;
1194
1195 EdbSegP *s = new EdbSegP (iseg, x, y, ax, ay, w, sflag);
1196 s->SetPID(gEDA->GetPID(ipl));
1197 s->SetZ(z);
1198 s->SetDZ(300);
1199 s->SetTrack(itrk);
1200 s->SetChi2(chi2);
1201
1202 for(i=0;i<trkarr->GetEntries();i++){
1203 EdbTrackP *t = (EdbTrackP *) trkarr->At(i);
1204 if(t->ID()==itrk&&t->Flag()==tflag) {
1205 t->AddSegment(s);
1206 }
1207 }
1208
1209 EdbPattern *pat = PVR->GetPattern( gEDA->GetPID(ipl));
1210 pat->AddSegment(*s);
1211 segarr->Add(s);
1212
1213 }
1214
1215 if(column1=="area"){
1216 int ipl;
1217 double xmin,xmax,ymin,ymax,z;
1218 ssline >> ipl >> xmin >> xmax >> ymin >> ymax >> z;
1219 EdbEDAArea *a = new EdbEDAArea(ipl, xmin, xmax, ymin, ymax, z);
1220 gEDA->GetTrackSet("TS")->GetAreaSet()->AddArea(a);
1221 }
1222
1223
1224 }
1225 fin.close();
1226
1227 printf("Vertices %4d\n", vtxarr->GetEntries());
1228 printf("Tracks %4d\n", trkarr->GetEntries());
1229 printf("Segments %4d\n", segarr->GetEntries());
1230
1231
1232 // Clear all track set
1233 for(i=0; i<gEDA->NTrackSets();i++){
1234 gEDA->GetTrackSet(i)->Clear();
1235 }
1236
1237 for(i=0;i<trkarr->GetEntries();i++){
1238 EdbTrackP *t = (EdbTrackP *) trkarr->At(i);
1239 t->SetCounters();
1240 EdbSegP *s = t->GetSegmentFirst();
1241 t->Set(t->ID(), s->X(), s->Y(), s->TX(), s->TY(), 0, t->Flag());
1242 t->SetZ(s->Z());
1243
1244 if( 0 <= t->Flag()) gEDA->GetTrackSet("TS")->AddTrack(t);
1245 if( -100 < t->Flag() && t->Flag() < 0) gEDA->GetTrackSet("MN")->AddTrack(t);
1246 if( -200 < t->Flag() && t->Flag() <= -99) gEDA->GetTrackSet("SB")->AddTrack(t);
1247 if(-1000 < t->Flag() && t->Flag() <= -200) gEDA->GetTrackSet("SF")->AddTrack(t);
1248 }
1249
1252 vset->Clear();
1253 for(i=0;i<vtxarr->GetEntries();i++){
1254 EdbVertex *v = (EdbVertex *) vtxarr->At(i);
1255
1256 for(j=0;j<v->N();j++){
1257 EdbVTA *vta = v->GetVTa(j);
1258 EdbTrackP *t = vta->GetTrack();
1259 EdbSegP *s = t->GetSegmentFirst();
1260 vta->SetZpos( s->Z()>v->Z() ? 1 : 0 );
1261 vta->SetImp(CalcIP(s,v));
1262 }
1263 vset->AddVertex(v);
1264 }
1265
1266 return PVR;
1267}
void a()
Definition: check_aligned.C:59
Double_t CalcIP(EdbSegP *s, EdbVertex *v)
Definition: ShowRec.cpp:8872
Definition: EdbEDASets.h:7
void ClearSelectedVertex()
Definition: EdbEDA.h:441
void AddTrack(EdbTrackP *t)
Definition: EdbEDATrackSet.h:277
EdbEDAAreaSet * GetAreaSet()
Definition: EdbEDATrackSet.h:602
void Clear()
Definition: EdbEDATrackSet.h:336
Definition: EdbEDASets.h:143
void AddVertex(EdbVertex *v)
Definition: EdbEDASets.h:183
void Clear()
Definition: EdbEDASets.h:214
int NTrackSets()
Definition: EdbEDA.h:616
void AddPattern(EdbPattern *pat)
Definition: EdbPattern.cxx:1707
EdbSegP * AddSegment(int i, EdbSegP &s)
Definition: EdbPattern.cxx:72
EdbTrackP * GetTrack() const
Definition: EdbVertex.h:51
EdbVTA * GetVTa(int i)
Definition: EdbVertex.h:139
Int_t N() const
Definition: EdbVertex.h:121
Float_t chi2
Definition: testBGReduction_By_ANN.C:14
void w(int rid=2, int nviews=2)
Definition: test.C:27

◆ WriteFeedbackFile()

void EdbEDAIO::WriteFeedbackFile ( char *  filename = NULL)

write feedback file
if eOutputFileMode == 0 (kBern), output comments with "//" escape sequence (default).
also track id is not 1,2,3,4, but id of track.
if eOutputFileMode == 1 (kOtherLabs). doesn't output comments.

to convert Bern format, use
gawk -f ConvertFeedback.awk tmp.feedback

– ConvertFeedback.awk –
BEGIN{ itrk=1;}
!/\/\//&&NF==27{print itrk++,$2,$3,$4,$5,$6,$7,$8,$9,$10,$11,$12,$13,$14,$15,$16,$18,$19,$20,$21,$22,$23,$24,$25,$26,$27;}
!/\/\//&&NF==9{print $1,$2,$3,$4,$5,$6,$7,$8,$9;}
!/\/\//&&NF==10{print $1,$2,$3,$4,$5,$6,$7,$8,$9,$10;}

308 {
322
323
324 if( gEDA->GetVertexSet()->N() == 0 ){
325 ErrorMessage("No vertex is found. stop.");
326 return;
327 }
328
329 if(filename==NULL){ // File dialog
330 TGFileInfo *fi=new TGFileInfo;
331 fi->fIniDir = StrDup(".");
332 const char *filetypes[] = { "Feedback file", "*.feedback", "All files","*",0,0};
333 fi->fFileTypes = filetypes;
334 new TGFileDialog(gClient->GetRoot(), gEve->GetMainWindow(), kFDSave, fi);
335 if(fi->fFilename!=NULL) filename = fi->fFilename;
336 }
337
338 FILE *fp;
339 if(NULL==filename||strlen(filename)==0) {
340 return;
341 }
342 else {
343 fp=fopen(filename,"wt");
344 if(NULL==fp){
345 printf("Couldn't open file : %s\n", filename);
346 return;
347 }
348 }
349
350 printf("--------- Feedback format---------------\n");
351 printf("VID X Y Z 1ry charm tau Ndw Nup OutOfBrick\n");
352 printf("trkid v1 v2 x y z tx ty ip1 ip2 ");
353 printf(" p pmin pmax mn pa sb dk ou ps ");
354 printf(" n RmaxT RmaxL rmsT rmsL pl1 pl2 res\n");
355 printf(" ipl x y z tx ty type irec ngrains\n");
356
358 fprintf(fp,"// --------- Feedback format---------------\n");
359 fprintf(fp,"// VID X Y Z 1ry charm tau Ndw Nup OutOfBrick\n");
360 fprintf(fp,"// OutOfBrick : 1= vertex in dead material, 0= others\n");
361 fprintf(fp,"// \n");
362 fprintf(fp,"// Trkid v1 v2 x y z tx ty ip1 ip2 p pmin pmax mn pa sb dk ou ps n RmaxT RmaxL rmsT rmsL pl1 pl2 res\n");
363 fprintf(fp,"// Manual : 0= auto, 1= manually added\n");
364 fprintf(fp,"// Particle : 1= muon, 2= charm, 3= electron, 4= e+e- pair\n");
365 fprintf(fp,"// Scanback : 1= scanback, 0=not scanback\n");
366 fprintf(fp,"// Darkness : 0= MIP, 1= BLACK, 2= GRAY\n");
367 fprintf(fp,"// OutOfBrick : 1= pass through up, 2= edge out, 0 others\n");
368 fprintf(fp,"// Last plate : 1= pass through down, n= edge out last plate number, 0 others\n");
369 fprintf(fp,"// Out_flag(res): 1= 1ry vertex track, 2= e+e- pair, 3= Low momentum (DTheta_RMS >20 mrad),\n");
370 fprintf(fp,"// 4= Scan Forth to be done, 5= Scan Forth done\n");
371 fprintf(fp,"// \n");
372 fprintf(fp,"// Ipl x y z tx ty type irec ngrains\n");
373 fprintf(fp,"// Type : 0 BT, 1 microtrack top, 2 microtrack bottom\n");
374 fprintf(fp,"// Irec : 0 Automatic, 1 SB, 2 Manual\n");
375 fprintf(fp,"// \n");
376 fprintf(fp,"///////////////////////////////////////////// \n");
377 fprintf(fp,"// ---------- Vertex Data ---------------- // \n");
378 fprintf(fp,"///////////////////////////////////////////// \n");
379 }
380 // Vertices
381
382 TObjArray *vertices = (gEDA->GetVertexSet())->GetVertices();
383
384 printf("%d vertices\n", vertices->GetEntries());
385 EdbVertex *v_1ry=NULL;
386 double z=1e9;
387 for(int i=0;i<vertices->GetEntries();i++){
388 EdbVertex *v = (EdbVertex *) vertices->At(i);
389 if(v->Z()<z) {
390 z = v->Z();
391 v_1ry=v;
392 }
393 }
394 int vtxpl = gEDA->GetIPLZ(v_1ry->Z());
395 printf("Vertex plate = %d. z=%7.1lf\n", vtxpl, gEDA->GetZ(vtxpl));
396 if(eOutputFileMode==kBern) fprintf(fp, "// primary vertex plate = %d. dz = %f\n", vtxpl, gEDA->GetZ(vtxpl)-v_1ry->Z());
397
398 if(eOutputFileMode==kBern) fprintf(fp, "// VID X Y Z 1ry charm tau Ndw Nup OutOfBrick\n");
399
400 for(int i=0; i<vertices->GetEntries(); i++){
401 EdbVertex *v = (EdbVertex *) vertices->At(i);
402 v->SetID(i+1);
403 int ISPRIMARY = v==v_1ry ? 1 : 0;
404 int ISCHARM = 0;
405 int ISTAU = 0;
406 int Nup=0, Ndown=0;
407 for(int j=0;j<v->N();j++){
408 EdbTrackP *t = v->GetTrack(j);
409 if(t->GetSegmentFirst()->Z()>v->Z()) Ndown++;
410 if(t->GetSegmentLast()->Z()<v->Z()) Nup++;
411 }
412 int OutOfBrick = v->Z()<gEDA->GetZ(1)? 1:0; // 1: if vertex is outside of the brick. 0: others.
413
415 fprintf(fp,"%3d %8.1lf %8.1lf %8.1lf %3d %5d %3d %3d %3d %3d\n",
416 v->ID(), v->X(), v->Y(), v->Z(),ISPRIMARY,ISCHARM,ISTAU, Ndown, Nup, OutOfBrick);
417 } else {
418 fprintf(fp,"%d %8.1lf %8.1lf %8.1lf %3d %5d %3d %3d %3d %3d\n",
419 v->ID(), v->X(), v->Y(), v->Z(),ISPRIMARY,ISCHARM,ISTAU, Ndown, Nup, OutOfBrick);
420 }
421 }
422
423 // Tracks + Segments
424 int itrk=1;
425 for(int i=0;i<gEDA->NTrackSets();i++){
427 if(set->N()==0) continue;
428 if(set->GetDraw()==0) continue;
429 printf("TrackSet %s : %3d tracks\n", set->GetName(), set->N());
431 EdbID& id=set->GetID();
432 fprintf(fp,"////////////////////////////////////////////////////////// \n");
433 fprintf(fp,"// ---------- TrackSet : %s %6d 0 %d %4d ----------- //\n", set->GetName(), id.eBrick, id.eMajor, id.eMinor);
434 fprintf(fp,"////////////////////////////////////////////////////////// \n");
435 }
436 for(int j=0;j<set->N();j++){
437 EdbTrackP *t = set->GetTrack(j);
438 EdbSegP *sf = t->GetSegmentFirst();
439 EdbSegP *sl = t->GetSegmentLast();
440 EdbVertex * v1 = NULL;
441 EdbVertex * v2 = NULL;
442 for(int k=0;k<vertices->GetEntries();k++){
443 EdbVertex *v = (EdbVertex *) vertices->At(k);
444 for(int l=0;l<v->N();l++){
445 EdbTrackP *tv = v->GetTrack(l);
446 if(t==tv){
447 if(NULL==v1) v1 = v;
448 else {
449 if(v1->Z()>v->Z()) { // v1 should 1ry.
450 v2 = v1;
451 v1 = v;
452 }
453 else v2 = v;
454 } } } }
455 if(v1==NULL) v1=v_1ry; // all tracks are attached to 1ry vertex at this moment.
456
457 // t : pointer to the track
458 // v1 : upstream vertex
459 // v2 : downstream vertex
460 // sf : first segment to be used for vertex calculation
461 // fp : file pointer to print data.
462
463 int ivtx1 = v1 ? v1->ID() : -1;
464 int ivtx2 = v2 ? v2->ID() : -1;
465
466 float ip1= v1? CalcIP(sf,v1) : -1.0; // v1->DistSeg(sf): -1.0;
467 float ip2= v2? CalcIP(sl,v2) : -1.0; // v2->DistSeg(sl): -1.0;
468
469 double p,pmin,pmax;
470 CalcP(t, p, pmin, pmax); // pmin = pmin, pmax = pmax
471
472 // Flags.... same as Feedback file.
473 int MAN = 0;
474 int SCANBACK=0, idscanback=-1;
475 int MUON =0;
476 int DARKNESS=0; // no function to define.
477
478 if(gEDA->GetTrackSet("MN")->FindTrack(t)) MAN=1;
479
480 EdbTrackP * sbt = set==gEDA->GetTrackSet("SB")? t : gEDA->CheckScanback(t); // check if the track is a scanback track
481 if(sbt!=NULL) idscanback = sbt->ID();
482 if (idscanback!=-1) SCANBACK = 1;
483 if (idscanback==gEDA->IdMuon()&&gEDA->IdMuon()!=-1) MUON =1;
484
485 int OutOfBrick=0; // OutOfBrick: 1= pass through, 2=edge out, 0 others
486 int LastPlate=0; // Last plate: 1= pass though, n=edge out track last plate, 0 others
487
488 int nseg = t->N();
489
490 // Small Kink search
491
492 // remove Micro-tracks.
493 double RmaxT=0.0, RmaxL=0.0, Rmax=0.0, dtmaxt=0.0, dtmaxl=0.0, rmstmax=-1.0, rmslmax=-1.0;
494 int iplRmax1=0, iplRmax2=0;
495
496 EdbTrackP *tt = new EdbTrackP;
497 tt->Copy(*t);
498 tt->SetCounters();
499
500 for(int k=0;k<tt->N()&&tt->N()>1;k++){
501 EdbSegP *s = tt->GetSegment(k);
502 if(s->Side()!=0) tt->RemoveSegment(s); // if it's microtrack, remove.
503 }
504
505 int n = tt->N();
506
507 if(n<=1) {
508 // cannot calculate rms.
509 RmaxT=RmaxL=-1.;
510 iplRmax1=iplRmax2=-1;
511 }
512 else if( n<4 ){
513 // if nseg <=3, cannot calculate kink, just calculate rms but not R.
514 RmaxT=RmaxL=-1.;
515 iplRmax1=iplRmax2=-1;
516 }
517 else {
518 // Calculate rms and R.
519 for(int k=0; k<n-1;k++){
520 EdbSegP *ss1 = tt->GetSegment(k);
521 EdbSegP *ss2 = tt->GetSegment(k+1);
522
523 int ipl1 = ss1->Plate();
524 int ipl2 = ss2->Plate();
525
526 // if ss2 is on pl 57, skip.
527 if(ipl2>=gEDA->GetLastPlate()) continue;
528 // if ss1 is more downstream than 3 plates from vertex, skip. search only kink between 1-x,2-x,3-x.
529 if(ipl1>vtxpl+2) continue;
530
531 double dtt,dtl;
532 CalcDTTransLongi(ss1, ss2, &dtt, &dtl);
533 dtt=fabs(dtt);
534 dtl=fabs(dtl);
535
536 double r,rt,rl;
537 int ndata;
538 DTRMSTLGiven1Kink(tt, k, &r, &rt, &rl, &ndata);
539
540 double RT =dtt/rt/sqrt((double)ipl2-ipl1);
541 double RL =dtl/rl/sqrt((double)ipl2-ipl1);
542 double R=sqrt(RT*RT+RL*RL);
543 printf("ipl %2d-%2d dtt=%7.4lf dtl=%7.4lf rmst=%7.4lf rmsl=%7.4lf RT=%5.2lf RL=%5.2lf n=%d\n",
544 ipl1, ipl2, dtt, dtl, rt, rl, RT,RL,ndata);
545 if(Rmax<R){
546 Rmax = R;
547 RmaxT = RT;
548 RmaxL = RL;
549 iplRmax1=ipl1;
550 iplRmax2=ipl2;
551 dtmaxt = dtt;
552 dtmaxl = dtl;
553 rmstmax = rt;
554 rmslmax = rl;
555 }
556 }
557 }
558 // nseg in 5 plate.
559 int nin5pl = 0;
560 for(int k=0;k<n;k++){
561 EdbSegP *s = tt->GetSegment(k);
562 if(s->Plate()>=vtxpl && s->Plate()<vtxpl+5) nin5pl++;
563 }
564 if( nin5pl<=2 ){
565 // nseg <=2 in first 5 plate is out of criteria. just calculate RMS but not R.
566 RmaxT=RmaxL=-1.;
567 iplRmax1=iplRmax2=-1;
568 }
569
570 double rms, rmst, rmsl;
571 if(n<=1) rms=rmst=rmsl=-1.;
572 else DTRMSTL(tt, &rms, &rmst, &rmsl);
573
574 delete tt;
575
576 int result=1;
577
578 const char *comment = set->GetComment(t);
580 // Bern format
582 fprintf(fp,"// tid v1 v2 x y z tx ty ip1 ip2 ");
583 fprintf(fp," p pmin pmax mn pa sb dk of ps ");
584 fprintf(fp," n RmaxT RmaxL rmsT rmsL pl1 pl2 res\n");
585 }
586 fprintf(fp, "%6d %2d %2d %8.1f %8.1f %8.1f %7.4f %7.4f %6.1f %6.1f ",
587 t->ID(), ivtx1, ivtx2, sf->X(), sf->Y(), sf->Z(), sf->TX(), sf->TY(), ip1, ip2);
588
589 } else {
590 // Official
591 fprintf(fp, "%d %2d %2d %8.1f %8.1f %8.1f %7.4f %7.4f %6.1f %6.1f ",
592 itrk++, ivtx1, ivtx2, sf->X(), sf->Y(), sf->Z(), sf->TX(), sf->TY(), ip1, ip2);
593 }
594 fprintf(fp, "%5.1lf %5.1lf %5.1lf %2d %2d %2d %2d %2d %2d ",
595 p,pmin,pmax, MAN, MUON, SCANBACK, DARKNESS, OutOfBrick, LastPlate);
596 fprintf(fp, "%2d %5.2lf %5.2lf %7.4lf %7.4lf %3d %3d %3d\n",
597 nseg, RmaxT, RmaxL, rmst, rmsl, iplRmax1, iplRmax2, result);
598
599 // print on the shell.
600 fprintf(stdout, "%5d %2d %2d %8.1f %8.1f %8.1f %7.4f %7.4f %6.1f %6.1f ",
601 t->ID(), ivtx1, ivtx2, sf->X(), sf->Y(), sf->Z(), sf->TX(), sf->TY(), ip1, ip2);
602 fprintf(stdout, "%5.1lf %5.1lf %5.1lf %2d %2d %2d %2d %2d %2d ",
603 p,pmin,pmax, MAN, MUON, SCANBACK, DARKNESS, OutOfBrick, LastPlate);
604 fprintf(stdout, "%2d %5.2lf %5.2lf %7.4lf %7.4lf %3d %3d %3d\n",
605 nseg, RmaxT, RmaxL, rmst, rmsl, iplRmax1, iplRmax2, result);
606
607 // Check if the rms is smaller than angular resolution or not.
608 // if it's smaller, then print a note in the feedback file.
609
610 // minimum value of rms (angular resolution).
611 double angle = sqrt(t->TX()*t->TX()+t->TY()*t->TY());
612 double resT = 0.0015*sqrt(2.0);
613 double resL = 0.0015*(1+angle*5)*sqrt(2.0);
614 double RresT = RmaxT>0 ? dtmaxt/resT : -1.;
615 double RresL = RmaxL>0 ? dtmaxl/resL : -1.;
616
617 if( eOutputFileMode==kBern ){
618 if( ( rmst<resT || rmsl<resL ) && rmst>0)
619 fprintf(fp, "// rms for R = %7.4lf %7.4lf , note : rms is smaller than angular resolution. "
620 "RresT RresL resT resL = %5.2lf %5.2lf %7.4lf %7.4lf\n",
621 rmstmax, rmslmax, RresT, RresL, resT, resL);
622 else fprintf(fp, "// rms for R = %7.4lf %7.4lf "
623 "RresT RresL resT resL = %5.2lf %5.2lf %7.4lf %7.4lf\n",
624 rmstmax, rmslmax, RresT, RresL, resT, resL);
625
626 if( nin5pl<3 ) fprintf(fp, "// nBT = %d. nBT in first five plate = %d.\n", n, nin5pl);
627 fprintf(fp, "// track comment = \"%s\"\n", comment);
628 }
629
630 // Fill microtrack information
631
632 int irec = 0; // 0:TS, 1:PredScan, 2: Manual
633 TString name = set->GetName();
634 if( name.BeginsWith("TS") ) irec=0;
635 if( name.BeginsWith("SB") ) irec=1;
636 if( name.BeginsWith("SF") ) irec=1;
637 if( name.BeginsWith("MN") ) irec=2;
638
639 for(int k=0;k<t->N();k++){
640 EdbSegP *s = t->GetSegment(k);
641
642 int type = s->Side(); // 0:BT, 1:MTtop, 2:MTbottom
643 EdbID idseg=s->ScanID();
644 EdbEDATrackSet *setseg = gEDA->GetTrackSet( idseg);
645 if(setseg){
646 TString name = setseg->GetName();
647 if( name.BeginsWith("TS") ) irec=0;
648 if( name.BeginsWith("SB") ) irec=1;
649 if( name.BeginsWith("SF") ) irec=1;
650 if( name.BeginsWith("MN") ) irec=2;
651 }
652
654 // Bern format
655 if(k==0) fprintf(fp,"// ipl x y z tx ty type irec ngrains\n");
656 fprintf(fp,"%12d %8.1f %8.1f %8.1f %7.4f %7.4f %4d %4d %4d\n",
657 s->Plate(), s->X(), s->Y(), s->Z(), s->TX(), s->TY(), type, irec, (int) s->W());
658 } else {
659 // Official
660 fprintf(fp,"%d %8.1f %8.1f %8.1f %7.4f %7.4f %4d %4d %4d\n",
661 s->Plate(), s->X(), s->Y(), s->Z(), s->TX(), s->TY(), type, irec, (int) s->W());
662 }
663 }
664
665 }
666 }
667
668 // Scanning area
669
671 EdbEDAAreaSet *aset = set->GetAreaSet();
672 if(aset->N()!=0){
673 fprintf(fp,"//////////////////////////////////////////\n");
674 fprintf(fp,"// ------ Area information : %s ------- //\n", set->GetName());
675 fprintf(fp,"//////////////////////////////////////////\n");
676 for(int i=0;i<aset->N();i++){
677 EdbEDAArea *a = aset->GetArea(i);
678 fprintf(fp,"// area : TS %2d %lf %lf %lf %lf %lf\n",
679 a->Plate(), a->Xmin(), a->Xmax(),
680 a->Ymin(), a->Ymax(), a->Z());
681 }
682 }
683
684 fclose(fp);
685
686 printf("file %s was written.\n", filename);
687
688}
FILE * stdout
Expr< UnaryOp< Fabs< T >, Expr< A, T, D >, T >, T, D > fabs(const Expr< A, T, D > &rhs)
Definition: UnaryOperators.hh:96
int N()
Definition: EdbEDASets.h:133
EdbEDAArea * GetArea(int i)
Definition: EdbEDASets.h:134
double GetZ(int ipl)
Definition: EdbEDA.h:300
int GetIPLZ(float z)
Definition: EdbEDA.h:267
int GetLastPlate()
Definition: EdbEDA.h:241
int IdMuon(char *filename="../cs_info.txt")
Definition: EdbEDA.C:1045
EdbTrackP * FindTrack(EdbSegP *s, double dx=10, double dt=0.02)
Definition: EdbEDATrackSet.h:369
EdbTrackP * CheckScanback(EdbTrackP *t)
Definition: EdbEDA.h:730
EdbID * GetID(Int_t i)
Definition: EdbScanSet.h:46
Float_t TX() const
tangens = deltaX/deltaZ
Definition: EdbSegP.h:175
Int_t ID() const
Definition: EdbSegP.h:147
Float_t X() const
Definition: EdbSegP.h:173
Int_t Plate() const
Definition: EdbSegP.h:159
Float_t Z() const
Definition: EdbSegP.h:153
Float_t Y() const
Definition: EdbSegP.h:174
Float_t TY() const
tangens = deltaY/deltaZ
Definition: EdbSegP.h:176
void RemoveSegment(EdbSegP *s)
Definition: EdbPattern.h:219
Int_t N() const
Definition: EdbPattern.h:177
EdbSegP * GetSegment(int i) const
Definition: EdbPattern.h:195
void Copy(const EdbTrackP &tr)
Definition: EdbPattern.cxx:473
EdbTrackP * GetTrack(int i)
Definition: EdbVertex.h:141
const char * name
Definition: merge_Energy_SytematicSources_Electron.C:24
double DTRMSTLGiven1Kink(EdbTrackP *t, int iKink, double *rmsspace, double *rmstransverse, double *rmslongitudinal, int *NKinkAngleUsed=NULL)
Definition: EdbEDAUtil.C:790
void CalcDTTransLongi(EdbSegP *s1, EdbSegP *s2, double *dtTransverse, double *dtLongitudinal)
Definition: EdbEDAUtil.C:648
EdbMomentumEstimator * CalcP(EdbTrackP *t, double &p, double &pmin, double &pmax, bool print=kTRUE)
Definition: EdbEDAUtil.C:369
double DTRMSTL(EdbTrackP *t, double *rmsspace, double *rmstransverse, double *rmslongitudinal, int *ndata=NULL)
Definition: EdbEDAUtil.C:686
void r(int rid=2)
Definition: test.C:201

◆ WriteFeedbackFile2008()

void EdbEDAIO::WriteFeedbackFile2008 ( EdbVertex v = NULL,
char *  filename = NULL 
)

write feedback file. original code from Frederic.
the flags are correct if you work in the correct directory.

690 {
691
694
695 // extract vertex position and multiplicity
696 // extract vertex tracks infos
697 if(v==NULL) v = gEDA->GetSelectedVertex();
698 if(filename==NULL){ // File dialog
699 TGFileInfo *fi=new TGFileInfo;
700 fi->fIniDir = StrDup(".");
701 const char *filetypes[] = { "Feedback file", "*.feedback", "All files","*",0,0};
702 fi->fFileTypes = filetypes;
703 new TGFileDialog(gClient->GetRoot(), gEve->GetMainWindow(), kFDSave, fi);
704 if(fi->fFilename!=NULL) filename = fi->fFilename;
705 }
706
707 if(NULL==filename||strlen(filename)==0) { printf("Cancel to write feedback\n"); return;}
708 FILE *fp=fopen(filename,"wt");
709 if(NULL==fp){
710 printf("Couldn't open file : %s\n", filename);
711 return;
712 }
713
714
715 // put first the primary vertex
716
717 fprintf(fp,"%i %.1f %.1f %.1f %i %i %i \n",1,v->X(),v->Y(),v->Z(),1,0,0);
718
719 // tracks from primary vertex
720 for(int j=0; j<v->N(); j++){
721 int MANUAL=0, MUON=0, SCANBACK=0, DARKNESS=0;
722 EdbTrackP *t = v->GetTrack(j);
723 EdbSegP *s = t->GetSegmentFirst();
724
725 // check if the track is a scanback track
726 int idscanback = -1;
727 EdbTrackP *tsb = gEDA->CheckScanback(t);
728 if(tsb) idscanback = tsb->ID();
729
730 int idmuon = IdMuon();
731 if (idscanback!=-1) SCANBACK = 1;
732 if (idscanback==idmuon&&idmuon!=-1) {
733 MUON =1;
734 printf(" Track: %d is muon\n", t->ID());
735 }
736
737 double p, pmin, pmax;
738 CalcP(t, p, pmin, pmax); // pmin = pmin, pmax = pmax
739
740 fprintf(fp,"%4i %i %i %8.1f %8.1f %8.1f %7.4f %7.4f %6.2f %3.1f %6.2lf %5.2lf %8.1lf %i %i %i %i \n",
741 t->ID(),1,-1,s->X(),s->Y(),s->Z(),s->TX(), s->TY(), v->Impact(j),-1.,p,pmin,pmax,MANUAL,MUON,SCANBACK,DARKNESS);
742
743 }
744 printf("-- FeedbackFile : %s--\n", filename);
745 fclose(fp);
746
747#ifdef _WINDOWS
748 char buf[256];
749 sprintf(buf, "type %s",filename);
750 for(int i=0;i<strlen(buf);i++) if(buf[i]=='/') buf[i]='\\';
751 gSystem->Exec(buf);
752#else
753 gSystem->Exec(Form("cat %s",filename));
754#endif
755}
EdbVertex * GetSelectedVertex(void)
Definition: EdbEDA.h:439
Float_t Impact(int i)
Definition: EdbVertex.h:149

◆ WriteFilteredText()

void EdbEDAIO::WriteFilteredText ( char *  filename = NULL)
1375 {
1376 int i,j,k,l;
1377 FILE *fp;
1378
1379 if(filename==NULL){ // File dialog
1380 TGFileInfo *fi=new TGFileInfo;
1381 fi->fIniDir = StrDup(".");
1382 const char *filetypes[] = { "filtered tracks", "*.txt", "All files","*",0,0};
1383 fi->fFileTypes = filetypes;
1384 new TGFileDialog(gClient->GetRoot(), gEve->GetMainWindow(), kFDSave, fi);
1385 if(fi->fFilename!=NULL) {
1386 filename = fi->fFilename;
1387 }
1388 }
1389 if(filename==NULL) return;
1390
1391 fp=fopen(filename,"wt");
1392
1393 if(fp==NULL) {
1394 printf("File open error : %s\n", filename);
1395 return;
1396 }
1397
1398 printf("Write Filtered Text File : %s\n", filename);
1399
1400 // Vertices
1401
1402 TObjArray *vertices = (gEDA->GetVertexSet())->GetVertices();
1403
1404 EdbVertex *v_1ry=NULL;
1405 double z=1e9;
1406 for(i=0;i<vertices->GetEntries();i++){
1407 EdbVertex *v = (EdbVertex *) vertices->At(i);
1408 if(v->Z()<z) {
1409 z = v->Z();
1410 v_1ry=v;
1411 }
1412 }
1413
1414 fprintf(fp,"// vtx id X Y Z 1ry charm tau\n");
1415 for( i=0; i<vertices->GetEntries(); i++){
1416 EdbVertex *v = (EdbVertex *) vertices->At(i);
1417 v->SetID(i);
1418 int ISPRIMARY = v==v_1ry ? 1 : 0;
1419 int ISCHARM = 0;
1420 int ISTAU = 0;
1421 fprintf(fp,"vertex %3d %8.1lf %8.1lf %8.1lf %3d %5d %3d\n",
1422 v->ID(), v->X(), v->Y(), v->Z(),ISPRIMARY,ISCHARM,ISTAU);
1423 }
1424
1425 printf("vertices %d\n", vertices->GetEntries());
1426
1427// ID_Track ID_UpVtx ID_DownVtx X Y Z SX SY IP_UpVtx IP_DownVtx P Pmin Pmax Manual Mu/e Scanback Darknes
1428
1429 fprintf(fp,"// trk id flag tx ty pl1 pl2 nseg ip1ry vtx1 ip1 vtx2 ip2 p pmin pmax man part sb dark\n");
1430
1431 // Tracks + Segments
1432
1433 for(i=0;i<gEDA->NTrackSets();i++){
1435 if(set->GetDraw()==0) continue;
1436 for(j=0;j<set->N();j++){
1437 EdbTrackP *t = set->GetTrack(j);
1438 EdbSegP *sf = t->GetSegmentFirst();
1439 EdbSegP *sl = t->GetSegmentLast();
1440 EdbVertex * vtx1st = NULL;
1441 EdbVertex * vtx2nd = NULL;
1442 for(k=0;k<vertices->GetEntries();k++){
1443 EdbVertex *v = (EdbVertex *) vertices->At(k);
1444 for(l=0;l<v->N();l++){
1445 EdbTrackP *tv = v->GetTrack(l);
1446 if(t==tv){
1447 if(NULL==vtx1st) vtx1st = v;
1448 else {
1449 if(vtx1st->Z()>v->Z()) { // vtx1st should 1ry.
1450 vtx2nd = vtx1st;
1451 vtx1st = v;
1452 }
1453 else vtx2nd = v;
1454 } } } }
1455
1456 // Flags.... same as Feedback file.
1457 int MAN = 0;
1458 int SCANBACK=0, idscanback=-1;
1459 int MUON =0;
1460 int DARKNESS=0; // no function to define.
1461
1462 if(-99<=t->Flag()&&t->Flag()<0) MAN=1;
1463
1464 EdbTrackP * sbt = gEDA->CheckScanback(t); // check if the track is a scanback track
1465 if(sbt!=NULL) idscanback = sbt->ID();
1466 if (idscanback!=-1) SCANBACK = 1;
1467 if (idscanback==IdMuon()&&IdMuon()!=-1) MUON =1;
1468
1469 double p,pmin,pmax;
1470 CalcP(t, p, pmin, pmax); // pmin = pmin, pmax = pmax
1471 int ivtx1 = vtx1st ? vtx1st->ID() : -1;
1472 int ivtx2 = vtx2nd ? vtx2nd->ID() : -1;
1473 float ip1=-1.0, ip2=-1.0;
1474 if(vtx1st){
1475 for(k=0;k<vtx1st->N();k++){
1476 EdbVTA *vta = vtx1st->GetVTa(k);
1477 if(t==vta->GetTrack()) ip1 = vta->Imp();
1478 }
1479 }
1480 if(vtx2nd){
1481 for(k=0;k<vtx2nd->N();k++){
1482 EdbVTA *vta = vtx2nd->GetVTa(k);
1483 if(t==vta->GetTrack()) ip2 = vta->Imp();
1484 }
1485 }
1486 double ip_1ry = CalcIP(sf, v_1ry);
1487
1488 fprintf(fp,"track %6d %4d %7.4f %7.4f %3d %3d %4d %7.1lf %4d %5.1f %4d %5.1f ",
1489 t->ID(), t->Flag(), t->TX(), t->TY(), sf->Plate(), sl->Plate(), t->N(),
1490 ip_1ry, ivtx1, ip1, ivtx2, ip2);
1491 fprintf(fp,"%6.2lf %6.2lf %10.2lf %3d %4d %2d %4d\n", p,pmin,pmax, MAN, MUON, SCANBACK, DARKNESS);
1492 }
1493 printf("TrackSet : %s %d tracks\n", set->GetName(), set->N());
1494 }
1495
1496 fprintf(fp,"// seg trkid tflag ipl segid x y z tx ty w chi2 flag\n");
1497
1498 for(i=0;i<gEDA->NTrackSets();i++){
1500 if(set->GetDraw()==0) continue;
1501 for(j=0;j<set->N();j++){
1502 EdbTrackP *t = set->GetTrack(j);
1503 for(k=0;k<t->N();k++){
1504 EdbSegP *s = t->GetSegment(k);
1505 fprintf(fp,"segment %6d %5d %3d %6d %8.1f %8.1f %8.1f %7.4f %7.4f %2d %6.4f %4d\n",
1506 t->ID(), t->Flag(), s->Plate(), s->ID(), s->X(),s->Y(),s->Z(),s->TX(),s->TY(),
1507 (int)s->W(),s->Chi2(),s->Flag());
1508 }
1509 }
1510 }
1511
1512
1513 // Scanning area
1514 fprintf(fp,"// area ipl xmin xmax ymin ymax z\n");
1515 EdbEDAAreaSet *aset = gEDA->GetTrackSet("TS")->GetAreaSet();
1516 for(i=0; i<aset->N(); i++){
1517 EdbEDAArea *a = aset->GetArea(i);
1518 fprintf(fp,"area %2d %8.1lf %8.1lf %8.1lf %8.1lf %8.1lf\n",
1519 a->Plate(), a->Xmin(), a->Xmax(), a->Ymin(), a->Ymax(), a->Z());
1520 }
1521
1522
1523
1524 fclose(fp);
1525}
Float_t Imp() const
Definition: EdbVertex.h:49

Member Data Documentation

◆ eCutCP

TCut EdbEDAIO::eCutCP

◆ eOutputFileMode

int EdbEDAIO::eOutputFileMode

The documentation for this class was generated from the following files: