/[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.2 by mocchiut, Mon Sep 22 20:08:25 2008 UTC revision 1.18 by mocchiut, Thu May 13 13:55:37 2010 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 80  CaloLat::CaloLat(PamLevel2 *l2p){   Line 281  CaloLat::CaloLat(PamLevel2 *l2p){  
281    PKT = 0;    PKT = 0;
282    atime = 0;    atime = 0;
283    //    //
284      suf = "";
285    debug = false;    debug = false;
286      usepl18x = false;
287    //    //
288  };  };
289    
# Line 96  Calo2D::Calo2D(PamLevel2 *l2p){   Line 299  Calo2D::Calo2D(PamLevel2 *l2p){  
299    PKT = 0;    PKT = 0;
300    atime = 0;    atime = 0;
301    //    //
302      suf = "";
303    debug = false;    debug = false;
304      usepl18x = false;
305    //    //
306  };  };
307    
# Line 174  void CaloLat::Draw(Int_t view,Int_t plan Line 379  void CaloLat::Draw(Int_t view,Int_t plan
379    //    //
380    for (Int_t v=minv; v<maxv;v++){    for (Int_t v=minv; v<maxv;v++){
381       for (Int_t p=minp; p<maxp;p++){       for (Int_t p=minp; p<maxp;p++){
382          TString hid = Form("clatv%ip%i",v,p);           TString hid = Form("clatv%ip%i%s",v,p,suf.Data());      
383          TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));          TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
384          if ( tc ){          if ( tc ){
385  //       tc->Clear();  //       tc->Clear();
# Line 182  void CaloLat::Draw(Int_t view,Int_t plan Line 387  void CaloLat::Draw(Int_t view,Int_t plan
387           tc = new TCanvas(hid,hid);           tc = new TCanvas(hid,hid);
388          };          };
389          //          //
390          TString thid = Form("hlatv%ip%i",v,p);            TString thid = Form("hlatv%ip%i%s",v,p,suf.Data());    
391          TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));          TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
392          if ( th ) th->Delete();          if ( th ) th->Delete();
393  //       th->Clear();  //       th->Clear();
# Line 225  void Calo2D::Draw(Int_t plane){ Line 430  void Calo2D::Draw(Int_t plane){
430    gStyle->SetNdivisions(510,"XY");    gStyle->SetNdivisions(510,"XY");
431    //    //
432    for (Int_t p=minp; p<maxp;p++){    for (Int_t p=minp; p<maxp;p++){
433      TString hid = Form("c2dp%i",p);          TString hid = Form("c2dp%i%s",p,suf.Data());        
434      TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));      TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
435      if ( tc ){      if ( tc ){
436        //         tc->Clear();        //         tc->Clear();
# Line 233  void Calo2D::Draw(Int_t plane){ Line 438  void Calo2D::Draw(Int_t plane){
438        tc = new TCanvas(hid,hid);        tc = new TCanvas(hid,hid);
439      };      };
440      //      //
441      TString thid = Form("h2dp%i",p);          TString thid = Form("h2dp%i%s",p,suf.Data());      
442      TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));      TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
443      if ( th ) th->Delete();      if ( th ) th->Delete();
444      //   th->Clear();      //   th->Clear();
# Line 324  void CaloLat::Process(){ Line 529  void CaloLat::Process(){
529    //    //
530    for (Int_t i=0; i<L2->GetCaloLevel1()->istrip ; i++){    for (Int_t i=0; i<L2->GetCaloLevel1()->istrip ; i++){
531      mip1 = L2->GetCaloLevel1()->DecodeEstrip(i,view1,plane1,strip1);      mip1 = L2->GetCaloLevel1()->DecodeEstrip(i,view1,plane1,strip1);
532        //
533        if ( !usepl18x && view1==0 && plane1==18 ) mip1 = 0.;
534        //
535      estrip[view1][plane1][strip1] = mip1;      estrip[view1][plane1][strip1] = mip1;
536    };    };
537    //    //
# Line 371  void Calo2D::Process(){ Line 579  void Calo2D::Process(){
579    //    //
580    for (Int_t i=0; i<L2->GetCaloLevel1()->istrip ; i++){    for (Int_t i=0; i<L2->GetCaloLevel1()->istrip ; i++){
581      mip1 = L2->GetCaloLevel1()->DecodeEstrip(i,view1,plane1,strip1);      mip1 = L2->GetCaloLevel1()->DecodeEstrip(i,view1,plane1,strip1);
582        //
583        if ( !usepl18x && view1==0 && plane1==18 ) mip1 = 0.;
584        //
585      es[view1][plane1][strip1] = mip1;      es[view1][plane1][strip1] = mip1;
586    };    };
587    //    //
# Line 451  CaloLong::CaloLong(PamLevel2 *l2p){   Line 662  CaloLong::CaloLong(PamLevel2 *l2p){  
662    PKT = 0;    PKT = 0;
663    atime = 0;    atime = 0;
664    //    //
665      clp = L2->GetCaloLevel2();
666      //
667    sel = true;    sel = true;
668    cont = false;    cont = false;
669    N = 0;    N = 0;
# Line 458  CaloLong::CaloLong(PamLevel2 *l2p){   Line 671  CaloLong::CaloLong(PamLevel2 *l2p){  
671    mask18b = -1;    mask18b = -1;
672    //    //
673    no18x = true;    no18x = true;
674      usepl18x = !no18x;
675    debug = false;    debug = false;
676      maskXE = false;
677      maskXO = false;
678      maskYE = false;
679      maskYO = false;
680      //
681      lmax = 0.;
682      umax = 100.;
683      slmax = "";
684      sumax = "";
685      suf = "";
686      xyaverage = true;
687      //
688      letmax = -1;
689      heavytail = false;
690      //  lmipth = 100.;
691      lmipth = 0.;
692    //    //
693  };  };
694    
695    TF1 *CaloLong::GetFit(){
696      TString fnam=Form("lfit%s",suf.Data());
697      TF1 *lfit  = dynamic_cast<TF1*>(gDirectory->FindObject(fnam));
698      if ( lfit ) return lfit;
699      return NULL;
700    }
701    
702    void CaloLong::MaskSection(TString sec){
703      sec.ToUpper();
704      if ( sec.Contains("XO") ) maskXO = true;
705      if ( sec.Contains("YO") ) maskYO = true;
706      if ( sec.Contains("XE") ) maskXE = true;
707      if ( sec.Contains("YE") ) maskYE = true;
708    }
709    
710    void CaloLong::UnMaskSections(){
711      this->UnMaskSection("XEXOYEYO");
712    }
713    
714    void CaloLong::UnMaskSection(TString sec){
715      sec.ToUpper();
716      if ( sec.Contains("XO") ) maskXO = false;
717      if ( sec.Contains("YO") ) maskYO = false;
718      if ( sec.Contains("XE") ) maskXE = false;
719      if ( sec.Contains("YE") ) maskYE = false;
720    }
721    
722  void CaloLong::Clear(){  void CaloLong::Clear(){
723    //    //
724      if ( debug ) printf(" Clear called \n");
725      //
726    memset(eplane,0, 2*22*sizeof(Float_t));    memset(eplane,0, 2*22*sizeof(Float_t));
727    //    //
728    chi2 = 0.;    chi2 = 0.;
729    ndf = 0.;    ndf = 0.;
730    E0 = 0.;    E0 = 0.;
731      defE0 = 0.;
732    a = 0.;    a = 0.;
733    b = 0.;    b = 0.;
734    errE0 = 0.;    errE0 = 0.;
# Line 484  void CaloLong::Clear(){ Line 744  void CaloLong::Clear(){
744    
745  void CaloLong::Print(){  void CaloLong::Print(){
746    //    //
747    Process();    Fit();
748    //    //
749    printf("==================== Calorimeter Longitudinal Profile =======================\n");    printf("==================== Calorimeter Longitudinal Profile =======================\n");
750    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 754  void CaloLong::Print(){
754    printf(" nchi2    :.. %f\n",chi2/ndf);    printf(" nchi2    :.. %f\n",chi2/ndf);
755    printf(" E0       :.. %f\n",E0);    printf(" E0       :.. %f\n",E0);
756    printf(" E0/260.  :.. %f\n",E0/260.);    printf(" E0/260.  :.. %f\n",E0/260.);
757      printf(" defE0    :.. %f\n",defE0);
758      printf(" lower    :.. %f\n",lmax);
759      printf(" upper    :.. %f\n",umax);
760      printf(" s lower  :.. %s\n",slmax.Data());
761      printf(" s upper  :.. %s\n",sumax.Data());
762    printf(" a        :.. %f\n",a);    printf(" a        :.. %f\n",a);
763    printf(" b        :.. %f\n",b);    printf(" b        :.. %f\n",b);
764    printf(" errE0    :.. %f\n",errE0);    printf(" errE0    :.. %f\n",errE0);
# Line 532  void CaloLong::SetNoWcalo(Int_t n){ Line 797  void CaloLong::SetNoWcalo(Int_t n){
797    };    };
798  }  }
799    
800    void CaloLong::SplitInto(Int_t NoWpreSampler, Int_t NoWcalo){
801      this->SetNoWpreSampler(0);
802      this->SetNoWcalo(0);
803      if ( NoWpreSampler < NoWcalo ){
804              this->SetNoWpreSampler(NoWpreSampler);
805              this->SetNoWcalo(NoWcalo);    
806      } else {
807              this->SetNoWcalo(NoWcalo);    
808              this->SetNoWpreSampler(NoWpreSampler);
809      };
810    }
811    
812  void CaloLong::Process(){  void CaloLong::Process(){
813    //      //  
814    if ( !L2 ){    if ( !L2 ){
# Line 585  void CaloLong::Process(){ Line 862  void CaloLong::Process(){
862    Int_t plane = 0;    Int_t plane = 0;
863    Int_t strip = 0;    Int_t strip = 0;
864    Float_t mip = 0.;    Float_t mip = 0.;
865      Bool_t gof = true;
866    for (Int_t i=0; i < L2->GetCaloLevel1()->istrip; i++){    for (Int_t i=0; i < L2->GetCaloLevel1()->istrip; i++){
867      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
868      eplane[view][plane] += mip;      //
869        if ( !usepl18x && view==0 && plane==18 ) mip = 0.;
870        //
871        gof = true;
872        if ( maskXE && (plane%2)==0 && view==1 ) gof = false;
873        if ( maskXO && (plane%2)!=0 && view==1 ) gof = false;
874        if ( maskYE && (plane%2)!=0 && view==0 ) gof = false;
875        if ( maskYO && (plane%2)==0 && view==0 ) gof = false;
876        if ( gof ) eplane[view][plane] += mip;
877    };    };
878    //    //
879    // inclination factor (steal from Daniele's code)    // inclination factor (stolen from Daniele's code)
880    //    //
881    Float_t ytgx = 0;    Float_t ytgx = 0.;
882    Float_t ytgy = 0;    Float_t ytgy = 0.;
883    ytgx = 0.76 * L2->GetCaloLevel2()->tanx[0];    ytgx = 0.76 * clp->tanx[0];
884    ytgy = 0.76 * L2->GetCaloLevel2()->tany[0];      ytgy = 0.76 * clp->tany[0];  
885    X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) );    X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) );
886    //    //
887    // Find experimental plane of maximum    // Find experimental plane of maximum
# Line 623  void CaloLong::Process(){ Line 909  void CaloLong::Process(){
909    //    //
910  };  };
911    
912    void CaloLong::SetEnergies(Float_t myene[][22]){
913      //  
914      if ( !L2 ){
915        printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
916        printf(" ERROR: CaloHough variables not filled \n");
917        return;
918      };
919      //
920      Bool_t newentry = false;
921      //
922      if ( L2->IsORB() ){
923        if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){
924          newentry = true;
925          OBT = L2->GetOrbitalInfo()->OBT;
926          PKT = L2->GetOrbitalInfo()->pkt_num;
927          atime = L2->GetOrbitalInfo()->absTime;
928        };
929      } else {
930        newentry = true;
931      };
932      //
933      if ( !newentry ) return;
934      //
935      if ( debug ) printf(" SET ENERGIES Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
936      //
937      Clear();
938      //
939      // let's start
940      //
941      if ( cont ){
942        for (Int_t i=0; i<22; i++){
943          if ( i == (18+N) ){
944            mask18b =  18 + N;
945            break;
946          };
947        };
948      };  
949      //  
950      if ( sel ){
951        for (Int_t i=0; i<22; i++){
952          if ( i == (18-N) ){
953            mask18b =  18 - N;
954            break;
955          };
956        };
957      };
958      //
959      //  if ( mask18b == 18 ) mask18b = -1;
960      //
961      Int_t view = 0;
962      Int_t plane = 0;
963      Bool_t gof = true;
964      for (view=0; view < 2; view++){
965        for (plane=0; plane < 22; plane++){      
966          gof = true;
967          if ( maskXE && (plane%2)==0 && view==1 ) gof = false;
968          if ( maskXO && (plane%2)!=0 && view==1 ) gof = false;
969          if ( maskYE && (plane%2)!=0 && view==0 ) gof = false;
970          if ( maskYO && (plane%2)==0 && view==0 ) gof = false;
971          if ( gof ) eplane[view][plane] = myene[view][plane];
972          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]);
973        };
974      };
975      //
976      if ( !clp && debug ) printf(" ERROR!! no CaloLevel2 obj!!\n");
977      //
978      // inclination factor (stolen from Daniele's code)
979      //
980      Float_t ytgx = 0.;
981      Float_t ytgy = 0.;
982      ytgx = 0.76 * clp->tanx[0];
983      ytgy = 0.76 * clp->tany[0];  
984      X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) );
985      if ( debug ) printf(" X0pl %f tanx %f tany %f \n",X0pl,ytgx,ytgy);
986      //
987      // Find experimental plane of maximum
988      //
989      Int_t pmax = 0;
990      Int_t vmax = 0;
991      Float_t emax = 0.;
992      for (Int_t v=0; v<2; v++){
993       for (Int_t i=0; i<22; i++){
994         if ( eplane[v][i] > emax ){
995           emax = eplane[v][i];
996           vmax = v;
997           pmax = i;
998         };
999       };
1000      };
1001      //
1002      //
1003      //
1004      if ( vmax == 0 ) pmax++;
1005      etmax = pmax * X0pl;
1006      letmax = etmax;
1007      //
1008      if ( debug ) this->Print();
1009      if ( debug ) printf(" exit \n");
1010      //
1011    };
1012    
1013  Double_t ccurve(Double_t *ti,Double_t *par){  Double_t ccurve(Double_t *ti,Double_t *par){
1014    //    //
# Line 642  void CaloLong::Fit(){ Line 1028  void CaloLong::Fit(){
1028    this->Fit(false);    this->Fit(false);
1029  };  };
1030    
1031    Float_t CaloLong::Evaluate(TString s, Float_t tmax, Float_t X0pl){
1032      /* SAMPLE OUTPUT:
1033         Enter Infix Expression : A + B + C / (E - F)
1034         Postfix Expression is : A B + C E F - / +  
1035      */  
1036      if ( !s.Contains("tmax") && !s.Contains("X0pl") ){
1037        printf(" ERROR, the input formula must contain \"t\"\n");
1038        return 0.;
1039      };
1040      if ( tmax != tmax ){
1041        printf(" ERROR, tmax is nan! \n");
1042        return 0.;
1043      };
1044      TString g=Form("%f",tmax);
1045      TString h=Form("%f",X0pl);
1046      TString *ts= new TString("");
1047      ts->Prepend(s.Data());
1048      ts->ReplaceAll("tmax",4,g.Data(),g.Capacity());
1049      ts->ReplaceAll("X0pl",4,h.Data(),h.Capacity());
1050      ts->Prepend("(");
1051      ts->Append(")");
1052      if ( debug )  printf(" ts %s tssize %i capac %i s %s g %s \n",ts->Data(),ts->Sizeof(),ts->Capacity(),s.Data(),g.Data());
1053      char in[50],post[50];  
1054      strcpy(&in[0],"                                                 ");
1055      strcpy(&in[0],ts->Data());
1056      strcpy(&post[0],"                                                 ");
1057      if ( debug ) printf("Infix Expression is : ##%s##\n",&in[0]);
1058      infix2postfix(&in[0],&post[0],1);
1059      if ( debug ) printf("Postfix Expression is : ##%s##\n",&post[0]);
1060      /* SAMPLE OUTPUT:
1061         Enter Postfix Expression : 3 5 + 2 /
1062         3 5 + 2 / EQUALS 4
1063      */
1064      Float_t res = evaluate(&post[0]);
1065      if ( debug ) printf("%s EQUALS %f\n",post,res);
1066      return res;
1067    }
1068    
1069  void CaloLong::Fit(Bool_t draw){  void CaloLong::Fit(Bool_t draw){
1070    //    //
1071    Process();    Process();
# Line 675  void CaloLong::Fit(Bool_t draw){ Line 1099  void CaloLong::Fit(Bool_t draw){
1099      gStyle->SetNdivisions(510,"XY");      gStyle->SetNdivisions(510,"XY");
1100    };    };
1101    //    //
1102    TString hid = Form("clongfit");          TString hid = Form("clongfit%s",suf.Data());  
1103    TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));    TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
1104    //  if ( tc ) tc->Delete();    //  if ( tc ) tc->Delete();
1105    //  if ( tc ) tc->Close();    //  if ( tc ) tc->Close();
# Line 685  void CaloLong::Fit(Bool_t draw){ Line 1109  void CaloLong::Fit(Bool_t draw){
1109      if ( tc ) tc->cd();      if ( tc ) tc->cd();
1110    };    };
1111    //    //
1112    TString thid = Form("hlongfit");          TString thid = Form("hlongfit%s",suf.Data());
1113    TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));    TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
1114    if ( th ) th->Delete();    if ( th ) th->Delete();
1115      //
1116      TString ghid = Form("glongfit%s",suf.Data());
1117      TGraphErrors *gh  = dynamic_cast<TGraphErrors*>(gDirectory->FindObject(ghid));
1118      if ( gh ) gh->Delete();
1119      //
1120    Float_t xpos = 0.;    Float_t xpos = 0.;
1121    Float_t enemip = 0.;    Float_t enemip = 0.;
1122    Float_t xmax = NC * X0pl + 0.2;    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);  
1123    //    //
1124    for (Int_t st=N;st<(N+NC);st++){    Double_t xxx[50];
1125      Double_t exx[50];
1126      Double_t yyy[50];
1127      Double_t eyy[50];
1128      Int_t numpo = 0;
1129      memset(xxx,0,50*sizeof(Double_t));
1130      memset(exx,0,50*sizeof(Double_t));
1131      memset(yyy,0,50*sizeof(Double_t));
1132      memset(eyy,0,50*sizeof(Double_t));
1133      //
1134      // AGH, BUG!
1135      //
1136      Int_t mmin = 0;
1137      Int_t mmax = 0;
1138      if ( cont ){
1139        mmin = N;
1140        mmax = NC+N;
1141      } else {
1142        mmin = 0;
1143        mmax = NC;
1144      };
1145      //
1146      Float_t emax = 0.;
1147      Float_t qtotparz = 0.;
1148      for (Int_t st=mmin;st<mmax+1;st++){
1149      enemip = 0.;      enemip = 0.;
1150      xpos = (st - N) * X0pl;      xpos = (st - mmin) * X0pl;
1151      if ( st > N && st < (N+NC-1) ){            //
1152        if ( no18x && ( st == 18+1 || st == mask18b+1 )){      if ( xyaverage ){
1153          enemip = 2. * eplane[1][st];        //
1154          if ( st > mmin && st < mmax ){      
1155            if ( no18x && ( st == 18+1 || st == mask18b+1 )){
1156              if ( !maskYO ){
1157                enemip = 2. * eplane[1][st];
1158              } else {
1159                enemip = eplane[1][st];
1160              };
1161            } else {
1162              enemip = eplane[0][st-1] + eplane[1][st];
1163            };
1164        } else {        } else {
1165          enemip = eplane[0][st-1] + eplane[1][st];          if ( st == mmin ){
1166              if ( !maskYE ){
1167                enemip = 2. * eplane[1][st];
1168              } else {
1169                enemip = eplane[1][st];
1170              };
1171            };
1172            if ( st == mmax ){
1173              if ( !maskXE ){
1174                enemip = 2. * eplane[0][st-1];
1175              } else {
1176                enemip = eplane[0][st-1];
1177              };
1178            };
1179          };
1180          //
1181          qtotparz += enemip;
1182          if ( enemip > emax ) emax = enemip;
1183          if ( enemip > 0. ){
1184            xxx[numpo] = xpos;
1185            exx[numpo] = 0.1;
1186            yyy[numpo] = enemip;
1187            eyy[numpo] = sqrt(enemip*3.)+sqrt(5.);
1188            //if ( xpos > letmax && enemip > lmipth && heavytail) eyy[numpo] = (sqrt(enemip*3.)+sqrt(5.))/numpo;
1189            if ( xpos > letmax && enemip > lmipth && heavytail) eyy[numpo] = sqrt(enemip)/5.;
1190            if ( xpos > letmax-1 && xpos < letmax+1 && heavytail ) eyy[numpo] /= 5.;
1191            //if ( xpos > letmax-2 && xpos < letmax+1 && heavytail ) eyy[numpo] /= 7.;
1192            if ( xpos < 3. && heavytail ) eyy[numpo] /= 5.;
1193            numpo++;
1194            //      th->Fill(xpos,enemip);
1195            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);
1196        };        };
1197      } else {      } else {
1198        if ( st == N ) enemip = 2. * eplane[1][st];        for (Int_t jj=0; jj<2; jj++){
1199        if ( st == (N+NC-1) ) enemip = 2. * eplane[0][st];          if ( st > mmin && st < mmax ){
1200      };            if ( jj == 0 && no18x && ( st == 18+1 || st == mask18b+1 )){
1201      //              enemip = eplane[1][st];
1202      if ( enemip > 0. ){            } else {
1203        th->Fill(xpos,enemip);              if ( jj == 0 ){
1204        if ( debug ) printf(" Filling: st %i xpos %f energy %f \n",st,xpos,enemip);                enemip = eplane[jj][st-1];
1205                } else {
1206                  enemip = eplane[jj][st];
1207                };        
1208              };
1209            } else {
1210              if ( st == mmin && jj == 1 ){
1211                enemip = eplane[jj][st];
1212              };
1213              if ( st == mmax && jj == 0){
1214                enemip = eplane[jj][st-1];
1215              };
1216            };
1217            //
1218            qtotparz += enemip;
1219          if ( enemip > emax ) emax = enemip;
1220            if ( enemip > 0. ){
1221              xxx[numpo] = xpos;
1222              exx[numpo] = 0.1;
1223              yyy[numpo] = enemip;
1224              eyy[numpo] = sqrt(enemip*3.)+sqrt(5.);
1225              //if ( xpos > letmax && enemip > lmipth && heavytail) eyy[numpo] = (sqrt(enemip*3.)+sqrt(5.))/numpo;
1226              //if ( xpos > letmax && enemip > lmipth && heavytail ) eyy[numpo] = (sqrt(enemip*3.)+sqrt(5.))/numpo;
1227              if ( xpos > letmax && enemip > lmipth && heavytail) eyy[numpo] = sqrt(enemip)/5.;
1228              if ( xpos > letmax-1 && xpos < letmax+1 && heavytail ) eyy[numpo] /= 5.;
1229            //if ( xpos > letmax-2 && xpos < letmax+1 && heavytail ) eyy[numpo] /= 7.;
1230              if ( xpos < 3. && heavytail ) eyy[numpo] /= 5.;        
1231              //      eyy[numpo] = sqrt(enemip)/(st*0.95);
1232              numpo++;
1233              //      th->Fill(xpos,enemip);
1234              if ( debug ) printf(" XY Filling: st %i xpos %f energy %f \n",st,xpos,enemip);
1235            };
1236          };        
1237      };      };
1238    
1239      //      //
1240      //    for (Int_t v=1; v>=0;v--)// {      //    for (Int_t v=1; v>=0;v--)// {
1241  //       //  //       //
# Line 737  void CaloLong::Fit(Bool_t draw){ Line 1261  void CaloLong::Fit(Bool_t draw){
1261  //     };  //     };
1262    };    };
1263    //    //
1264    TF1 *lfit = new TF1("lfit",ccurve,0.,xmax,3);    //  th = new TH2F(thid,thid,int(NC*1.5),-0.2,xmax);
1265    E0 = L2->GetCaloLevel2()->qtot;    th = new TH2F(thid,thid,1000,-0.2,xmax,1000,0.,emax*1.2);
1266      gh = new TGraphErrors(numpo,xxx,yyy,exx,eyy);
1267      TString fnam=Form("lfit%s",suf.Data());
1268      TF1 *lfit  = dynamic_cast<TF1*>(gDirectory->FindObject(fnam));
1269      if ( lfit ) lfit->Delete();
1270      lfit = new TF1(fnam,ccurve,0.,xmax,3);
1271    //   if ( !lfit ){
1272    //     lfit = new TF1(fnam,ccurve,0.,xmax,3);
1273    //   };
1274    //   lfit->Clear();
1275    //   lfit->SetName(fnam);
1276    //   lfit->SetRange(0.,xmax);
1277    //  lfit->Print();
1278      //
1279      if ( debug ) printf("qtot %f qtotparz %f \n",clp->qtot,qtotparz);
1280      E0 = qtotparz;
1281      //  E0 = clp->qtot;
1282    a = 5.;    a = 5.;
1283    b = 0.5;    b = 0.5;
1284    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 746  void CaloLong::Fit(Bool_t draw){ Line 1286  void CaloLong::Fit(Bool_t draw){
1286    //  lfit->SetParLimits(0,0.,1000.);    //  lfit->SetParLimits(0,0.,1000.);
1287    //  lfit->SetParLimits(1,-1.,80.);    //  lfit->SetParLimits(1,-1.,80.);
1288    //  lfit->SetParLimits(2,-1.,10.);    //  lfit->SetParLimits(2,-1.,10.);
1289    TString optio;    // TString optio = "ROW"; // "RO"
1290      TString optio = "RO"; // "RO"
1291    if ( debug ){    if ( debug ){
1292      //    optio = "MERBOV";      optio += "NV";
1293      //    optio = "MEROV";      if ( draw ) optio += "V";
     //    optio = "EROV";  
     optio = "RNOV";  
     if ( draw ) optio = "ROV";  
1294    } else {    } else {
1295      //    optio = "MERNOQ";      optio += "NQ";
1296      //    optio = "ERNOQ";      if ( draw ) optio += "Q";
     optio = "RNOQ";  
     if ( draw ) optio = "ROQ";  
1297    };    };
1298    //    //
1299    if ( debug ) printf(" OK, start the fitting procedure...\n");    if ( debug ) printf(" OK, start the fitting procedure... options: %s \n",optio.Data());
1300    //    //
1301    fitresult = th->Fit("lfit",optio);    //  fitresult = th->Fit("lfit",optio);
1302      //  lfit->Update();
1303      fitresult = gh->Fit("lfit",optio);
1304      //  fitresult = gh->Fit("lfit","ROQW");
1305    //    //
1306    if ( debug ) printf(" the fit is done! result: %i \n",fitresult);    if ( debug ) printf(" the fit is done! result: %i \n",fitresult);
1307    //    //
# Line 781  void CaloLong::Fit(Bool_t draw){ Line 1320  void CaloLong::Fit(Bool_t draw){
1320    if ( fitresult != 0 ){    if ( fitresult != 0 ){
1321      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");
1322      asymm = -1.;      asymm = -1.;
1323        defE0 = -1.;
1324    } else {    } else {
1325        if ( slmax.MaybeRegexp() ) lmax = this->Evaluate(slmax,tmax,X0pl);
1326        if ( sumax.MaybeRegexp() ) umax = this->Evaluate(sumax,tmax,X0pl);
1327      Int_t npp = 1000;      Int_t npp = 1000;
1328        double *xpp=new double[npp];
1329        double *wpp=new double[npp];
1330        lfit->CalcGaussLegendreSamplingPoints(npp,xpp,wpp,1e-12);
1331        defE0 = lfit->IntegralFast(npp,xpp,wpp,lmax,umax);
1332        //
1333      double *xp=new double[npp];      double *xp=new double[npp];
1334      double *wp=new double[npp];      double *wp=new double[npp];
1335      lfit->CalcGaussLegendreSamplingPoints(npp,xp,wp,1e-12);      lfit->CalcGaussLegendreSamplingPoints(npp,xp,wp,1e-12);
# Line 798  void CaloLong::Fit(Bool_t draw){ Line 1345  void CaloLong::Fit(Bool_t draw){
1345      delete w;      delete w;
1346      delete xp;      delete xp;
1347      delete wp;      delete wp;
1348        delete xpp;
1349        delete wpp;
1350      //    Float_t i10max = lfit->Integral(0.,10.*tmax);      //    Float_t i10max = lfit->Integral(0.,10.*tmax);
1351      if ( debug ) printf(" Integral: %f \n",i10max);      if ( debug ) printf(" Integral: %f \n",i10max);
1352      //      //
# Line 807  void CaloLong::Fit(Bool_t draw){ Line 1356  void CaloLong::Fit(Bool_t draw){
1356        if ( debug ) printf(" i10max == imax, asymm undefined\n");        if ( debug ) printf(" i10max == imax, asymm undefined\n");
1357        asymm = -2.;        asymm = -2.;
1358      };      };
1359        if ( asymm != asymm ){
1360          if ( debug ) printf(" asymm is nan \n");      
1361          asymm = -3.;
1362        };
1363      //lfit->Integral(0.,tmax)/(lfit->Integral(0.,10.*tmax)-lfit->Integral(0.,tmax));      //lfit->Integral(0.,tmax)/(lfit->Integral(0.,10.*tmax)-lfit->Integral(0.,tmax));
1364      if ( debug ) printf(" Asymmetry has been calculated \n");      if ( debug ) printf(" Asymmetry has been calculated \n");
1365    };    };
1366    //    //
1367      if ( asymm < 0. || ndf <= 0. || chi2 < 0. || tmax < 0. ){
1368        if ( debug ) printf(" Funny asymm||ndf||chi2||tmax values, fit failed \n");
1369        fitresult = 100;
1370      };
1371      //
1372    if ( draw ){    if ( draw ){
1373      //      //
1374      tc->cd();          tc->cd();    
# Line 822  void CaloLong::Fit(Bool_t draw){ Line 1380  void CaloLong::Fit(Bool_t draw){
1380      // axis titles      // axis titles
1381      th->SetXTitle("Depth [X0]");      th->SetXTitle("Depth [X0]");
1382      th->SetYTitle("Energy [MIP]");      th->SetYTitle("Energy [MIP]");
1383      th->DrawCopy("Perror");      //    th->DrawCopy("Perror");
1384        th->DrawCopy();
1385        gh->SetMarkerSize(1.);
1386        gh->SetMarkerStyle(21);
1387        //    gh->SetMarkerColor(kRed);
1388        gh->Draw("Psame");
1389      lfit->Draw("same");      lfit->Draw("same");
1390      tc->Modified();      tc->Modified();
1391      tc->Update();      tc->Update();
# Line 830  void CaloLong::Fit(Bool_t draw){ Line 1393  void CaloLong::Fit(Bool_t draw){
1393      gStyle->SetLabelSize(0);      gStyle->SetLabelSize(0);
1394      gStyle->SetNdivisions(1,"XY");      gStyle->SetNdivisions(1,"XY");
1395      //      //
1396      } else {
1397        if ( th ) th->Delete();
1398    };    };
1399    //    //
1400    delete lfit;  //   gMinuit->SetPrintLevel(1);
1401    //   gMinuit->mnamin();
1402    //   //  gMinuit->mnfree(0);
1403    //   gMinuit->mncler();
1404    //   gMinuit->mnrset(1);
1405      //  gMinuit->Delete();
1406      //ROOT::Math::Minimizer::Clear();
1407      //  gMinuit->mnhelp("*");
1408      //
1409      //  delete lfit;
1410    //    //
1411  };  };
1412    
# Line 861  void CaloLong::Draw(Int_t view){ Line 1435  void CaloLong::Draw(Int_t view){
1435    //    //
1436    if ( maxv != -1 ){    if ( maxv != -1 ){
1437      for (Int_t v=minv; v<maxv;v++){      for (Int_t v=minv; v<maxv;v++){
1438        TString hid = Form("clongv%i",v);        TString hid = Form("clongv%i%s",v,suf.Data());    
1439        TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));        TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
1440        if ( tc ){        if ( tc ){
1441          //       tc->Clear();          //       tc->Clear();
# Line 869  void CaloLong::Draw(Int_t view){ Line 1443  void CaloLong::Draw(Int_t view){
1443          tc = new TCanvas(hid,hid);          tc = new TCanvas(hid,hid);
1444        };        };
1445        //        //
1446        TString thid = Form("hlongv%i",v);                TString thid = Form("hlongv%i%s",v,suf.Data());  
1447        TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));        TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
1448        if ( th ) th->Delete();        if ( th ) th->Delete();
1449        //         th->Clear();        //         th->Clear();
# Line 888  void CaloLong::Draw(Int_t view){ Line 1462  void CaloLong::Draw(Int_t view){
1462      };      };
1463    } else {    } else {
1464      //      //
1465      TString hid = Form("clongvyvx");          TString hid = Form("clongvyvx%s",suf.Data());      
1466      TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));      TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
1467      if ( tc ){      if ( tc ){
1468      } else {      } else {
1469        tc = new TCanvas(hid,hid);        tc = new TCanvas(hid,hid);
1470      };      };
1471      //      //
1472      TString thid = Form("hlongvyvx");        TString thid = Form("hlongvyvx%s",suf.Data());      
1473      TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));      TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
1474      if ( th ) th->Delete();      if ( th ) th->Delete();
1475      th = new TH1F(thid,thid,44,-0.5,43.5);      th = new TH1F(thid,thid,44,-0.5,43.5);
# Line 923  void CaloLong::Delete(){ Line 1497  void CaloLong::Delete(){
1497    Clear();    Clear();
1498    //delete this;    //delete this;
1499  };  };
1500    
1501    
1502    
1503    
1504    
1505    
1506    
1507    
1508    
1509    
1510    
1511    
1512    
1513    
1514    
1515    
1516    
1517    
1518    
1519    
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.18
  ViewVC Help
Powered by ViewVC 1.1.23