#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; // } //////////////////////////////////////////////////////////////////////// /** * 1-dimension function describing lateral distribution of the * shower as viewed by calorimeter * (projection of 3d function in one direction) * * xi[0] = x or y coordinate relative to shower axis * parmin[0] = rt * parmin[1] = p * parmin[2] = rc * */ //////////////////////////////////////////////////////////////////////// Double_t cfradx(Double_t *xi, Double_t *parmin) { double fradxmin2,p,rt,rc,es,x,pig,norm,c; x=*xi; pig = acos(-1.); rt=parmin[0]; p=parmin[1]; rc=parmin[2]; norm=parmin[3]; c=parmin[4]; x=x-c; es=1.5; fradxmin2=p*pig*pow(rc,2)/pow((pow(x,2)+pow(rc,2)),es); fradxmin2=fradxmin2+(1-p)*pig*pow(rt,2)/pow((pow(x,2)+pow(rt,2)),es); fradxmin2=norm*fradxmin2/(2*pig); //cout<<"x,fradxmin2 "<< x<<" "<IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n"); // OBT = 0; PKT = 0; atime = 0; // debug = false; // }; Calo2D::Calo2D(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; // debug = false; // }; void CaloLat::Clear(){ // }; void Calo2D::Clear(){ // }; void CaloLat::Print(){ // Process(); // printf("==================== Calorimeter Lateral 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); printf("========================================================================\n"); // }; void Calo2D::Print(){ // Process(); // 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); printf("========================================================================\n"); // }; void CaloLat::Draw(){ // Process(); Draw(-1,-1); }; void Calo2D::Draw(){ // Process(); Draw(-1); }; void CaloLat::Draw(Int_t view,Int_t plane){ // Int_t minv = 0; Int_t maxv = 0; Int_t minp = 0; Int_t maxp = 0; // if ( view == -1 ){ minv = 0; maxv = 2; } else { minv = view; maxv = view+1; }; if ( plane == -1 ){ minp = 0; maxp = 22; } else { minp = plane; maxp = plane+1; }; // Process(); // gStyle->SetLabelSize(0.04); gStyle->SetNdivisions(510,"XY"); // for (Int_t v=minv; v(gDirectory->FindObject(hid)); if ( tc ){ // tc->Clear(); } else { tc = new TCanvas(hid,hid); }; // TString thid = Form("hlatv%ip%i",v,p); TH1F *th = dynamic_cast(gDirectory->FindObject(thid)); if ( th ) th->Delete(); // th->Clear(); // th->Reset(); // } else { th = new TH1F(thid,thid,96,-0.5,95.5); // }; tc->cd(); // for (Int_t st=0;st<96;st++){ th->Fill(st,estrip[v][p][st]); }; th->Draw(); tc->Modified(); tc->Update(); }; }; // gStyle->SetLabelSize(0); gStyle->SetNdivisions(1,"XY"); // }; void Calo2D::Draw(Int_t plane){ // Int_t minp = 0; Int_t maxp = 0; // if ( plane == -1 ){ minp = 0; maxp = 23; } else { minp = plane; maxp = plane+1; }; // Process(); // gStyle->SetLabelSize(0.04); gStyle->SetNdivisions(510,"XY"); // 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); gStyle->SetNdivisions(1,"XY"); // }; void CaloLat::Delete(){ Clear(); //delete this; }; void Calo2D::Delete(){ Clear(); //delete this; }; void CaloLat::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 // 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); estrip[view1][plane1][strip1] = mip1; }; // if ( debug ) this->Print(); if ( debug ) printf(" exit \n"); // }; void Calo2D::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 // 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 = ""; xyaverage = true; // }; 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++){ 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"); TH2F *th = dynamic_cast(gDirectory->FindObject(thid)); if ( th ) th->Delete(); // TString ghid = Form("glongfit"); 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.); numpo++; // th->Fill(xpos,enemip); if ( debug ) printf(" AVE Filling: st %i xpos %f energy %f \n",st,xpos,enemip); }; } 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.); // 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); TF1 *lfit = dynamic_cast(gDirectory->FindObject("lfit")); if ( lfit ) lfit->Delete(); 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 = "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...\n"); // // fitresult = th->Fit("lfit",optio); fitresult = gh->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"); 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(); }; // // 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; };