--- calo/flight/CaloProfile/src/CaloProfile.cpp 2007/11/09 09:11:24 1.1 +++ calo/flight/CaloProfile/src/CaloProfile.cpp 2009/08/12 14:57:00 1.10 @@ -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(); }; @@ -80,7 +285,7 @@ // }; -CaloLong::CaloLong(PamLevel2 *l2p){ +Calo2D::Calo2D(PamLevel2 *l2p){ // Clear(); // @@ -100,7 +305,7 @@ // }; -void CaloLong::Clear(){ +void Calo2D::Clear(){ // }; @@ -116,11 +321,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 +339,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; @@ -203,16 +407,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 +425,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",p); + 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 +482,12 @@ //delete this; }; -void CaloLong::Delete(){ +void Calo2D::Delete(){ Clear(); //delete this; }; + void CaloLat::Process(){ // if ( !L2 ){ @@ -335,7 +533,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 +562,823 @@ // // 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); + 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; + debug = false; + maskXE = false; + maskXO = false; + maskYE = false; + maskYO = false; + // + lmax = 0.; + umax = 100.; + slmax = ""; + sumax = ""; + // +}; + +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(){ + // 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); + 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++){ + 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]); + }; + }; + // + // 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"); // }; +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){ + /* SAMPLE OUTPUT: + Enter Infix Expression : A + B + C / (E - F) + Postfix Expression is : A B + C E F - / + + */ + if ( !s.Contains("t") ){ + 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 *ts= new TString(""); + ts->Prepend(s.Data()); + ts->ReplaceAll("t",1,g.Data(),g.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"); + 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"); + TH1F *th = dynamic_cast(gDirectory->FindObject(thid)); + if ( th ) th->Delete(); + Float_t xpos = 0.; + Float_t enemip = 0.; + Float_t xmax = NC * X0pl + 0.2; + // th = new TH1F(thid,thid,int(NC*1.5),-0.2,xmax); + th = new TH1F(thid,thid,100,-0.2,xmax); + // + // AGH, BUG! + // + Int_t mmin = 0; + Int_t mmax = 0; + if ( cont ){ + mmin = N; + mmax = NC+N; + } else { + mmin = 0; + mmax = NC; + }; + // + 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 > 0. ){ + th->Fill(xpos,enemip); + if ( debug ) printf(" 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]); +// }; +// }; +// }; +// // +// }; + }; + // + TF1 *lfit = new TF1("lfit",ccurve,0.,xmax,3); + 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; + if ( debug ){ + // optio = "MERBOV"; + // optio = "MEROV"; + // optio = "EROV"; + optio = "RNOV"; + if ( draw ) optio = "ROV"; + } else { + // optio = "MERNOQ"; + // optio = "ERNOQ"; + optio = "RNOQ"; + if ( draw ) optio = "ROQ"; + }; + // + if ( debug ) printf(" OK, start the fitting procedure...\n"); + // + fitresult = th->Fit("lfit",optio); + // + 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); + if ( sumax.MaybeRegexp() ) umax = this->Evaluate(sumax,tmax); + 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; + // 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"); + lfit->Draw("same"); + tc->Modified(); + tc->Update(); + // + gStyle->SetLabelSize(0); + gStyle->SetNdivisions(1,"XY"); + // + } else { + if ( th ) th->Delete(); + }; + // + 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",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(); + }; + // + gStyle->SetLabelSize(0); + gStyle->SetNdivisions(1,"XY"); + // +}; + +void CaloLong::Delete(){ + Clear(); + //delete this; +}; + + + + + + + + + + + + + + + + + + + +