/[PAMELA software]/calo/flight/CaloProfile/src/CaloProfile.cpp
ViewVC logotype

Diff of /calo/flight/CaloProfile/src/CaloProfile.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.4 by mocchiut, Wed Jun 10 13:00:23 2009 UTC revision 1.10 by mocchiut, Wed Aug 12 14:57:00 2009 UTC
# Line 4  ClassImp(CaloLat); Line 4  ClassImp(CaloLat);
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
# Line 451  CaloLong::CaloLong(PamLevel2 *l2p){   Line 652  CaloLong::CaloLong(PamLevel2 *l2p){  
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;
# Line 459  CaloLong::CaloLong(PamLevel2 *l2p){   Line 662  CaloLong::CaloLong(PamLevel2 *l2p){  
662    //    //
663    no18x = true;    no18x = true;
664    debug = false;    debug = false;
665      maskXE = false;
666      maskXO = false;
667      maskYE = false;
668      maskYO = false;
669      //
670      lmax = 0.;
671      umax = 100.;
672      slmax = "";
673      sumax = "";
674    //    //
675  };  };
676    
677    void CaloLong::MaskSection(TString sec){
678      sec.ToUpper();
679      if ( sec.Contains("XO") ) maskXO = true;
680      if ( sec.Contains("YO") ) maskYO = true;
681      if ( sec.Contains("XE") ) maskXE = true;
682      if ( sec.Contains("YE") ) maskYE = true;
683    }
684    
685    void CaloLong::UnMaskSections(){
686      this->UnMaskSection("XEXOYEYO");
687    }
688    
689    void CaloLong::UnMaskSection(TString sec){
690      sec.ToUpper();
691      if ( sec.Contains("XO") ) maskXO = false;
692      if ( sec.Contains("YO") ) maskYO = false;
693      if ( sec.Contains("XE") ) maskXE = false;
694      if ( sec.Contains("YE") ) maskYE = false;
695    }
696    
697  void CaloLong::Clear(){  void CaloLong::Clear(){
698    //    //
699    memset(eplane,0, 2*22*sizeof(Float_t));    memset(eplane,0, 2*22*sizeof(Float_t));
# Line 469  void CaloLong::Clear(){ Line 701  void CaloLong::Clear(){
701    chi2 = 0.;    chi2 = 0.;
702    ndf = 0.;    ndf = 0.;
703    E0 = 0.;    E0 = 0.;
704      defE0 = 0.;
705    a = 0.;    a = 0.;
706    b = 0.;    b = 0.;
707    errE0 = 0.;    errE0 = 0.;
# Line 484  void CaloLong::Clear(){ Line 717  void CaloLong::Clear(){
717    
718  void CaloLong::Print(){  void CaloLong::Print(){
719    //    //
720    Process();    Fit();
721    //    //
722    printf("==================== Calorimeter Longitudinal Profile =======================\n");    printf("==================== Calorimeter Longitudinal Profile =======================\n");
723    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
# Line 494  void CaloLong::Print(){ Line 727  void CaloLong::Print(){
727    printf(" nchi2    :.. %f\n",chi2/ndf);    printf(" nchi2    :.. %f\n",chi2/ndf);
728    printf(" E0       :.. %f\n",E0);    printf(" E0       :.. %f\n",E0);
729    printf(" E0/260.  :.. %f\n",E0/260.);    printf(" E0/260.  :.. %f\n",E0/260.);
730      printf(" defE0    :.. %f\n",defE0);
731      printf(" lower    :.. %f\n",lmax);
732      printf(" upper    :.. %f\n",umax);
733      printf(" s lower  :.. %s\n",slmax.Data());
734      printf(" s upper  :.. %s\n",sumax.Data());
735    printf(" a        :.. %f\n",a);    printf(" a        :.. %f\n",a);
736    printf(" b        :.. %f\n",b);    printf(" b        :.. %f\n",b);
737    printf(" errE0    :.. %f\n",errE0);    printf(" errE0    :.. %f\n",errE0);
# Line 597  void CaloLong::Process(){ Line 835  void CaloLong::Process(){
835    Int_t plane = 0;    Int_t plane = 0;
836    Int_t strip = 0;    Int_t strip = 0;
837    Float_t mip = 0.;    Float_t mip = 0.;
838      Bool_t gof = true;
839    for (Int_t i=0; i < L2->GetCaloLevel1()->istrip; i++){    for (Int_t i=0; i < L2->GetCaloLevel1()->istrip; i++){
840      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
841      eplane[view][plane] += mip;      gof = true;
842        if ( maskXE && (plane%2)==0 && view==1 ) gof = false;
843        if ( maskXO && (plane%2)!=0 && view==1 ) gof = false;
844        if ( maskYE && (plane%2)!=0 && view==0 ) gof = false;
845        if ( maskYO && (plane%2)==0 && view==0 ) gof = false;
846        if ( gof ) eplane[view][plane] += mip;
847    };    };
848    //    //
849    // inclination factor (stolen from Daniele's code)    // inclination factor (stolen from Daniele's code)
850    //    //
851    Float_t ytgx = 0;    Float_t ytgx = 0;
852    Float_t ytgy = 0;    Float_t ytgy = 0;
853    ytgx = 0.76 * L2->GetCaloLevel2()->tanx[0];    ytgx = 0.76 * clp->tanx[0];
854    ytgy = 0.76 * L2->GetCaloLevel2()->tany[0];      ytgy = 0.76 * clp->tany[0];  
855    X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) );    X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) );
856    //    //
857    // Find experimental plane of maximum    // Find experimental plane of maximum
# Line 635  void CaloLong::Process(){ Line 879  void CaloLong::Process(){
879    //    //
880  };  };
881    
882    void CaloLong::SetEnergies(Float_t myene[][22]){
883      //  
884      if ( !L2 ){
885        printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
886        printf(" ERROR: CaloHough variables not filled \n");
887        return;
888      };
889      //
890      Bool_t newentry = false;
891      //
892      if ( L2->IsORB() ){
893        if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){
894          newentry = true;
895          OBT = L2->GetOrbitalInfo()->OBT;
896          PKT = L2->GetOrbitalInfo()->pkt_num;
897          atime = L2->GetOrbitalInfo()->absTime;
898        };
899      } else {
900        newentry = true;
901      };
902      //
903      if ( !newentry ) return;
904      //
905      if ( debug ) printf(" SET ENERGIES Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
906      //
907      Clear();
908      //
909      // let's start
910      //
911      if ( cont ){
912        for (Int_t i=0; i<22; i++){
913          if ( i == (18+N) ){
914            mask18b =  18 + N;
915            break;
916          };
917        };
918      };  
919      //  
920      if ( sel ){
921        for (Int_t i=0; i<22; i++){
922          if ( i == (18-N) ){
923            mask18b =  18 - N;
924            break;
925          };
926        };
927      };
928      //
929      //  if ( mask18b == 18 ) mask18b = -1;
930      //
931      Int_t view = 0;
932      Int_t plane = 0;
933      Bool_t gof = true;
934      for (view=0; view < 2; view++){
935        for (plane=0; plane < 22; plane++){      
936          gof = true;
937          if ( maskXE && (plane%2)==0 && view==1 ) gof = false;
938          if ( maskXO && (plane%2)!=0 && view==1 ) gof = false;
939          if ( maskYE && (plane%2)!=0 && view==0 ) gof = false;
940          if ( maskYO && (plane%2)==0 && view==0 ) gof = false;
941          if ( gof ) eplane[view][plane] = myene[view][plane];
942          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]);
943        };
944      };
945      //
946      // inclination factor (stolen from Daniele's code)
947      //
948      Float_t ytgx = 0;
949      Float_t ytgy = 0;
950      ytgx = 0.76 * clp->tanx[0];
951      ytgy = 0.76 * clp->tany[0];  
952      X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) );
953      //
954      // Find experimental plane of maximum
955      //
956      Int_t pmax = 0;
957      Int_t vmax = 0;
958      Float_t emax = 0.;
959      for (Int_t v=0; v<2; v++){
960       for (Int_t i=0; i<22; i++){
961         if ( eplane[v][i] > emax ){
962           emax = eplane[v][i];
963           vmax = v;
964           pmax = i;
965         };
966       };
967      };
968      //
969      //
970      //
971      if ( vmax == 0 ) pmax++;
972      etmax = pmax * X0pl;
973      //
974      if ( debug ) this->Print();
975      if ( debug ) printf(" exit \n");
976      //
977    };
978    
979  Double_t ccurve(Double_t *ti,Double_t *par){  Double_t ccurve(Double_t *ti,Double_t *par){
980    //    //
# Line 654  void CaloLong::Fit(){ Line 994  void CaloLong::Fit(){
994    this->Fit(false);    this->Fit(false);
995  };  };
996    
997    Float_t CaloLong::Evaluate(TString s, Float_t tmax){
998      /* SAMPLE OUTPUT:
999         Enter Infix Expression : A + B + C / (E - F)
1000         Postfix Expression is : A B + C E F - / +  
1001      */  
1002      if ( !s.Contains("t") ){
1003        printf(" ERROR, the input formula must contain \"t\"\n");
1004        return 0.;
1005      };
1006      if ( tmax != tmax ){
1007        printf(" ERROR, tmax is nan! \n");
1008        return 0.;
1009      };
1010      TString g=Form("%f",tmax);
1011      TString *ts= new TString("");
1012      ts->Prepend(s.Data());
1013      ts->ReplaceAll("t",1,g.Data(),g.Capacity());
1014      ts->Prepend("(");
1015      ts->Append(")");
1016      if ( debug )  printf(" ts %s tssize %i capac %i s %s g %s \n",ts->Data(),ts->Sizeof(),ts->Capacity(),s.Data(),g.Data());
1017      char in[50],post[50];  
1018      strcpy(&in[0],"                                                 ");
1019      strcpy(&in[0],ts->Data());
1020      strcpy(&post[0],"                                                 ");
1021      if ( debug ) printf("Infix Expression is : ##%s##\n",&in[0]);
1022      infix2postfix(&in[0],&post[0],1);
1023      if ( debug ) printf("Postfix Expression is : ##%s##\n",&post[0]);
1024      /* SAMPLE OUTPUT:
1025         Enter Postfix Expression : 3 5 + 2 /
1026         3 5 + 2 / EQUALS 4
1027      */
1028      Float_t res = evaluate(&post[0]);
1029      if ( debug ) printf("%s EQUALS %f\n",post,res);
1030      return res;
1031    }
1032    
1033  void CaloLong::Fit(Bool_t draw){  void CaloLong::Fit(Bool_t draw){
1034    //    //
1035    Process();    Process();
# Line 724  void CaloLong::Fit(Bool_t draw){ Line 1100  void CaloLong::Fit(Bool_t draw){
1100      xpos = (st - mmin) * X0pl;      xpos = (st - mmin) * X0pl;
1101      if ( st > mmin && st < mmax ){            if ( st > mmin && st < mmax ){      
1102        if ( no18x && ( st == 18+1 || st == mask18b+1 )){        if ( no18x && ( st == 18+1 || st == mask18b+1 )){
1103          enemip = 2. * eplane[1][st];          if ( !maskYO ){
1104              enemip = 2. * eplane[1][st];
1105            } else {
1106              enemip = eplane[1][st];
1107            };
1108        } else {        } else {
1109          enemip = eplane[0][st-1] + eplane[1][st];          enemip = eplane[0][st-1] + eplane[1][st];
1110        };        };
1111      } else {      } else {
1112        if ( st == mmin ) enemip = 2. * eplane[1][st];        if ( st == mmin ){
1113        if ( st == mmax ) enemip = 2. * eplane[0][st-1];          if ( !maskYE ){
1114              enemip = 2. * eplane[1][st];
1115            } else {
1116              enemip = eplane[1][st];
1117            };
1118          };
1119          if ( st == mmax ){
1120            if ( !maskXE ){
1121              enemip = 2. * eplane[0][st-1];
1122            } else {
1123              enemip = eplane[0][st-1];
1124            };
1125          };
1126      };      };
1127      //      //
1128      qtotparz += enemip;      qtotparz += enemip;
# Line 764  void CaloLong::Fit(Bool_t draw){ Line 1156  void CaloLong::Fit(Bool_t draw){
1156    };    };
1157    //    //
1158    TF1 *lfit = new TF1("lfit",ccurve,0.,xmax,3);    TF1 *lfit = new TF1("lfit",ccurve,0.,xmax,3);
1159    if ( debug ) printf("qtot %f qtotparz %f \n",L2->GetCaloLevel2()->qtot,qtotparz);    if ( debug ) printf("qtot %f qtotparz %f \n",clp->qtot,qtotparz);
1160    E0 = qtotparz;    E0 = qtotparz;
1161    //  E0 = L2->GetCaloLevel2()->qtot;    //  E0 = clp->qtot;
1162    a = 5.;    a = 5.;
1163    b = 0.5;    b = 0.5;
1164    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);
# Line 809  void CaloLong::Fit(Bool_t draw){ Line 1201  void CaloLong::Fit(Bool_t draw){
1201    if ( fitresult != 0 ){    if ( fitresult != 0 ){
1202      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");
1203      asymm = -1.;      asymm = -1.;
1204        defE0 = -1.;
1205    } else {    } else {
1206        if ( slmax.MaybeRegexp() ) lmax = this->Evaluate(slmax,tmax);
1207        if ( sumax.MaybeRegexp() ) umax = this->Evaluate(sumax,tmax);
1208      Int_t npp = 1000;      Int_t npp = 1000;
1209        double *xpp=new double[npp];
1210        double *wpp=new double[npp];
1211        lfit->CalcGaussLegendreSamplingPoints(npp,xpp,wpp,1e-12);
1212        defE0 = lfit->IntegralFast(npp,xpp,wpp,lmax,umax);
1213        //
1214      double *xp=new double[npp];      double *xp=new double[npp];
1215      double *wp=new double[npp];      double *wp=new double[npp];
1216      lfit->CalcGaussLegendreSamplingPoints(npp,xp,wp,1e-12);      lfit->CalcGaussLegendreSamplingPoints(npp,xp,wp,1e-12);
# Line 962  void CaloLong::Delete(){ Line 1362  void CaloLong::Delete(){
1362    Clear();    Clear();
1363    //delete this;    //delete this;
1364  };  };
1365    
1366    
1367    
1368    
1369    
1370    
1371    
1372    
1373    
1374    
1375    
1376    
1377    
1378    
1379    
1380    
1381    
1382    
1383    
1384    
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.23