--- calo/flight/CaloProfile/src/CaloProfile.cpp 2007/11/09 09:11:24 1.1 +++ calo/flight/CaloProfile/src/CaloProfile.cpp 2010/05/13 13:55:37 1.18 @@ -1,4 +1,209 @@ #include +// +ClassImp(CaloLat); +ClassImp(CaloLong); +ClassImp(Calo2D); + + +int isempty(struct stack *s){ + return (s->top == EMPTY) ? 1 : 0; +} + +void emptystack(struct stack* s){ + s->top=EMPTY; + strcpy(s->data," "); +} + +void push(struct stack* s,int item){ + if(s->top == (MAX-1)){ + printf("\n ERROR! STACK FULL (too many digits in SetLower/UpperLimit)\n"); + } else { + s->data[++s->top]=(char)item; + }; +} + +char pop(struct stack* s){ + char ret=(char)EMPTY; + if(!isempty(s)){ + ret= s->data[s->top--]; + }; + return ret; +} + + +int ipop(struct stack* s){ + int ret=EMPTY; + if(s->top == EMPTY) + printf("\n ERROR! STACK EMPTY (too few digits in SetLower/UpperLimit)\n"); + else { + ret= s->data[s->top--]; + }; + return ret; +} + + +void display(struct stack s){ + int ss = s.top; + while(s.top != EMPTY){ + printf("\n%d",s.data[s.top--]); + }; + printf(" s.top %i \n",ss); +} + +int isoperator(char e){ + if(e == '+' || e == '-' || e == '*' || e == '/' || e == '%') + return 1; + else + return 0; +} + + +int priority(char e){ + int pri = 0; + if(e == '*' || e == '/' || e =='%') + pri = 2; + else { + if(e == '+' || e == '-') pri = 1; + }; + return pri; +} + +void infix2postfix(char* infix, char * postfix, int insertspace){ + char *i,*p; + struct stack X; + char n1; + emptystack(&X); + i = &infix[0]; + p = &postfix[0]; + while(*i){ + while(*i == ' ' || *i == '\t' ){ + i++; + }; + TString c=i; + if( isdigit(*i) || isalpha(*i) || c.BeginsWith(".") ){ + while( isdigit(*i) || isalpha(*i) || c.BeginsWith(".") ){ + *p = *i; + p++; + i++; + c=i; + }; + /*SPACE CODE*/ + if(insertspace){ + *p = ' '; + p++; + }; + /*END SPACE CODE*/ + }; + if( *i == '(' ){ + push(&X,*i); + i++; + }; + if( *i == ')'){ + n1 = pop(&X); + while( n1 != '(' ){ + *p = n1; + p++; + /*SPACE CODE*/ + if(insertspace){ + *p = ' '; + p++; + }; + /*END SPACE CODE*/ + n1 = pop(&X); + }; + i++; + }; + if( isoperator(*i) ){ + if(isempty(&X)){ + push(&X,*i); + } else { + n1 = pop(&X); + while( priority(n1) >= priority(*i) ){ + *p = n1; + p++; + /*SPACE CODE*/ + if(insertspace){ + *p = ' '; + p++; + }; + /*END SPACE CODE*/ + n1 = pop(&X); + }; + push(&X,n1); + push(&X,*i); + }; + i++; + }; + }; + while(!isempty(&X)){ + n1 = pop(&X); + *p = n1; + p++; + /*SPACE CODE*/ + if(insertspace){ + *p = ' '; + p++; + }; + /*END SPACE CODE*/ + }; + *p = '\0'; +} + + +Float_t evaluate(char *postfix){ + // + Float_t op1 = 0.; + Float_t op2 = 0.; + Float_t result = 0.; + // + TString e = postfix; + Float_t st[50]; + memset(st,0,50*sizeof(Float_t)); + TObjArray *ae = e.Tokenize(" "); + Int_t o = 0; + Int_t a = 0; + Int_t ap = 0; + // + while ( o < ae->GetEntries() ){ + if ( ((TString)(((TObjString*)ae->At(o))->GetString())).IsFloat() ){ + st[a] = (Float_t)((TString)(((TObjString*)ae->At(o))->GetString())).Atof();; + a++; + } else { + ap = a-1; + op1 = st[ap--]; + op2 = st[ap]; + const char *p=((TString)(((TObjString*)ae->At(o))->GetString())).Data(); + switch(*p){ + case '+': + result = op2 + op1; + break; + case '-': + result = op2 - op1; + break; + case '/': + result = op2 / op1; + break; + case '*': + result = op2 * op1; + break; + case '%': + result = (Int_t)round(op2) % (Int_t)round(op1); + break; + default: + printf("\nInvalid Operator: %s \n",((TString)(((TObjString*)ae->At(o))->GetString())).Data()); + return 0; + }; + st[ap] = result; + a = ap+1; + // + }; + o++; + }; + return result; + // +} + + //////////////////////////////////////////////////////////////////////// /** @@ -60,7 +265,7 @@ /** * Default constructor */ -CaloLong::CaloLong(){ +Calo2D::Calo2D(){ Clear(); }; @@ -76,11 +281,13 @@ PKT = 0; atime = 0; // + suf = ""; debug = false; + usepl18x = false; // }; -CaloLong::CaloLong(PamLevel2 *l2p){ +Calo2D::Calo2D(PamLevel2 *l2p){ // Clear(); // @@ -92,7 +299,9 @@ PKT = 0; atime = 0; // + suf = ""; debug = false; + usepl18x = false; // }; @@ -100,7 +309,7 @@ // }; -void CaloLong::Clear(){ +void Calo2D::Clear(){ // }; @@ -116,11 +325,11 @@ // }; -void CaloLong::Print(){ +void Calo2D::Print(){ // Process(); // - printf("==================== Calorimeter Lateral Profile =======================\n"); + printf("==================== Calorimeter 2D Profile =======================\n"); printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime); // printf(" nx [number of X combination]:.. %i\n",nx); // printf(" ny [number of Y combination]:.. %i\n",ny); @@ -134,13 +343,12 @@ Draw(-1,-1); }; -void CaloLong::Draw(){ +void Calo2D::Draw(){ // Process(); Draw(-1); }; - void CaloLat::Draw(Int_t view,Int_t plane){ // Int_t minv = 0; @@ -171,7 +379,7 @@ // for (Int_t v=minv; v(gDirectory->FindObject(hid)); if ( tc ){ // tc->Clear(); @@ -179,7 +387,7 @@ tc = new TCanvas(hid,hid); }; // - TString thid = Form("hlatv%ip%i",v,p); + TString thid = Form("hlatv%ip%i%s",v,p,suf.Data()); TH1F *th = dynamic_cast(gDirectory->FindObject(thid)); if ( th ) th->Delete(); // th->Clear(); @@ -203,16 +411,17 @@ // }; -void CaloLong::Draw(Int_t view){ +void Calo2D::Draw(Int_t plane){ // - Int_t minv = 0; - Int_t maxv = 0; + Int_t minp = 0; + Int_t maxp = 0; // - if ( view == -1 ){ - maxv = -1; + if ( plane == -1 ){ + minp = 0; + maxp = 23; } else { - minv = view; - maxv = view+1; + minp = plane; + maxp = plane+1; }; // Process(); @@ -220,59 +429,51 @@ gStyle->SetLabelSize(0.04); gStyle->SetNdivisions(510,"XY"); // - if ( maxv != -1 ){ - for (Int_t v=minv; v(gDirectory->FindObject(hid)); - if ( tc ){ -// tc->Clear(); - } else { - tc = new TCanvas(hid,hid); - }; - // - TString thid = Form("hlongv%i",v); - TH1F *th = dynamic_cast(gDirectory->FindObject(thid)); - if ( th ) th->Delete(); -// th->Clear(); -// th->Reset(); -// } else { - th = new TH1F(thid,thid,22,-0.5,21.5); -// }; - tc->cd(); - // - for (Int_t st=0;st<22;st++){ - th->Fill(st,eplane[v][st]); - }; - th->Draw(); - tc->Modified(); - tc->Update(); - }; - } else { - // - TString hid = Form("clongvyvx"); - TCanvas *tc = dynamic_cast(gDirectory->FindObject(hid)); - if ( tc ){ - } else { - tc = new TCanvas(hid,hid); - }; - // - TString thid = Form("hlongvyvx"); - TH1F *th = dynamic_cast(gDirectory->FindObject(thid)); - if ( th ) th->Delete(); - th = new TH1F(thid,thid,44,-0.5,43.5); - tc->cd(); - Int_t pp=0; - for (Int_t st=0;st<22;st++){ - for (Int_t v=1; v>=0;v--){ - // - th->Fill(pp,eplane[v][st]); - // - pp++; - }; - }; - th->Draw(); - tc->Modified(); - tc->Update(); + for (Int_t p=minp; p(gDirectory->FindObject(hid)); + if ( tc ){ + // tc->Clear(); + } else { + tc = new TCanvas(hid,hid); + }; + // + TString thid = Form("h2dp%i%s",p,suf.Data()); + TH2F *th = dynamic_cast(gDirectory->FindObject(thid)); + if ( th ) th->Delete(); + // th->Clear(); + // th->Reset(); + // } else { + Int_t minx = smax[p] - 10; + if ( minx < 0 ) minx = 0; + Int_t maxx = minx + 20; + if ( maxx > 95 ){ + maxx = 95; + minx = 75; + }; + Int_t miny = smay[p] - 10; + if ( miny < 0 ) miny = 0; + Int_t maxy = miny + 20; + if ( maxy > 95 ){ + maxy = 95; + miny = 75; + }; + th = new TH2F(thid,thid,20,(Float_t)minx-0.5,(Float_t)maxx-0.5,20,(Float_t)miny-0.5,(Float_t)maxy-0.5); + // th = new TH2F(thid,thid,96,-0.5,95.5,96,-0.5,95.5); + // }; + tc->cd(); + // + for (Int_t stx=minx;stxFill(stx,sty,estrip[p][stx][sty]); + }; + }; + gStyle->SetPalette(1); + // tc->SetLogz(); + // th->Draw("colbox"); + th->Draw("cont4"); + tc->Modified(); + tc->Update(); }; // gStyle->SetLabelSize(0); @@ -285,11 +486,12 @@ //delete this; }; -void CaloLong::Delete(){ +void Calo2D::Delete(){ Clear(); //delete this; }; + void CaloLat::Process(){ // if ( !L2 ){ @@ -327,6 +529,9 @@ // for (Int_t i=0; iGetCaloLevel1()->istrip ; i++){ mip1 = L2->GetCaloLevel1()->DecodeEstrip(i,view1,plane1,strip1); + // + if ( !usepl18x && view1==0 && plane1==18 ) mip1 = 0.; + // estrip[view1][plane1][strip1] = mip1; }; // @@ -335,7 +540,7 @@ // }; -void CaloLong::Process(){ +void Calo2D::Process(){ // if ( !L2 ){ printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n"); @@ -364,16 +569,951 @@ // // let's start // + Float_t es[2][22][96]; + memset(es,0, 4224*sizeof(Float_t)); + memset(estrip,0, 4224*sizeof(Float_t)); + Float_t mip1 = 0.; + Int_t view1 = 0; + Int_t plane1 = 0; + Int_t strip1 = 0; + // + for (Int_t i=0; iGetCaloLevel1()->istrip ; i++){ + mip1 = L2->GetCaloLevel1()->DecodeEstrip(i,view1,plane1,strip1); + // + if ( !usepl18x && view1==0 && plane1==18 ) mip1 = 0.; + // + es[view1][plane1][strip1] = mip1; + }; + // + Int_t plane2 = 0; + Float_t emax[23]; + memset(emax,0,sizeof(Float_t)*23); + memset(smax,0,sizeof(Int_t)*23); + memset(smay,0,sizeof(Int_t)*23); + // + // + for (Int_t p=0; p < 23 ; p++){ + // + plane1 = p-1; + plane2 = p; + // + if ( p == 0 ){ + plane1 = -1; + plane2 = 0; + }; + if ( p == 22 ){ + plane1 = 21; + plane2 = -1; + }; + // + for (Int_t s=0; s < 96 ; s++){ // x + for (Int_t ss=0; ss < 96 ; ss++){ // y + if ( p == 0 ){ + estrip[p][s][ss] += es[1][plane2][ss]; + if ( (es[1][plane2][ss]) > emax[p] ){ + smax[p] = 45; + smay[p] = ss; + emax[p] = es[1][plane2][ss] ; + }; + }; + if ( p > 0 && p < 22 ){ + estrip[p][s][ss] += es[0][plane1][s] + es[1][plane2][ss]; + if ( (es[0][plane1][s] + es[1][plane2][ss]) > emax[p] ){ + smax[p] = s; + smay[p] = ss; + emax[p] = es[0][plane1][s] + es[1][plane2][ss] ; + }; + }; + if ( p == 22 ){ + estrip[p][s][ss] += es[0][plane1][s]; + if ( (es[1][plane2][s]) > emax[p] ){ + smax[p] = s; + smay[p] = 45; + emax[p] = es[1][plane2][s] ; + }; + }; + }; + }; + // + }; + // + if ( debug ) this->Print(); + if ( debug ) printf(" exit \n"); + // +}; + + +/** + * Default constructor + */ +CaloLong::CaloLong(){ + Clear(); +}; + +CaloLong::CaloLong(PamLevel2 *l2p){ + // + Clear(); + // + L2 = l2p; + // + if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n"); + // + OBT = 0; + PKT = 0; + atime = 0; + // + clp = L2->GetCaloLevel2(); + // + sel = true; + cont = false; + N = 0; + NC = 22; + mask18b = -1; + // + no18x = true; + usepl18x = !no18x; + debug = false; + maskXE = false; + maskXO = false; + maskYE = false; + maskYO = false; + // + lmax = 0.; + umax = 100.; + slmax = ""; + sumax = ""; + suf = ""; + xyaverage = true; + // + letmax = -1; + heavytail = false; + // lmipth = 100.; + lmipth = 0.; + // +}; + +TF1 *CaloLong::GetFit(){ + TString fnam=Form("lfit%s",suf.Data()); + TF1 *lfit = dynamic_cast(gDirectory->FindObject(fnam)); + if ( lfit ) return lfit; + return NULL; +} + +void CaloLong::MaskSection(TString sec){ + sec.ToUpper(); + if ( sec.Contains("XO") ) maskXO = true; + if ( sec.Contains("YO") ) maskYO = true; + if ( sec.Contains("XE") ) maskXE = true; + if ( sec.Contains("YE") ) maskYE = true; +} + +void CaloLong::UnMaskSections(){ + this->UnMaskSection("XEXOYEYO"); +} + +void CaloLong::UnMaskSection(TString sec){ + sec.ToUpper(); + if ( sec.Contains("XO") ) maskXO = false; + if ( sec.Contains("YO") ) maskYO = false; + if ( sec.Contains("XE") ) maskXE = false; + if ( sec.Contains("YE") ) maskYE = false; +} + +void CaloLong::Clear(){ + // + if ( debug ) printf(" Clear called \n"); + // memset(eplane,0, 2*22*sizeof(Float_t)); // + chi2 = 0.; + ndf = 0.; + E0 = 0.; + defE0 = 0.; + a = 0.; + b = 0.; + errE0 = 0.; + erra = 0.; + errb = 0.; + etmax = 0.; + asymm = 0.; + fitresult = 0; + // + X0pl = 0.76; + // +}; + +void CaloLong::Print(){ + // + Fit(); + // + printf("==================== Calorimeter Longitudinal Profile =======================\n"); + printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime); + printf(" fitresult:.. %i\n",fitresult); + printf(" chi2 :.. %f\n",chi2); + printf(" ndf :.. %f\n",ndf); + printf(" nchi2 :.. %f\n",chi2/ndf); + printf(" E0 :.. %f\n",E0); + printf(" E0/260. :.. %f\n",E0/260.); + printf(" defE0 :.. %f\n",defE0); + printf(" lower :.. %f\n",lmax); + printf(" upper :.. %f\n",umax); + printf(" s lower :.. %s\n",slmax.Data()); + printf(" s upper :.. %s\n",sumax.Data()); + printf(" a :.. %f\n",a); + printf(" b :.. %f\n",b); + printf(" errE0 :.. %f\n",errE0); + printf(" erra :.. %f\n",erra); + printf(" errb :.. %f\n",errb); + printf(" asymm :.. %f\n",asymm); + printf(" tmax :.. %f\n",((a-1.)/b)); + printf(" etmax :.. %f\n",etmax); + printf(" X0pl :.. %f\n",X0pl); + printf("========================================================================\n"); + // +}; + +void CaloLong::SetNoWpreSampler(Int_t n){ + // + if ( NC+n <= 22 && NC+n >= 0 ){ + N = n; + } else { + printf(" ERROR! Calorimeter is made of 22 W planes\n"); + printf(" you are giving N presampler = %i and N calo = %i \n",n,NC); + printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n"); + NC = 22; + N = 0; + }; +} + +void CaloLong::SetNoWcalo(Int_t n){ + if ( N+n <= 22 && N+n >= 0 ){ + NC = n; + } else { + printf(" ERROR! Calorimeter is made of 22 W planes\n"); + printf(" you are giving N W presampler = %i and N W calo = %i \n",N,n); + printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n"); + NC = 22; + N = 0; + }; +} + +void CaloLong::SplitInto(Int_t NoWpreSampler, Int_t NoWcalo){ + this->SetNoWpreSampler(0); + this->SetNoWcalo(0); + if ( NoWpreSampler < NoWcalo ){ + this->SetNoWpreSampler(NoWpreSampler); + this->SetNoWcalo(NoWcalo); + } else { + this->SetNoWcalo(NoWcalo); + this->SetNoWpreSampler(NoWpreSampler); + }; +} + +void CaloLong::Process(){ + // + if ( !L2 ){ + printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n"); + printf(" ERROR: CaloHough variables not filled \n"); + return; + }; + // + Bool_t newentry = false; + // + if ( L2->IsORB() ){ + if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){ + newentry = true; + OBT = L2->GetOrbitalInfo()->OBT; + PKT = L2->GetOrbitalInfo()->pkt_num; + atime = L2->GetOrbitalInfo()->absTime; + }; + } else { + newentry = true; + }; + // + if ( !newentry ) return; + // + if ( debug ) printf(" Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime); + // + Clear(); + // + // let's start + // + if ( cont ){ + for (Int_t i=0; i<22; i++){ + if ( i == (18+N) ){ + mask18b = 18 + N; + break; + }; + }; + }; + // + if ( sel ){ + for (Int_t i=0; i<22; i++){ + if ( i == (18-N) ){ + mask18b = 18 - N; + break; + }; + }; + }; + // + // if ( mask18b == 18 ) mask18b = -1; + // + Int_t view = 0; + Int_t plane = 0; + Int_t strip = 0; + Float_t mip = 0.; + Bool_t gof = true; + for (Int_t i=0; i < L2->GetCaloLevel1()->istrip; i++){ + mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); + // + if ( !usepl18x && view==0 && plane==18 ) mip = 0.; + // + gof = true; + if ( maskXE && (plane%2)==0 && view==1 ) gof = false; + if ( maskXO && (plane%2)!=0 && view==1 ) gof = false; + if ( maskYE && (plane%2)!=0 && view==0 ) gof = false; + if ( maskYO && (plane%2)==0 && view==0 ) gof = false; + if ( gof ) eplane[view][plane] += mip; + }; + // + // inclination factor (stolen from Daniele's code) + // + Float_t ytgx = 0.; + Float_t ytgy = 0.; + ytgx = 0.76 * clp->tanx[0]; + ytgy = 0.76 * clp->tany[0]; + X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) ); + // + // Find experimental plane of maximum + // + Int_t pmax = 0; + Int_t vmax = 0; + Float_t emax = 0.; for (Int_t v=0; v<2; v++){ for (Int_t i=0; i<22; i++){ - eplane[v][i] = L2->GetCaloLevel1()->qtotpl(v,i); + if ( eplane[v][i] > emax ){ + emax = eplane[v][i]; + vmax = v; + pmax = i; + }; }; }; // + // + // + if ( vmax == 0 ) pmax++; + etmax = pmax * X0pl; + // if ( debug ) this->Print(); if ( debug ) printf(" exit \n"); // }; +void CaloLong::SetEnergies(Float_t myene[][22]){ + // + if ( !L2 ){ + printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n"); + printf(" ERROR: CaloHough variables not filled \n"); + return; + }; + // + Bool_t newentry = false; + // + if ( L2->IsORB() ){ + if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){ + newentry = true; + OBT = L2->GetOrbitalInfo()->OBT; + PKT = L2->GetOrbitalInfo()->pkt_num; + atime = L2->GetOrbitalInfo()->absTime; + }; + } else { + newentry = true; + }; + // + if ( !newentry ) return; + // + if ( debug ) printf(" SET ENERGIES Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime); + // + Clear(); + // + // let's start + // + if ( cont ){ + for (Int_t i=0; i<22; i++){ + if ( i == (18+N) ){ + mask18b = 18 + N; + break; + }; + }; + }; + // + if ( sel ){ + for (Int_t i=0; i<22; i++){ + if ( i == (18-N) ){ + mask18b = 18 - N; + break; + }; + }; + }; + // + // if ( mask18b == 18 ) mask18b = -1; + // + Int_t view = 0; + Int_t plane = 0; + Bool_t gof = true; + for (view=0; view < 2; view++){ + for (plane=0; plane < 22; plane++){ + gof = true; + if ( maskXE && (plane%2)==0 && view==1 ) gof = false; + if ( maskXO && (plane%2)!=0 && view==1 ) gof = false; + if ( maskYE && (plane%2)!=0 && view==0 ) gof = false; + if ( maskYO && (plane%2)==0 && view==0 ) gof = false; + if ( gof ) eplane[view][plane] = myene[view][plane]; + 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]); + }; + }; + // + if ( !clp && debug ) printf(" ERROR!! no CaloLevel2 obj!!\n"); + // + // inclination factor (stolen from Daniele's code) + // + Float_t ytgx = 0.; + Float_t ytgy = 0.; + ytgx = 0.76 * clp->tanx[0]; + ytgy = 0.76 * clp->tany[0]; + X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) ); + if ( debug ) printf(" X0pl %f tanx %f tany %f \n",X0pl,ytgx,ytgy); + // + // Find experimental plane of maximum + // + Int_t pmax = 0; + Int_t vmax = 0; + Float_t emax = 0.; + for (Int_t v=0; v<2; v++){ + for (Int_t i=0; i<22; i++){ + if ( eplane[v][i] > emax ){ + emax = eplane[v][i]; + vmax = v; + pmax = i; + }; + }; + }; + // + // + // + if ( vmax == 0 ) pmax++; + etmax = pmax * X0pl; + letmax = etmax; + // + if ( debug ) this->Print(); + if ( debug ) printf(" exit \n"); + // +}; + +Double_t ccurve(Double_t *ti,Double_t *par){ + // + Double_t t = *ti; + Double_t cE0 = par[0]; + Double_t ca = par[1]; + Double_t cb = par[2]; + Double_t gammaa = TMath::Gamma(ca); + // + Double_t value = cE0 * cb * ( pow((cb * t),(ca - 1)) * exp( -cb * t ) ) / gammaa; + // + return value; + // +} + +void CaloLong::Fit(){ + this->Fit(false); +}; + +Float_t CaloLong::Evaluate(TString s, Float_t tmax, Float_t X0pl){ + /* SAMPLE OUTPUT: + Enter Infix Expression : A + B + C / (E - F) + Postfix Expression is : A B + C E F - / + + */ + if ( !s.Contains("tmax") && !s.Contains("X0pl") ){ + printf(" ERROR, the input formula must contain \"t\"\n"); + return 0.; + }; + if ( tmax != tmax ){ + printf(" ERROR, tmax is nan! \n"); + return 0.; + }; + TString g=Form("%f",tmax); + TString h=Form("%f",X0pl); + TString *ts= new TString(""); + ts->Prepend(s.Data()); + ts->ReplaceAll("tmax",4,g.Data(),g.Capacity()); + ts->ReplaceAll("X0pl",4,h.Data(),h.Capacity()); + ts->Prepend("("); + ts->Append(")"); + if ( debug ) printf(" ts %s tssize %i capac %i s %s g %s \n",ts->Data(),ts->Sizeof(),ts->Capacity(),s.Data(),g.Data()); + char in[50],post[50]; + strcpy(&in[0]," "); + strcpy(&in[0],ts->Data()); + strcpy(&post[0]," "); + if ( debug ) printf("Infix Expression is : ##%s##\n",&in[0]); + infix2postfix(&in[0],&post[0],1); + if ( debug ) printf("Postfix Expression is : ##%s##\n",&post[0]); + /* SAMPLE OUTPUT: + Enter Postfix Expression : 3 5 + 2 / + 3 5 + 2 / EQUALS 4 + */ + Float_t res = evaluate(&post[0]); + if ( debug ) printf("%s EQUALS %f\n",post,res); + return res; +} + +void CaloLong::Fit(Bool_t draw){ + // + Process(); + // + // + if ( !L2 ){ + printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n"); + printf(" ERROR: CaloHough variables not filled \n"); + return; + }; + // + Bool_t newentry = false; + // + if ( L2->IsORB() ){ + if ( L2->GetOrbitalInfo()->pkt_num != fPKT || L2->GetOrbitalInfo()->OBT != fOBT || L2->GetOrbitalInfo()->absTime != fatime ){ + newentry = true; + fOBT = L2->GetOrbitalInfo()->OBT; + fPKT = L2->GetOrbitalInfo()->pkt_num; + fatime = L2->GetOrbitalInfo()->absTime; + }; + } else { + newentry = true; + }; + // + if ( !newentry ) return; + // + if ( debug ) printf(" Start fitting event at OBT %u PKT %u time %u \n",fOBT,fPKT,fatime); + // + if ( draw ){ + gStyle->SetLabelSize(0.04); + gStyle->SetNdivisions(510,"XY"); + }; + // + TString hid = Form("clongfit%s",suf.Data()); + TCanvas *tc = dynamic_cast(gDirectory->FindObject(hid)); + // if ( tc ) tc->Delete(); + // if ( tc ) tc->Close(); + if ( !tc && draw ){ + tc = new TCanvas(hid,hid); + } else { + if ( tc ) tc->cd(); + }; + // + TString thid = Form("hlongfit%s",suf.Data()); + TH2F *th = dynamic_cast(gDirectory->FindObject(thid)); + if ( th ) th->Delete(); + // + TString ghid = Form("glongfit%s",suf.Data()); + TGraphErrors *gh = dynamic_cast(gDirectory->FindObject(ghid)); + if ( gh ) gh->Delete(); + // + Float_t xpos = 0.; + Float_t enemip = 0.; + Float_t xmax = NC * X0pl + 0.2; + // + Double_t xxx[50]; + Double_t exx[50]; + Double_t yyy[50]; + Double_t eyy[50]; + Int_t numpo = 0; + memset(xxx,0,50*sizeof(Double_t)); + memset(exx,0,50*sizeof(Double_t)); + memset(yyy,0,50*sizeof(Double_t)); + memset(eyy,0,50*sizeof(Double_t)); + // + // AGH, BUG! + // + Int_t mmin = 0; + Int_t mmax = 0; + if ( cont ){ + mmin = N; + mmax = NC+N; + } else { + mmin = 0; + mmax = NC; + }; + // + Float_t emax = 0.; + Float_t qtotparz = 0.; + for (Int_t st=mmin;st mmin && st < mmax ){ + if ( no18x && ( st == 18+1 || st == mask18b+1 )){ + if ( !maskYO ){ + enemip = 2. * eplane[1][st]; + } else { + enemip = eplane[1][st]; + }; + } else { + enemip = eplane[0][st-1] + eplane[1][st]; + }; + } else { + if ( st == mmin ){ + if ( !maskYE ){ + enemip = 2. * eplane[1][st]; + } else { + enemip = eplane[1][st]; + }; + }; + if ( st == mmax ){ + if ( !maskXE ){ + enemip = 2. * eplane[0][st-1]; + } else { + enemip = eplane[0][st-1]; + }; + }; + }; + // + qtotparz += enemip; + if ( enemip > emax ) emax = enemip; + if ( enemip > 0. ){ + xxx[numpo] = xpos; + exx[numpo] = 0.1; + yyy[numpo] = enemip; + eyy[numpo] = sqrt(enemip*3.)+sqrt(5.); + //if ( xpos > letmax && enemip > lmipth && heavytail) eyy[numpo] = (sqrt(enemip*3.)+sqrt(5.))/numpo; + if ( xpos > letmax && enemip > lmipth && heavytail) eyy[numpo] = sqrt(enemip)/5.; + if ( xpos > letmax-1 && xpos < letmax+1 && heavytail ) eyy[numpo] /= 5.; + //if ( xpos > letmax-2 && xpos < letmax+1 && heavytail ) eyy[numpo] /= 7.; + if ( xpos < 3. && heavytail ) eyy[numpo] /= 5.; + numpo++; + // th->Fill(xpos,enemip); + if ( debug ) printf(" AVE Filling: st %i xpos %f energy %f error %f letmax %f lmpith %f \n",st,xpos,enemip,eyy[numpo-1],letmax,lmipth); + }; + } else { + for (Int_t jj=0; jj<2; jj++){ + if ( st > mmin && st < mmax ){ + if ( jj == 0 && no18x && ( st == 18+1 || st == mask18b+1 )){ + enemip = eplane[1][st]; + } else { + if ( jj == 0 ){ + enemip = eplane[jj][st-1]; + } else { + enemip = eplane[jj][st]; + }; + }; + } else { + if ( st == mmin && jj == 1 ){ + enemip = eplane[jj][st]; + }; + if ( st == mmax && jj == 0){ + enemip = eplane[jj][st-1]; + }; + }; + // + qtotparz += enemip; + if ( enemip > emax ) emax = enemip; + if ( enemip > 0. ){ + xxx[numpo] = xpos; + exx[numpo] = 0.1; + yyy[numpo] = enemip; + eyy[numpo] = sqrt(enemip*3.)+sqrt(5.); + //if ( xpos > letmax && enemip > lmipth && heavytail) eyy[numpo] = (sqrt(enemip*3.)+sqrt(5.))/numpo; + //if ( xpos > letmax && enemip > lmipth && heavytail ) eyy[numpo] = (sqrt(enemip*3.)+sqrt(5.))/numpo; + if ( xpos > letmax && enemip > lmipth && heavytail) eyy[numpo] = sqrt(enemip)/5.; + if ( xpos > letmax-1 && xpos < letmax+1 && heavytail ) eyy[numpo] /= 5.; + //if ( xpos > letmax-2 && xpos < letmax+1 && heavytail ) eyy[numpo] /= 7.; + if ( xpos < 3. && heavytail ) eyy[numpo] /= 5.; + // eyy[numpo] = sqrt(enemip)/(st*0.95); + numpo++; + // th->Fill(xpos,enemip); + if ( debug ) printf(" XY Filling: st %i xpos %f energy %f \n",st,xpos,enemip); + }; + }; + }; + + // + // for (Int_t v=1; v>=0;v--)// { +// // +// if ( v == 1 ){ +// xpos = (st - N) * X0pl; +// } else { +// xpos = (st + 1 - N) * X0pl; +// }; +// // +// if ( no18x && st == 18 && v == 0 ){ +// // skip plane 18x +// } else { +// if ( v == 1 && st == mask18b ){ +// // emulate plane 18x +// } else { +// if ( eplane[v][st] > 0. ){ +// th->Fill(xpos,eplane[v][st]); +// if ( debug ) printf(" Filling: st %i v %i xpos %f energy %f \n",st,v,xpos,eplane[v][st]); +// }; +// }; +// }; +// // +// }; + }; + // + // th = new TH2F(thid,thid,int(NC*1.5),-0.2,xmax); + th = new TH2F(thid,thid,1000,-0.2,xmax,1000,0.,emax*1.2); + gh = new TGraphErrors(numpo,xxx,yyy,exx,eyy); + TString fnam=Form("lfit%s",suf.Data()); + TF1 *lfit = dynamic_cast(gDirectory->FindObject(fnam)); + if ( lfit ) lfit->Delete(); + lfit = new TF1(fnam,ccurve,0.,xmax,3); +// if ( !lfit ){ +// lfit = new TF1(fnam,ccurve,0.,xmax,3); +// }; +// lfit->Clear(); +// lfit->SetName(fnam); +// lfit->SetRange(0.,xmax); +// lfit->Print(); + // + if ( debug ) printf("qtot %f qtotparz %f \n",clp->qtot,qtotparz); + E0 = qtotparz; + // E0 = clp->qtot; + a = 5.; + b = 0.5; + if ( debug ) printf(" STARTING PARAMETERS: E0 %f a %f b %f \n",E0,a,b); + lfit->SetParameters(E0,a,b); + // lfit->SetParLimits(0,0.,1000.); + // lfit->SetParLimits(1,-1.,80.); + // lfit->SetParLimits(2,-1.,10.); + // TString optio = "ROW"; // "RO" + TString optio = "RO"; // "RO" + if ( debug ){ + optio += "NV"; + if ( draw ) optio += "V"; + } else { + optio += "NQ"; + if ( draw ) optio += "Q"; + }; + // + if ( debug ) printf(" OK, start the fitting procedure... options: %s \n",optio.Data()); + // + // fitresult = th->Fit("lfit",optio); + // lfit->Update(); + fitresult = gh->Fit("lfit",optio); + // fitresult = gh->Fit("lfit","ROQW"); + // + if ( debug ) printf(" the fit is done! result: %i \n",fitresult); + // + E0 = lfit->GetParameter(0); + a = lfit->GetParameter(1); + b = lfit->GetParameter(2); + errE0 = lfit->GetParError(0); + erra = lfit->GetParError(1); + errb = lfit->GetParError(2); + chi2 = lfit->GetChisquare(); + ndf = lfit->GetNDF(); + Float_t tmax = 0.; + if ( debug ) printf(" Parameters are retrieved \n"); + if ( b != 0 ) tmax = (a - 1.)/b; + // + if ( fitresult != 0 ){ + if ( debug ) printf(" The fit failed, no integrals calculation and asymm is set to -1. \n"); + asymm = -1.; + defE0 = -1.; + } else { + if ( slmax.MaybeRegexp() ) lmax = this->Evaluate(slmax,tmax,X0pl); + if ( sumax.MaybeRegexp() ) umax = this->Evaluate(sumax,tmax,X0pl); + Int_t npp = 1000; + double *xpp=new double[npp]; + double *wpp=new double[npp]; + lfit->CalcGaussLegendreSamplingPoints(npp,xpp,wpp,1e-12); + defE0 = lfit->IntegralFast(npp,xpp,wpp,lmax,umax); + // + double *xp=new double[npp]; + double *wp=new double[npp]; + lfit->CalcGaussLegendreSamplingPoints(npp,xp,wp,1e-12); + Float_t imax = lfit->IntegralFast(npp,xp,wp,0.,tmax); + // Float_t imax = lfit->Integral(0.,tmax); + if ( debug ) printf(" Integral till maximum (%f): %f \n",tmax,imax); + Int_t np = 1000; + double *x=new double[np]; + double *w=new double[np]; + lfit->CalcGaussLegendreSamplingPoints(np,x,w,1e-12); + Float_t i10max = lfit->IntegralFast(np,x,w,0.,10.*tmax); + delete x; + delete w; + delete xp; + delete wp; + delete xpp; + delete wpp; + // Float_t i10max = lfit->Integral(0.,10.*tmax); + if ( debug ) printf(" Integral: %f \n",i10max); + // + if ( i10max != imax ){ + asymm = imax / (i10max-imax); + } else { + if ( debug ) printf(" i10max == imax, asymm undefined\n"); + asymm = -2.; + }; + if ( asymm != asymm ){ + if ( debug ) printf(" asymm is nan \n"); + asymm = -3.; + }; + //lfit->Integral(0.,tmax)/(lfit->Integral(0.,10.*tmax)-lfit->Integral(0.,tmax)); + if ( debug ) printf(" Asymmetry has been calculated \n"); + }; + // + if ( asymm < 0. || ndf <= 0. || chi2 < 0. || tmax < 0. ){ + if ( debug ) printf(" Funny asymm||ndf||chi2||tmax values, fit failed \n"); + fitresult = 100; + }; + // + if ( draw ){ + // + tc->cd(); + // gStyle->SetOptStat(11111); + tc->SetTitle(); + th->SetTitle(""); + th->SetName(""); + th->SetMarkerStyle(20); + // axis titles + th->SetXTitle("Depth [X0]"); + th->SetYTitle("Energy [MIP]"); + // th->DrawCopy("Perror"); + th->DrawCopy(); + gh->SetMarkerSize(1.); + gh->SetMarkerStyle(21); + // gh->SetMarkerColor(kRed); + gh->Draw("Psame"); + lfit->Draw("same"); + tc->Modified(); + tc->Update(); + // + gStyle->SetLabelSize(0); + gStyle->SetNdivisions(1,"XY"); + // + } else { + if ( th ) th->Delete(); + }; + // +// gMinuit->SetPrintLevel(1); +// gMinuit->mnamin(); +// // gMinuit->mnfree(0); +// gMinuit->mncler(); +// gMinuit->mnrset(1); + // gMinuit->Delete(); + //ROOT::Math::Minimizer::Clear(); + // gMinuit->mnhelp("*"); + // + // delete lfit; + // +}; + +void CaloLong::Draw(){ + // + Process(); + Draw(-1); +}; + +void CaloLong::Draw(Int_t view){ + // + Int_t minv = 0; + Int_t maxv = 0; + // + if ( view == -1 ){ + maxv = -1; + } else { + minv = view; + maxv = view+1; + }; + // + Process(); + // + gStyle->SetLabelSize(0.04); + gStyle->SetNdivisions(510,"XY"); + // + if ( maxv != -1 ){ + for (Int_t v=minv; v(gDirectory->FindObject(hid)); + if ( tc ){ + // tc->Clear(); + } else { + tc = new TCanvas(hid,hid); + }; + // + TString thid = Form("hlongv%i%s",v,suf.Data()); + TH1F *th = dynamic_cast(gDirectory->FindObject(thid)); + if ( th ) th->Delete(); + // th->Clear(); + // th->Reset(); + // } else { + th = new TH1F(thid,thid,22,-0.5,21.5); + // }; + tc->cd(); + // + for (Int_t st=0;st<22;st++){ + th->Fill(st,eplane[v][st]); + }; + th->Draw(); + tc->Modified(); + tc->Update(); + }; + } else { + // + TString hid = Form("clongvyvx%s",suf.Data()); + TCanvas *tc = dynamic_cast(gDirectory->FindObject(hid)); + if ( tc ){ + } else { + tc = new TCanvas(hid,hid); + }; + // + TString thid = Form("hlongvyvx%s",suf.Data()); + TH1F *th = dynamic_cast(gDirectory->FindObject(thid)); + if ( th ) th->Delete(); + th = new TH1F(thid,thid,44,-0.5,43.5); + tc->cd(); + Int_t pp=0; + for (Int_t st=0;st<22;st++){ + for (Int_t v=1; v>=0;v--){ + // + th->Fill(pp,eplane[v][st]); + // + pp++; + }; + }; + th->Draw(); + tc->Modified(); + tc->Update(); + }; + // + gStyle->SetLabelSize(0); + gStyle->SetNdivisions(1,"XY"); + // +}; + +void CaloLong::Delete(){ + Clear(); + //delete this; +}; + + + + + + + + + + + + + + + + + + + +