4 |
ClassImp(CaloLong); |
ClassImp(CaloLong); |
5 |
ClassImp(Calo2D); |
ClassImp(Calo2D); |
6 |
|
|
7 |
|
|
8 |
|
int isempty(struct stack *s){ |
9 |
|
return (s->top == EMPTY) ? 1 : 0; |
10 |
|
} |
11 |
|
|
12 |
|
void emptystack(struct stack* s){ |
13 |
|
s->top=EMPTY; |
14 |
|
strcpy(s->data," "); |
15 |
|
} |
16 |
|
|
17 |
|
void push(struct stack* s,int item){ |
18 |
|
if(s->top == (MAX-1)){ |
19 |
|
printf("\n ERROR! STACK FULL (too many digits in SetLower/UpperLimit)\n"); |
20 |
|
} else { |
21 |
|
s->data[++s->top]=(char)item; |
22 |
|
}; |
23 |
|
} |
24 |
|
|
25 |
|
char pop(struct stack* s){ |
26 |
|
char ret=(char)EMPTY; |
27 |
|
if(!isempty(s)){ |
28 |
|
ret= s->data[s->top--]; |
29 |
|
}; |
30 |
|
return ret; |
31 |
|
} |
32 |
|
|
33 |
|
|
34 |
|
int ipop(struct stack* s){ |
35 |
|
int ret=EMPTY; |
36 |
|
if(s->top == EMPTY) |
37 |
|
printf("\n ERROR! STACK EMPTY (too few digits in SetLower/UpperLimit)\n"); |
38 |
|
else { |
39 |
|
ret= s->data[s->top--]; |
40 |
|
}; |
41 |
|
return ret; |
42 |
|
} |
43 |
|
|
44 |
|
|
45 |
|
void display(struct stack s){ |
46 |
|
int ss = s.top; |
47 |
|
while(s.top != EMPTY){ |
48 |
|
printf("\n%d",s.data[s.top--]); |
49 |
|
}; |
50 |
|
printf(" s.top %i \n",ss); |
51 |
|
} |
52 |
|
|
53 |
|
int isoperator(char e){ |
54 |
|
if(e == '+' || e == '-' || e == '*' || e == '/' || e == '%') |
55 |
|
return 1; |
56 |
|
else |
57 |
|
return 0; |
58 |
|
} |
59 |
|
|
60 |
|
|
61 |
|
int priority(char e){ |
62 |
|
int pri = 0; |
63 |
|
if(e == '*' || e == '/' || e =='%') |
64 |
|
pri = 2; |
65 |
|
else { |
66 |
|
if(e == '+' || e == '-') pri = 1; |
67 |
|
}; |
68 |
|
return pri; |
69 |
|
} |
70 |
|
|
71 |
|
void infix2postfix(char* infix, char * postfix, int insertspace){ |
72 |
|
char *i,*p; |
73 |
|
struct stack X; |
74 |
|
char n1; |
75 |
|
emptystack(&X); |
76 |
|
i = &infix[0]; |
77 |
|
p = &postfix[0]; |
78 |
|
while(*i){ |
79 |
|
while(*i == ' ' || *i == '\t' ){ |
80 |
|
i++; |
81 |
|
}; |
82 |
|
TString c=i; |
83 |
|
if( isdigit(*i) || isalpha(*i) || c.BeginsWith(".") ){ |
84 |
|
while( isdigit(*i) || isalpha(*i) || c.BeginsWith(".") ){ |
85 |
|
*p = *i; |
86 |
|
p++; |
87 |
|
i++; |
88 |
|
c=i; |
89 |
|
}; |
90 |
|
/*SPACE CODE*/ |
91 |
|
if(insertspace){ |
92 |
|
*p = ' '; |
93 |
|
p++; |
94 |
|
}; |
95 |
|
/*END SPACE CODE*/ |
96 |
|
}; |
97 |
|
if( *i == '(' ){ |
98 |
|
push(&X,*i); |
99 |
|
i++; |
100 |
|
}; |
101 |
|
if( *i == ')'){ |
102 |
|
n1 = pop(&X); |
103 |
|
while( n1 != '(' ){ |
104 |
|
*p = n1; |
105 |
|
p++; |
106 |
|
/*SPACE CODE*/ |
107 |
|
if(insertspace){ |
108 |
|
*p = ' '; |
109 |
|
p++; |
110 |
|
}; |
111 |
|
/*END SPACE CODE*/ |
112 |
|
n1 = pop(&X); |
113 |
|
}; |
114 |
|
i++; |
115 |
|
}; |
116 |
|
if( isoperator(*i) ){ |
117 |
|
if(isempty(&X)){ |
118 |
|
push(&X,*i); |
119 |
|
} else { |
120 |
|
n1 = pop(&X); |
121 |
|
while( priority(n1) >= priority(*i) ){ |
122 |
|
*p = n1; |
123 |
|
p++; |
124 |
|
/*SPACE CODE*/ |
125 |
|
if(insertspace){ |
126 |
|
*p = ' '; |
127 |
|
p++; |
128 |
|
}; |
129 |
|
/*END SPACE CODE*/ |
130 |
|
n1 = pop(&X); |
131 |
|
}; |
132 |
|
push(&X,n1); |
133 |
|
push(&X,*i); |
134 |
|
}; |
135 |
|
i++; |
136 |
|
}; |
137 |
|
}; |
138 |
|
while(!isempty(&X)){ |
139 |
|
n1 = pop(&X); |
140 |
|
*p = n1; |
141 |
|
p++; |
142 |
|
/*SPACE CODE*/ |
143 |
|
if(insertspace){ |
144 |
|
*p = ' '; |
145 |
|
p++; |
146 |
|
}; |
147 |
|
/*END SPACE CODE*/ |
148 |
|
}; |
149 |
|
*p = '\0'; |
150 |
|
} |
151 |
|
|
152 |
|
|
153 |
|
Float_t evaluate(char *postfix){ |
154 |
|
// |
155 |
|
Float_t op1 = 0.; |
156 |
|
Float_t op2 = 0.; |
157 |
|
Float_t result = 0.; |
158 |
|
// |
159 |
|
TString e = postfix; |
160 |
|
Float_t st[50]; |
161 |
|
memset(st,0,50*sizeof(Float_t)); |
162 |
|
TObjArray *ae = e.Tokenize(" "); |
163 |
|
Int_t o = 0; |
164 |
|
Int_t a = 0; |
165 |
|
Int_t ap = 0; |
166 |
|
// |
167 |
|
while ( o < ae->GetEntries() ){ |
168 |
|
if ( ((TString)(((TObjString*)ae->At(o))->GetString())).IsFloat() ){ |
169 |
|
st[a] = (Float_t)((TString)(((TObjString*)ae->At(o))->GetString())).Atof();; |
170 |
|
a++; |
171 |
|
} else { |
172 |
|
ap = a-1; |
173 |
|
op1 = st[ap--]; |
174 |
|
op2 = st[ap]; |
175 |
|
const char *p=((TString)(((TObjString*)ae->At(o))->GetString())).Data(); |
176 |
|
switch(*p){ |
177 |
|
case '+': |
178 |
|
result = op2 + op1; |
179 |
|
break; |
180 |
|
case '-': |
181 |
|
result = op2 - op1; |
182 |
|
break; |
183 |
|
case '/': |
184 |
|
result = op2 / op1; |
185 |
|
break; |
186 |
|
case '*': |
187 |
|
result = op2 * op1; |
188 |
|
break; |
189 |
|
case '%': |
190 |
|
result = (Int_t)round(op2) % (Int_t)round(op1); |
191 |
|
break; |
192 |
|
default: |
193 |
|
printf("\nInvalid Operator: %s \n",((TString)(((TObjString*)ae->At(o))->GetString())).Data()); |
194 |
|
return 0; |
195 |
|
}; |
196 |
|
st[ap] = result; |
197 |
|
a = ap+1; |
198 |
|
// |
199 |
|
}; |
200 |
|
o++; |
201 |
|
}; |
202 |
|
return result; |
203 |
|
// |
204 |
|
} |
205 |
|
|
206 |
|
|
207 |
|
|
208 |
//////////////////////////////////////////////////////////////////////// |
//////////////////////////////////////////////////////////////////////// |
209 |
/** |
/** |
210 |
* 1-dimension function describing lateral distribution of the |
* 1-dimension function describing lateral distribution of the |
652 |
PKT = 0; |
PKT = 0; |
653 |
atime = 0; |
atime = 0; |
654 |
// |
// |
655 |
|
clp = L2->GetCaloLevel2(); |
656 |
|
// |
657 |
sel = true; |
sel = true; |
658 |
cont = false; |
cont = false; |
659 |
N = 0; |
N = 0; |
667 |
maskYE = false; |
maskYE = false; |
668 |
maskYO = false; |
maskYO = false; |
669 |
// |
// |
670 |
|
lmax = 0.; |
671 |
|
umax = 100.; |
672 |
|
slmax = ""; |
673 |
|
sumax = ""; |
674 |
|
xyaverage = true; |
675 |
|
// |
676 |
}; |
}; |
677 |
|
|
678 |
void CaloLong::MaskSection(TString sec){ |
void CaloLong::MaskSection(TString sec){ |
683 |
if ( sec.Contains("YE") ) maskYE = true; |
if ( sec.Contains("YE") ) maskYE = true; |
684 |
} |
} |
685 |
|
|
686 |
|
void CaloLong::UnMaskSections(){ |
687 |
|
this->UnMaskSection("XEXOYEYO"); |
688 |
|
} |
689 |
|
|
690 |
|
void CaloLong::UnMaskSection(TString sec){ |
691 |
|
sec.ToUpper(); |
692 |
|
if ( sec.Contains("XO") ) maskXO = false; |
693 |
|
if ( sec.Contains("YO") ) maskYO = false; |
694 |
|
if ( sec.Contains("XE") ) maskXE = false; |
695 |
|
if ( sec.Contains("YE") ) maskYE = false; |
696 |
|
} |
697 |
|
|
698 |
void CaloLong::Clear(){ |
void CaloLong::Clear(){ |
699 |
// |
// |
700 |
memset(eplane,0, 2*22*sizeof(Float_t)); |
memset(eplane,0, 2*22*sizeof(Float_t)); |
702 |
chi2 = 0.; |
chi2 = 0.; |
703 |
ndf = 0.; |
ndf = 0.; |
704 |
E0 = 0.; |
E0 = 0.; |
705 |
|
defE0 = 0.; |
706 |
a = 0.; |
a = 0.; |
707 |
b = 0.; |
b = 0.; |
708 |
errE0 = 0.; |
errE0 = 0.; |
728 |
printf(" nchi2 :.. %f\n",chi2/ndf); |
printf(" nchi2 :.. %f\n",chi2/ndf); |
729 |
printf(" E0 :.. %f\n",E0); |
printf(" E0 :.. %f\n",E0); |
730 |
printf(" E0/260. :.. %f\n",E0/260.); |
printf(" E0/260. :.. %f\n",E0/260.); |
731 |
|
printf(" defE0 :.. %f\n",defE0); |
732 |
|
printf(" lower :.. %f\n",lmax); |
733 |
|
printf(" upper :.. %f\n",umax); |
734 |
|
printf(" s lower :.. %s\n",slmax.Data()); |
735 |
|
printf(" s upper :.. %s\n",sumax.Data()); |
736 |
printf(" a :.. %f\n",a); |
printf(" a :.. %f\n",a); |
737 |
printf(" b :.. %f\n",b); |
printf(" b :.. %f\n",b); |
738 |
printf(" errE0 :.. %f\n",errE0); |
printf(" errE0 :.. %f\n",errE0); |
851 |
// |
// |
852 |
Float_t ytgx = 0; |
Float_t ytgx = 0; |
853 |
Float_t ytgy = 0; |
Float_t ytgy = 0; |
854 |
ytgx = 0.76 * L2->GetCaloLevel2()->tanx[0]; |
ytgx = 0.76 * clp->tanx[0]; |
855 |
ytgy = 0.76 * L2->GetCaloLevel2()->tany[0]; |
ytgy = 0.76 * clp->tany[0]; |
856 |
X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) ); |
X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) ); |
857 |
// |
// |
858 |
// Find experimental plane of maximum |
// Find experimental plane of maximum |
880 |
// |
// |
881 |
}; |
}; |
882 |
|
|
883 |
|
void CaloLong::SetEnergies(Float_t myene[][22]){ |
884 |
|
// |
885 |
|
if ( !L2 ){ |
886 |
|
printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n"); |
887 |
|
printf(" ERROR: CaloHough variables not filled \n"); |
888 |
|
return; |
889 |
|
}; |
890 |
|
// |
891 |
|
Bool_t newentry = false; |
892 |
|
// |
893 |
|
if ( L2->IsORB() ){ |
894 |
|
if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){ |
895 |
|
newentry = true; |
896 |
|
OBT = L2->GetOrbitalInfo()->OBT; |
897 |
|
PKT = L2->GetOrbitalInfo()->pkt_num; |
898 |
|
atime = L2->GetOrbitalInfo()->absTime; |
899 |
|
}; |
900 |
|
} else { |
901 |
|
newentry = true; |
902 |
|
}; |
903 |
|
// |
904 |
|
if ( !newentry ) return; |
905 |
|
// |
906 |
|
if ( debug ) printf(" SET ENERGIES Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime); |
907 |
|
// |
908 |
|
Clear(); |
909 |
|
// |
910 |
|
// let's start |
911 |
|
// |
912 |
|
if ( cont ){ |
913 |
|
for (Int_t i=0; i<22; i++){ |
914 |
|
if ( i == (18+N) ){ |
915 |
|
mask18b = 18 + N; |
916 |
|
break; |
917 |
|
}; |
918 |
|
}; |
919 |
|
}; |
920 |
|
// |
921 |
|
if ( sel ){ |
922 |
|
for (Int_t i=0; i<22; i++){ |
923 |
|
if ( i == (18-N) ){ |
924 |
|
mask18b = 18 - N; |
925 |
|
break; |
926 |
|
}; |
927 |
|
}; |
928 |
|
}; |
929 |
|
// |
930 |
|
// if ( mask18b == 18 ) mask18b = -1; |
931 |
|
// |
932 |
|
Int_t view = 0; |
933 |
|
Int_t plane = 0; |
934 |
|
Bool_t gof = true; |
935 |
|
for (view=0; view < 2; view++){ |
936 |
|
for (plane=0; plane < 22; plane++){ |
937 |
|
gof = true; |
938 |
|
if ( maskXE && (plane%2)==0 && view==1 ) gof = false; |
939 |
|
if ( maskXO && (plane%2)!=0 && view==1 ) gof = false; |
940 |
|
if ( maskYE && (plane%2)!=0 && view==0 ) gof = false; |
941 |
|
if ( maskYO && (plane%2)==0 && view==0 ) gof = false; |
942 |
|
if ( gof ) eplane[view][plane] = myene[view][plane]; |
943 |
|
if ( debug ) printf(" XE %i XO %i YE %i YO %i eplane %i %i = %f ........... myene %i %i = %f\n", maskXE, maskXO, maskYE, maskYO,view,plane,eplane[view][plane],view,plane,myene[view][plane]); |
944 |
|
}; |
945 |
|
}; |
946 |
|
// |
947 |
|
// inclination factor (stolen from Daniele's code) |
948 |
|
// |
949 |
|
Float_t ytgx = 0; |
950 |
|
Float_t ytgy = 0; |
951 |
|
ytgx = 0.76 * clp->tanx[0]; |
952 |
|
ytgy = 0.76 * clp->tany[0]; |
953 |
|
X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) ); |
954 |
|
// |
955 |
|
// Find experimental plane of maximum |
956 |
|
// |
957 |
|
Int_t pmax = 0; |
958 |
|
Int_t vmax = 0; |
959 |
|
Float_t emax = 0.; |
960 |
|
for (Int_t v=0; v<2; v++){ |
961 |
|
for (Int_t i=0; i<22; i++){ |
962 |
|
if ( eplane[v][i] > emax ){ |
963 |
|
emax = eplane[v][i]; |
964 |
|
vmax = v; |
965 |
|
pmax = i; |
966 |
|
}; |
967 |
|
}; |
968 |
|
}; |
969 |
|
// |
970 |
|
// |
971 |
|
// |
972 |
|
if ( vmax == 0 ) pmax++; |
973 |
|
etmax = pmax * X0pl; |
974 |
|
// |
975 |
|
if ( debug ) this->Print(); |
976 |
|
if ( debug ) printf(" exit \n"); |
977 |
|
// |
978 |
|
}; |
979 |
|
|
980 |
Double_t ccurve(Double_t *ti,Double_t *par){ |
Double_t ccurve(Double_t *ti,Double_t *par){ |
981 |
// |
// |
995 |
this->Fit(false); |
this->Fit(false); |
996 |
}; |
}; |
997 |
|
|
998 |
|
Float_t CaloLong::Evaluate(TString s, Float_t tmax, Float_t X0pl){ |
999 |
|
/* SAMPLE OUTPUT: |
1000 |
|
Enter Infix Expression : A + B + C / (E - F) |
1001 |
|
Postfix Expression is : A B + C E F - / + |
1002 |
|
*/ |
1003 |
|
if ( !s.Contains("tmax") && !s.Contains("X0pl") ){ |
1004 |
|
printf(" ERROR, the input formula must contain \"t\"\n"); |
1005 |
|
return 0.; |
1006 |
|
}; |
1007 |
|
if ( tmax != tmax ){ |
1008 |
|
printf(" ERROR, tmax is nan! \n"); |
1009 |
|
return 0.; |
1010 |
|
}; |
1011 |
|
TString g=Form("%f",tmax); |
1012 |
|
TString h=Form("%f",X0pl); |
1013 |
|
TString *ts= new TString(""); |
1014 |
|
ts->Prepend(s.Data()); |
1015 |
|
ts->ReplaceAll("tmax",4,g.Data(),g.Capacity()); |
1016 |
|
ts->ReplaceAll("X0pl",4,h.Data(),h.Capacity()); |
1017 |
|
ts->Prepend("("); |
1018 |
|
ts->Append(")"); |
1019 |
|
if ( debug ) printf(" ts %s tssize %i capac %i s %s g %s \n",ts->Data(),ts->Sizeof(),ts->Capacity(),s.Data(),g.Data()); |
1020 |
|
char in[50],post[50]; |
1021 |
|
strcpy(&in[0]," "); |
1022 |
|
strcpy(&in[0],ts->Data()); |
1023 |
|
strcpy(&post[0]," "); |
1024 |
|
if ( debug ) printf("Infix Expression is : ##%s##\n",&in[0]); |
1025 |
|
infix2postfix(&in[0],&post[0],1); |
1026 |
|
if ( debug ) printf("Postfix Expression is : ##%s##\n",&post[0]); |
1027 |
|
/* SAMPLE OUTPUT: |
1028 |
|
Enter Postfix Expression : 3 5 + 2 / |
1029 |
|
3 5 + 2 / EQUALS 4 |
1030 |
|
*/ |
1031 |
|
Float_t res = evaluate(&post[0]); |
1032 |
|
if ( debug ) printf("%s EQUALS %f\n",post,res); |
1033 |
|
return res; |
1034 |
|
} |
1035 |
|
|
1036 |
void CaloLong::Fit(Bool_t draw){ |
void CaloLong::Fit(Bool_t draw){ |
1037 |
// |
// |
1038 |
Process(); |
Process(); |
1077 |
}; |
}; |
1078 |
// |
// |
1079 |
TString thid = Form("hlongfit"); |
TString thid = Form("hlongfit"); |
1080 |
TH1F *th = dynamic_cast<TH1F*>(gDirectory->FindObject(thid)); |
TH2F *th = dynamic_cast<TH2F*>(gDirectory->FindObject(thid)); |
1081 |
if ( th ) th->Delete(); |
if ( th ) th->Delete(); |
1082 |
|
// |
1083 |
|
TString ghid = Form("glongfit"); |
1084 |
|
TGraphErrors *gh = dynamic_cast<TGraphErrors*>(gDirectory->FindObject(ghid)); |
1085 |
|
if ( gh ) gh->Delete(); |
1086 |
|
// |
1087 |
Float_t xpos = 0.; |
Float_t xpos = 0.; |
1088 |
Float_t enemip = 0.; |
Float_t enemip = 0.; |
1089 |
Float_t xmax = NC * X0pl + 0.2; |
Float_t xmax = NC * X0pl + 0.2; |
1090 |
// th = new TH1F(thid,thid,int(NC*1.5),-0.2,xmax); |
// |
1091 |
th = new TH1F(thid,thid,100,-0.2,xmax); |
Double_t xxx[50]; |
1092 |
|
Double_t exx[50]; |
1093 |
|
Double_t yyy[50]; |
1094 |
|
Double_t eyy[50]; |
1095 |
|
Int_t numpo = 0; |
1096 |
|
memset(xxx,0,50*sizeof(Double_t)); |
1097 |
|
memset(exx,0,50*sizeof(Double_t)); |
1098 |
|
memset(yyy,0,50*sizeof(Double_t)); |
1099 |
|
memset(eyy,0,50*sizeof(Double_t)); |
1100 |
// |
// |
1101 |
// AGH, BUG! |
// AGH, BUG! |
1102 |
// |
// |
1110 |
mmax = NC; |
mmax = NC; |
1111 |
}; |
}; |
1112 |
// |
// |
1113 |
|
Float_t emax = 0.; |
1114 |
Float_t qtotparz = 0.; |
Float_t qtotparz = 0.; |
1115 |
for (Int_t st=mmin;st<mmax+1;st++){ |
for (Int_t st=mmin;st<mmax+1;st++){ |
1116 |
enemip = 0.; |
enemip = 0.; |
1117 |
xpos = (st - mmin) * X0pl; |
xpos = (st - mmin) * X0pl; |
1118 |
if ( st > mmin && st < mmax ){ |
// |
1119 |
if ( no18x && ( st == 18+1 || st == mask18b+1 )){ |
if ( xyaverage ){ |
1120 |
enemip = 2. * eplane[1][st]; |
// |
1121 |
|
if ( st > mmin && st < mmax ){ |
1122 |
|
if ( no18x && ( st == 18+1 || st == mask18b+1 )){ |
1123 |
|
if ( !maskYO ){ |
1124 |
|
enemip = 2. * eplane[1][st]; |
1125 |
|
} else { |
1126 |
|
enemip = eplane[1][st]; |
1127 |
|
}; |
1128 |
|
} else { |
1129 |
|
enemip = eplane[0][st-1] + eplane[1][st]; |
1130 |
|
}; |
1131 |
} else { |
} else { |
1132 |
enemip = eplane[0][st-1] + eplane[1][st]; |
if ( st == mmin ){ |
1133 |
|
if ( !maskYE ){ |
1134 |
|
enemip = 2. * eplane[1][st]; |
1135 |
|
} else { |
1136 |
|
enemip = eplane[1][st]; |
1137 |
|
}; |
1138 |
|
}; |
1139 |
|
if ( st == mmax ){ |
1140 |
|
if ( !maskXE ){ |
1141 |
|
enemip = 2. * eplane[0][st-1]; |
1142 |
|
} else { |
1143 |
|
enemip = eplane[0][st-1]; |
1144 |
|
}; |
1145 |
|
}; |
1146 |
|
}; |
1147 |
|
// |
1148 |
|
qtotparz += enemip; |
1149 |
|
if ( enemip > emax ) emax = enemip; |
1150 |
|
if ( enemip > 0. ){ |
1151 |
|
xxx[numpo] = xpos; |
1152 |
|
exx[numpo] = 0.1; |
1153 |
|
yyy[numpo] = enemip; |
1154 |
|
eyy[numpo] = sqrt(enemip*3.)+sqrt(5.); |
1155 |
|
numpo++; |
1156 |
|
// th->Fill(xpos,enemip); |
1157 |
|
if ( debug ) printf(" AVE Filling: st %i xpos %f energy %f \n",st,xpos,enemip); |
1158 |
}; |
}; |
1159 |
} else { |
} else { |
1160 |
if ( st == mmin ) enemip = 2. * eplane[1][st]; |
for (Int_t jj=0; jj<2; jj++){ |
1161 |
if ( st == mmax ) enemip = 2. * eplane[0][st-1]; |
if ( st > mmin && st < mmax ){ |
1162 |
}; |
if ( jj == 0 && no18x && ( st == 18+1 || st == mask18b+1 )){ |
1163 |
// |
enemip = eplane[1][st]; |
1164 |
qtotparz += enemip; |
} else { |
1165 |
if ( enemip > 0. ){ |
if ( jj == 0 ){ |
1166 |
th->Fill(xpos,enemip); |
enemip = eplane[jj][st-1]; |
1167 |
if ( debug ) printf(" Filling: st %i xpos %f energy %f \n",st,xpos,enemip); |
} else { |
1168 |
|
enemip = eplane[jj][st]; |
1169 |
|
}; |
1170 |
|
}; |
1171 |
|
} else { |
1172 |
|
if ( st == mmin && jj == 1 ){ |
1173 |
|
enemip = eplane[jj][st]; |
1174 |
|
}; |
1175 |
|
if ( st == mmax && jj == 0){ |
1176 |
|
enemip = eplane[jj][st-1]; |
1177 |
|
}; |
1178 |
|
}; |
1179 |
|
// |
1180 |
|
qtotparz += enemip; |
1181 |
|
if ( enemip > emax ) emax = enemip; |
1182 |
|
if ( enemip > 0. ){ |
1183 |
|
xxx[numpo] = xpos; |
1184 |
|
exx[numpo] = 0.1; |
1185 |
|
yyy[numpo] = enemip; |
1186 |
|
eyy[numpo] = sqrt(enemip*3.)+sqrt(5.); |
1187 |
|
// eyy[numpo] = sqrt(enemip)/(st*0.95); |
1188 |
|
numpo++; |
1189 |
|
// th->Fill(xpos,enemip); |
1190 |
|
if ( debug ) printf(" XY Filling: st %i xpos %f energy %f \n",st,xpos,enemip); |
1191 |
|
}; |
1192 |
|
}; |
1193 |
}; |
}; |
1194 |
|
|
1195 |
// |
// |
1196 |
// for (Int_t v=1; v>=0;v--)// { |
// for (Int_t v=1; v>=0;v--)// { |
1197 |
// // |
// // |
1217 |
// }; |
// }; |
1218 |
}; |
}; |
1219 |
// |
// |
1220 |
TF1 *lfit = new TF1("lfit",ccurve,0.,xmax,3); |
// th = new TH2F(thid,thid,int(NC*1.5),-0.2,xmax); |
1221 |
if ( debug ) printf("qtot %f qtotparz %f \n",L2->GetCaloLevel2()->qtot,qtotparz); |
th = new TH2F(thid,thid,1000,-0.2,xmax,1000,0.,emax*1.2); |
1222 |
|
gh = new TGraphErrors(numpo,xxx,yyy,exx,eyy); |
1223 |
|
TF1 *lfit = dynamic_cast<TF1*>(gDirectory->FindObject("lfit")); |
1224 |
|
if ( lfit ) lfit->Delete(); |
1225 |
|
lfit = new TF1("lfit",ccurve,0.,xmax,3); |
1226 |
|
if ( debug ) printf("qtot %f qtotparz %f \n",clp->qtot,qtotparz); |
1227 |
E0 = qtotparz; |
E0 = qtotparz; |
1228 |
// E0 = L2->GetCaloLevel2()->qtot; |
// E0 = clp->qtot; |
1229 |
a = 5.; |
a = 5.; |
1230 |
b = 0.5; |
b = 0.5; |
1231 |
if ( debug ) printf(" STARTING PARAMETERS: E0 %f a %f b %f \n",E0,a,b); |
if ( debug ) printf(" STARTING PARAMETERS: E0 %f a %f b %f \n",E0,a,b); |
1233 |
// lfit->SetParLimits(0,0.,1000.); |
// lfit->SetParLimits(0,0.,1000.); |
1234 |
// lfit->SetParLimits(1,-1.,80.); |
// lfit->SetParLimits(1,-1.,80.); |
1235 |
// lfit->SetParLimits(2,-1.,10.); |
// lfit->SetParLimits(2,-1.,10.); |
1236 |
TString optio; |
// TString optio = "ROW"; // "RO" |
1237 |
|
TString optio = "RO"; // "RO" |
1238 |
if ( debug ){ |
if ( debug ){ |
1239 |
// optio = "MERBOV"; |
optio += "NV"; |
1240 |
// optio = "MEROV"; |
if ( draw ) optio += "V"; |
1241 |
// optio = "EROV"; |
} else { |
1242 |
optio = "RNOV"; |
optio += "NQ"; |
1243 |
if ( draw ) optio = "ROV"; |
if ( draw ) optio += "Q"; |
|
} else { |
|
|
// optio = "MERNOQ"; |
|
|
// optio = "ERNOQ"; |
|
|
optio = "RNOQ"; |
|
|
if ( draw ) optio = "ROQ"; |
|
1244 |
}; |
}; |
1245 |
// |
// |
1246 |
if ( debug ) printf(" OK, start the fitting procedure...\n"); |
if ( debug ) printf(" OK, start the fitting procedure...\n"); |
1247 |
// |
// |
1248 |
fitresult = th->Fit("lfit",optio); |
// fitresult = th->Fit("lfit",optio); |
1249 |
|
fitresult = gh->Fit("lfit",optio); |
1250 |
// |
// |
1251 |
if ( debug ) printf(" the fit is done! result: %i \n",fitresult); |
if ( debug ) printf(" the fit is done! result: %i \n",fitresult); |
1252 |
// |
// |
1265 |
if ( fitresult != 0 ){ |
if ( fitresult != 0 ){ |
1266 |
if ( debug ) printf(" The fit failed, no integrals calculation and asymm is set to -1. \n"); |
if ( debug ) printf(" The fit failed, no integrals calculation and asymm is set to -1. \n"); |
1267 |
asymm = -1.; |
asymm = -1.; |
1268 |
|
defE0 = -1.; |
1269 |
} else { |
} else { |
1270 |
|
if ( slmax.MaybeRegexp() ) lmax = this->Evaluate(slmax,tmax,X0pl); |
1271 |
|
if ( sumax.MaybeRegexp() ) umax = this->Evaluate(sumax,tmax,X0pl); |
1272 |
Int_t npp = 1000; |
Int_t npp = 1000; |
1273 |
|
double *xpp=new double[npp]; |
1274 |
|
double *wpp=new double[npp]; |
1275 |
|
lfit->CalcGaussLegendreSamplingPoints(npp,xpp,wpp,1e-12); |
1276 |
|
defE0 = lfit->IntegralFast(npp,xpp,wpp,lmax,umax); |
1277 |
|
// |
1278 |
double *xp=new double[npp]; |
double *xp=new double[npp]; |
1279 |
double *wp=new double[npp]; |
double *wp=new double[npp]; |
1280 |
lfit->CalcGaussLegendreSamplingPoints(npp,xp,wp,1e-12); |
lfit->CalcGaussLegendreSamplingPoints(npp,xp,wp,1e-12); |
1323 |
// axis titles |
// axis titles |
1324 |
th->SetXTitle("Depth [X0]"); |
th->SetXTitle("Depth [X0]"); |
1325 |
th->SetYTitle("Energy [MIP]"); |
th->SetYTitle("Energy [MIP]"); |
1326 |
th->DrawCopy("Perror"); |
// th->DrawCopy("Perror"); |
1327 |
|
th->DrawCopy(); |
1328 |
|
gh->SetMarkerSize(1.); |
1329 |
|
gh->SetMarkerStyle(21); |
1330 |
|
// gh->SetMarkerColor(kRed); |
1331 |
|
gh->Draw("Psame"); |
1332 |
lfit->Draw("same"); |
lfit->Draw("same"); |
1333 |
tc->Modified(); |
tc->Modified(); |
1334 |
tc->Update(); |
tc->Update(); |
1340 |
if ( th ) th->Delete(); |
if ( th ) th->Delete(); |
1341 |
}; |
}; |
1342 |
// |
// |
1343 |
delete lfit; |
// delete lfit; |
1344 |
// |
// |
1345 |
}; |
}; |
1346 |
|
|
1431 |
Clear(); |
Clear(); |
1432 |
//delete this; |
//delete this; |
1433 |
}; |
}; |
1434 |
|
|
1435 |
|
|
1436 |
|
|
1437 |
|
|
1438 |
|
|
1439 |
|
|
1440 |
|
|
1441 |
|
|
1442 |
|
|
1443 |
|
|
1444 |
|
|
1445 |
|
|
1446 |
|
|
1447 |
|
|
1448 |
|
|
1449 |
|
|
1450 |
|
|
1451 |
|
|
1452 |
|
|
1453 |
|
|