/[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.1 by mocchiut, Fri Nov 9 09:11:24 2007 UTC revision 1.11 by mocchiut, Thu Aug 13 15:56:39 2009 UTC
# Line 1  Line 1 
1  #include <CaloProfile.h>  #include <CaloProfile.h>
2    //
3    ClassImp(CaloLat);
4    ClassImp(CaloLong);
5    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  /**  /**
# Line 60  CaloLat::CaloLat(){ Line 265  CaloLat::CaloLat(){
265  /**  /**
266   * Default constructor   * Default constructor
267   */   */
268  CaloLong::CaloLong(){  Calo2D::Calo2D(){
269    Clear();    Clear();
270  };  };
271    
# Line 80  CaloLat::CaloLat(PamLevel2 *l2p){   Line 285  CaloLat::CaloLat(PamLevel2 *l2p){  
285    //    //
286  };  };
287    
288  CaloLong::CaloLong(PamLevel2 *l2p){    Calo2D::Calo2D(PamLevel2 *l2p){  
289    //    //
290    Clear();    Clear();
291    //    //
# Line 100  void CaloLat::Clear(){ Line 305  void CaloLat::Clear(){
305    //    //
306  };  };
307    
308  void CaloLong::Clear(){  void Calo2D::Clear(){
309    //    //
310  };  };
311    
# Line 116  void CaloLat::Print(){ Line 321  void CaloLat::Print(){
321    //    //
322  };  };
323    
324  void CaloLong::Print(){  void Calo2D::Print(){
325    //    //
326    Process();    Process();
327    //    //
328    printf("==================== Calorimeter Lateral Profile =======================\n");    printf("==================== Calorimeter 2D Profile =======================\n");
329    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
330  //  printf(" nx [number of X combination]:.. %i\n",nx);  //  printf(" nx [number of X combination]:.. %i\n",nx);
331  //  printf(" ny [number of Y combination]:.. %i\n",ny);  //  printf(" ny [number of Y combination]:.. %i\n",ny);
# Line 134  void CaloLat::Draw(){ Line 339  void CaloLat::Draw(){
339    Draw(-1,-1);    Draw(-1,-1);
340  };  };
341    
342  void CaloLong::Draw(){  void Calo2D::Draw(){
343    //    //
344    Process();    Process();
345    Draw(-1);    Draw(-1);
346  };  };
347    
   
348  void CaloLat::Draw(Int_t view,Int_t plane){  void CaloLat::Draw(Int_t view,Int_t plane){
349    //    //
350    Int_t minv = 0;    Int_t minv = 0;
# Line 203  void CaloLat::Draw(Int_t view,Int_t plan Line 407  void CaloLat::Draw(Int_t view,Int_t plan
407    //    //
408  };  };
409    
410  void CaloLong::Draw(Int_t view){  void Calo2D::Draw(Int_t plane){
411    //    //
412    Int_t minv = 0;    Int_t minp = 0;
413    Int_t maxv = 0;    Int_t maxp = 0;
414    //    //
415    if ( view == -1 ){    if ( plane == -1 ){
416          maxv = -1;          minp = 0;
417            maxp = 23;
418    } else {    } else {
419          minv = view;          minp = plane;
420          maxv = view+1;          maxp = plane+1;
421    };    };
422    //    //
423    Process();    Process();
# Line 220  void CaloLong::Draw(Int_t view){ Line 425  void CaloLong::Draw(Int_t view){
425    gStyle->SetLabelSize(0.04);    gStyle->SetLabelSize(0.04);
426    gStyle->SetNdivisions(510,"XY");    gStyle->SetNdivisions(510,"XY");
427    //    //
428    if ( maxv != -1 ){    for (Int_t p=minp; p<maxp;p++){
429    for (Int_t v=minv; v<maxv;v++){      TString hid = Form("c2dp%i",p);    
430          TString hid = Form("clongv%i",v);            TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
431          TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));      if ( tc ){
432          if ( tc ){        //         tc->Clear();
433  //       tc->Clear();      } else {
434          } else {        tc = new TCanvas(hid,hid);
435           tc = new TCanvas(hid,hid);      };
436          };      //
437          //      TString thid = Form("h2dp%i",p);    
438          TString thid = Form("hlongv%i",v);            TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
439          TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));      if ( th ) th->Delete();
440          if ( th ) th->Delete();      //   th->Clear();
441  //       th->Clear();      //   th->Reset();
442  //       th->Reset();      //  } else {
443  //      } else {      Int_t minx = smax[p] - 10;
444           th = new TH1F(thid,thid,22,-0.5,21.5);      if ( minx < 0 ) minx = 0;
445  //      };      Int_t maxx = minx + 20;
446          tc->cd();      if ( maxx > 95 ){
447          //        maxx = 95;
448          for (Int_t st=0;st<22;st++){        minx = 75;
449           th->Fill(st,eplane[v][st]);      };
450          };      Int_t miny = smay[p] - 10;
451          th->Draw();      if ( miny < 0 ) miny = 0;
452          tc->Modified();      Int_t maxy = miny + 20;
453          tc->Update();      if ( maxy > 95 ){
454    };        maxy = 95;
455    } else {        miny = 75;
456          //      };
457          TString hid = Form("clongvyvx");              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);
458          TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));      //    th = new TH2F(thid,thid,96,-0.5,95.5,96,-0.5,95.5);
459          if ( tc ){      //  };
460          } else {      tc->cd();
461           tc = new TCanvas(hid,hid);      //
462          };      for (Int_t stx=minx;stx<maxx+1;stx++){
463          //        for (Int_t sty=miny;sty<maxy+1;sty++){
464          TString thid = Form("hlongvyvx");                th->Fill(stx,sty,estrip[p][stx][sty]);
465          TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));        };
466          if ( th ) th->Delete();      };
467          th = new TH1F(thid,thid,44,-0.5,43.5);      gStyle->SetPalette(1);
468          tc->cd();      //    tc->SetLogz();    
469          Int_t pp=0;      //    th->Draw("colbox");
470          for (Int_t st=0;st<22;st++){      th->Draw("cont4");
471            for (Int_t v=1; v>=0;v--){      tc->Modified();
472                  //      tc->Update();
                 th->Fill(pp,eplane[v][st]);  
                 //  
                 pp++;  
           };  
         };  
         th->Draw();  
         tc->Modified();  
         tc->Update();  
473    };    };
474    //    //
475    gStyle->SetLabelSize(0);    gStyle->SetLabelSize(0);
# Line 285  void CaloLat::Delete(){ Line 482  void CaloLat::Delete(){
482    //delete this;    //delete this;
483  };  };
484    
485  void CaloLong::Delete(){  void Calo2D::Delete(){
486    Clear();    Clear();
487    //delete this;    //delete this;
488  };  };
489    
490    
491  void CaloLat::Process(){  void CaloLat::Process(){
492    //      //  
493    if ( !L2 ){    if ( !L2 ){
# Line 335  void CaloLat::Process(){ Line 533  void CaloLat::Process(){
533    //    //
534  };  };
535    
536  void CaloLong::Process(){  void Calo2D::Process(){
537    //      //  
538    if ( !L2 ){    if ( !L2 ){
539      printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");      printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
# Line 364  void CaloLong::Process(){ Line 562  void CaloLong::Process(){
562    //    //
563    // let's start    // let's start
564    //    //
565      Float_t es[2][22][96];
566      memset(es,0, 4224*sizeof(Float_t));
567      memset(estrip,0, 4224*sizeof(Float_t));
568      Float_t mip1 = 0.;
569      Int_t view1 = 0;
570      Int_t plane1 = 0;
571      Int_t strip1 = 0;
572      //
573      for (Int_t i=0; i<L2->GetCaloLevel1()->istrip ; i++){
574        mip1 = L2->GetCaloLevel1()->DecodeEstrip(i,view1,plane1,strip1);
575        es[view1][plane1][strip1] = mip1;
576      };
577      //
578      Int_t plane2 = 0;
579      Float_t emax[23];
580      memset(emax,0,sizeof(Float_t)*23);
581      memset(smax,0,sizeof(Int_t)*23);
582      memset(smay,0,sizeof(Int_t)*23);
583      //
584      //
585      for (Int_t p=0; p < 23 ; p++){
586        //
587        plane1 = p-1;
588        plane2 = p;
589        //
590        if ( p == 0 ){
591          plane1 = -1;
592          plane2 = 0;
593        };
594        if ( p == 22 ){
595          plane1 = 21;
596          plane2 = -1;
597        };
598        //
599        for (Int_t s=0; s < 96 ; s++){ // x
600          for (Int_t ss=0; ss < 96 ; ss++){ // y
601            if ( p == 0 ){
602              estrip[p][s][ss] += es[1][plane2][ss];
603              if ( (es[1][plane2][ss]) > emax[p] ){
604                smax[p] = 45;
605                smay[p] = ss;
606                emax[p] = es[1][plane2][ss] ;
607              };
608            };
609            if ( p > 0 && p < 22 ){
610              estrip[p][s][ss] += es[0][plane1][s] + es[1][plane2][ss];
611              if ( (es[0][plane1][s] + es[1][plane2][ss]) > emax[p] ){
612                smax[p] = s;
613                smay[p] = ss;
614                emax[p] = es[0][plane1][s] + es[1][plane2][ss] ;
615              };
616            };      
617            if ( p == 22 ){
618              estrip[p][s][ss] += es[0][plane1][s];
619              if ( (es[1][plane2][s]) > emax[p] ){
620                smax[p] = s;
621                smay[p] = 45;
622                emax[p] = es[1][plane2][s] ;
623              };
624            };
625          };
626        };
627        //
628      };
629      //
630      if ( debug ) this->Print();
631      if ( debug ) printf(" exit \n");
632      //
633    };
634    
635    
636    /**
637     * Default constructor
638     */
639    CaloLong::CaloLong(){
640      Clear();
641    };
642    
643    CaloLong::CaloLong(PamLevel2 *l2p){  
644      //
645      Clear();
646      //
647      L2 = l2p;
648      //
649      if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
650      //
651      OBT = 0;
652      PKT = 0;
653      atime = 0;
654      //
655      clp = L2->GetCaloLevel2();
656      //
657      sel = true;
658      cont = false;
659      N = 0;
660      NC = 22;
661      mask18b = -1;
662      //
663      no18x = true;
664      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      xyaverage = true;
675      //
676    };
677    
678    void CaloLong::MaskSection(TString sec){
679      sec.ToUpper();
680      if ( sec.Contains("XO") ) maskXO = true;
681      if ( sec.Contains("YO") ) maskYO = true;
682      if ( sec.Contains("XE") ) maskXE = true;
683      if ( sec.Contains("YE") ) maskYE = true;
684    }
685    
686    void CaloLong::UnMaskSections(){
687      this->UnMaskSection("XEXOYEYO");
688    }
689    
690    void CaloLong::UnMaskSection(TString sec){
691      sec.ToUpper();
692      if ( sec.Contains("XO") ) maskXO = false;
693      if ( sec.Contains("YO") ) maskYO = false;
694      if ( sec.Contains("XE") ) maskXE = false;
695      if ( sec.Contains("YE") ) maskYE = false;
696    }
697    
698    void CaloLong::Clear(){
699      //
700    memset(eplane,0, 2*22*sizeof(Float_t));    memset(eplane,0, 2*22*sizeof(Float_t));
701    //    //
702      chi2 = 0.;
703      ndf = 0.;
704      E0 = 0.;
705      defE0 = 0.;
706      a = 0.;
707      b = 0.;
708      errE0 = 0.;
709      erra = 0.;
710      errb = 0.;
711      etmax = 0.;
712      asymm = 0.;
713      fitresult = 0;
714      //
715      X0pl = 0.76;
716      //
717    };
718    
719    void CaloLong::Print(){
720      //
721      Fit();
722      //
723      printf("==================== Calorimeter Longitudinal Profile =======================\n");
724      printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
725      printf(" fitresult:.. %i\n",fitresult);
726      printf(" chi2     :.. %f\n",chi2);
727      printf(" ndf      :.. %f\n",ndf);
728      printf(" nchi2    :.. %f\n",chi2/ndf);
729      printf(" E0       :.. %f\n",E0);
730      printf(" E0/260.  :.. %f\n",E0/260.);
731      printf(" defE0    :.. %f\n",defE0);
732      printf(" lower    :.. %f\n",lmax);
733      printf(" upper    :.. %f\n",umax);
734      printf(" s lower  :.. %s\n",slmax.Data());
735      printf(" s upper  :.. %s\n",sumax.Data());
736      printf(" a        :.. %f\n",a);
737      printf(" b        :.. %f\n",b);
738      printf(" errE0    :.. %f\n",errE0);
739      printf(" erra     :.. %f\n",erra);
740      printf(" errb     :.. %f\n",errb);
741      printf(" asymm    :.. %f\n",asymm);
742      printf(" tmax     :.. %f\n",((a-1.)/b));
743      printf(" etmax    :.. %f\n",etmax);
744      printf(" X0pl     :.. %f\n",X0pl);
745      printf("========================================================================\n");
746      //
747    };
748    
749    void CaloLong::SetNoWpreSampler(Int_t n){
750      //
751      if ( NC+n <= 22 && NC+n >= 0 ){
752        N = n;
753      } else {
754        printf(" ERROR! Calorimeter is made of 22 W planes\n");
755        printf(" you are giving N presampler = %i and N calo = %i \n",n,NC);
756        printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
757        NC = 22;
758        N = 0;
759      };
760    }
761    
762    void CaloLong::SetNoWcalo(Int_t n){
763      if ( N+n <= 22 && N+n >= 0 ){
764        NC = n;
765      } else {
766        printf(" ERROR! Calorimeter is made of 22 W planes\n");
767        printf(" you are giving N W presampler = %i and N W calo = %i \n",N,n);
768        printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
769        NC = 22;
770        N = 0;
771      };
772    }
773    
774    void CaloLong::SplitInto(Int_t NoWpreSampler, Int_t NoWcalo){
775      this->SetNoWpreSampler(0);
776      this->SetNoWcalo(0);
777      if ( NoWpreSampler < NoWcalo ){
778              this->SetNoWpreSampler(NoWpreSampler);
779              this->SetNoWcalo(NoWcalo);    
780      } else {
781              this->SetNoWcalo(NoWcalo);    
782              this->SetNoWpreSampler(NoWpreSampler);
783      };
784    }
785    
786    void CaloLong::Process(){
787      //  
788      if ( !L2 ){
789        printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
790        printf(" ERROR: CaloHough variables not filled \n");
791        return;
792      };
793      //
794      Bool_t newentry = false;
795      //
796      if ( L2->IsORB() ){
797        if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){
798          newentry = true;
799          OBT = L2->GetOrbitalInfo()->OBT;
800          PKT = L2->GetOrbitalInfo()->pkt_num;
801          atime = L2->GetOrbitalInfo()->absTime;
802        };
803      } else {
804        newentry = true;
805      };
806      //
807      if ( !newentry ) return;
808      //
809      if ( debug ) printf(" Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
810      //
811      Clear();
812      //
813      // let's start
814      //
815      if ( cont ){
816        for (Int_t i=0; i<22; i++){
817          if ( i == (18+N) ){
818            mask18b =  18 + N;
819            break;
820          };
821        };
822      };  
823      //  
824      if ( sel ){
825        for (Int_t i=0; i<22; i++){
826          if ( i == (18-N) ){
827            mask18b =  18 - N;
828            break;
829          };
830        };
831      };
832      //
833      //  if ( mask18b == 18 ) mask18b = -1;
834      //
835      Int_t view = 0;
836      Int_t plane = 0;
837      Int_t strip = 0;
838      Float_t mip = 0.;
839      Bool_t gof = true;
840      for (Int_t i=0; i < L2->GetCaloLevel1()->istrip; i++){
841        mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
842        gof = true;
843        if ( maskXE && (plane%2)==0 && view==1 ) gof = false;
844        if ( maskXO && (plane%2)!=0 && view==1 ) gof = false;
845        if ( maskYE && (plane%2)!=0 && view==0 ) gof = false;
846        if ( maskYO && (plane%2)==0 && view==0 ) gof = false;
847        if ( gof ) eplane[view][plane] += mip;
848      };
849      //
850      // inclination factor (stolen from Daniele's code)
851      //
852      Float_t ytgx = 0;
853      Float_t ytgy = 0;
854      ytgx = 0.76 * clp->tanx[0];
855      ytgy = 0.76 * clp->tany[0];  
856      X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) );
857      //
858      // Find experimental plane of maximum
859      //
860      Int_t pmax = 0;
861      Int_t vmax = 0;
862      Float_t emax = 0.;
863      for (Int_t v=0; v<2; v++){
864       for (Int_t i=0; i<22; i++){
865         if ( eplane[v][i] > emax ){
866           emax = eplane[v][i];
867           vmax = v;
868           pmax = i;
869         };
870       };
871      };
872      //
873      //
874      //
875      if ( vmax == 0 ) pmax++;
876      etmax = pmax * X0pl;
877      //
878      if ( debug ) this->Print();
879      if ( debug ) printf(" exit \n");
880      //
881    };
882    
883    void CaloLong::SetEnergies(Float_t myene[][22]){
884      //  
885      if ( !L2 ){
886        printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
887        printf(" ERROR: CaloHough variables not filled \n");
888        return;
889      };
890      //
891      Bool_t newentry = false;
892      //
893      if ( L2->IsORB() ){
894        if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){
895          newentry = true;
896          OBT = L2->GetOrbitalInfo()->OBT;
897          PKT = L2->GetOrbitalInfo()->pkt_num;
898          atime = L2->GetOrbitalInfo()->absTime;
899        };
900      } else {
901        newentry = true;
902      };
903      //
904      if ( !newentry ) return;
905      //
906      if ( debug ) printf(" SET ENERGIES Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
907      //
908      Clear();
909      //
910      // let's start
911      //
912      if ( cont ){
913        for (Int_t i=0; i<22; i++){
914          if ( i == (18+N) ){
915            mask18b =  18 + N;
916            break;
917          };
918        };
919      };  
920      //  
921      if ( sel ){
922        for (Int_t i=0; i<22; i++){
923          if ( i == (18-N) ){
924            mask18b =  18 - N;
925            break;
926          };
927        };
928      };
929      //
930      //  if ( mask18b == 18 ) mask18b = -1;
931      //
932      Int_t view = 0;
933      Int_t plane = 0;
934      Bool_t gof = true;
935      for (view=0; view < 2; view++){
936        for (plane=0; plane < 22; plane++){      
937          gof = true;
938          if ( maskXE && (plane%2)==0 && view==1 ) gof = false;
939          if ( maskXO && (plane%2)!=0 && view==1 ) gof = false;
940          if ( maskYE && (plane%2)!=0 && view==0 ) gof = false;
941          if ( maskYO && (plane%2)==0 && view==0 ) gof = false;
942          if ( gof ) eplane[view][plane] = myene[view][plane];
943          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]);
944        };
945      };
946      //
947      // inclination factor (stolen from Daniele's code)
948      //
949      Float_t ytgx = 0;
950      Float_t ytgy = 0;
951      ytgx = 0.76 * clp->tanx[0];
952      ytgy = 0.76 * clp->tany[0];  
953      X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) );
954      //
955      // Find experimental plane of maximum
956      //
957      Int_t pmax = 0;
958      Int_t vmax = 0;
959      Float_t emax = 0.;
960    for (Int_t v=0; v<2; v++){    for (Int_t v=0; v<2; v++){
961     for (Int_t i=0; i<22; i++){     for (Int_t i=0; i<22; i++){
962       eplane[v][i] = L2->GetCaloLevel1()->qtotpl(v,i);       if ( eplane[v][i] > emax ){
963           emax = eplane[v][i];
964           vmax = v;
965           pmax = i;
966         };
967     };     };
968    };    };
969    //    //
970      //
971      //
972      if ( vmax == 0 ) pmax++;
973      etmax = pmax * X0pl;
974      //
975    if ( debug ) this->Print();    if ( debug ) this->Print();
976    if ( debug ) printf(" exit \n");    if ( debug ) printf(" exit \n");
977    //    //
978  };  };
979    
980    Double_t ccurve(Double_t *ti,Double_t *par){
981      //
982      Double_t t = *ti;
983      Double_t cE0 = par[0];
984      Double_t ca = par[1];
985      Double_t cb = par[2];
986      Double_t gammaa = TMath::Gamma(ca);
987      //
988      Double_t value = cE0 * cb * ( pow((cb * t),(ca - 1)) * exp( -cb * t ) ) / gammaa;
989      //
990      return value;
991      //
992    }
993    
994    void CaloLong::Fit(){
995      this->Fit(false);
996    };
997    
998    Float_t CaloLong::Evaluate(TString s, Float_t tmax){
999      /* SAMPLE OUTPUT:
1000         Enter Infix Expression : A + B + C / (E - F)
1001         Postfix Expression is : A B + C E F - / +  
1002      */  
1003      if ( !s.Contains("t") ){
1004        printf(" ERROR, the input formula must contain \"t\"\n");
1005        return 0.;
1006      };
1007      if ( tmax != tmax ){
1008        printf(" ERROR, tmax is nan! \n");
1009        return 0.;
1010      };
1011      TString g=Form("%f",tmax);
1012      TString *ts= new TString("");
1013      ts->Prepend(s.Data());
1014      ts->ReplaceAll("t",1,g.Data(),g.Capacity());
1015      ts->Prepend("(");
1016      ts->Append(")");
1017      if ( debug )  printf(" ts %s tssize %i capac %i s %s g %s \n",ts->Data(),ts->Sizeof(),ts->Capacity(),s.Data(),g.Data());
1018      char in[50],post[50];  
1019      strcpy(&in[0],"                                                 ");
1020      strcpy(&in[0],ts->Data());
1021      strcpy(&post[0],"                                                 ");
1022      if ( debug ) printf("Infix Expression is : ##%s##\n",&in[0]);
1023      infix2postfix(&in[0],&post[0],1);
1024      if ( debug ) printf("Postfix Expression is : ##%s##\n",&post[0]);
1025      /* SAMPLE OUTPUT:
1026         Enter Postfix Expression : 3 5 + 2 /
1027         3 5 + 2 / EQUALS 4
1028      */
1029      Float_t res = evaluate(&post[0]);
1030      if ( debug ) printf("%s EQUALS %f\n",post,res);
1031      return res;
1032    }
1033    
1034    void CaloLong::Fit(Bool_t draw){
1035      //
1036      Process();
1037      //
1038      //  
1039      if ( !L2 ){
1040        printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
1041        printf(" ERROR: CaloHough variables not filled \n");
1042        return;
1043      };
1044      //
1045      Bool_t newentry = false;
1046      //
1047      if ( L2->IsORB() ){
1048        if ( L2->GetOrbitalInfo()->pkt_num != fPKT || L2->GetOrbitalInfo()->OBT != fOBT || L2->GetOrbitalInfo()->absTime != fatime ){
1049          newentry = true;
1050          fOBT = L2->GetOrbitalInfo()->OBT;
1051          fPKT = L2->GetOrbitalInfo()->pkt_num;
1052          fatime = L2->GetOrbitalInfo()->absTime;
1053        };
1054      } else {
1055        newentry = true;
1056      };
1057      //
1058      if ( !newentry ) return;
1059      //
1060      if ( debug ) printf(" Start fitting event at OBT %u PKT %u time %u \n",fOBT,fPKT,fatime);
1061      //
1062      if ( draw ){
1063        gStyle->SetLabelSize(0.04);
1064        gStyle->SetNdivisions(510,"XY");
1065      };
1066      //
1067      TString hid = Form("clongfit");      
1068      TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
1069      //  if ( tc ) tc->Delete();
1070      //  if ( tc ) tc->Close();
1071      if ( !tc && draw ){
1072        tc = new TCanvas(hid,hid);
1073      } else {
1074        if ( tc ) tc->cd();
1075      };
1076      //
1077      TString thid = Form("hlongfit");      
1078      TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
1079      if ( th ) th->Delete();
1080      //
1081      TString ghid = Form("glongfit");      
1082      TGraphErrors *gh  = dynamic_cast<TGraphErrors*>(gDirectory->FindObject(ghid));
1083      if ( gh ) gh->Delete();
1084      //
1085      Float_t xpos = 0.;
1086      Float_t enemip = 0.;
1087      Float_t xmax = NC * X0pl + 0.2;
1088      //
1089      Double_t xxx[50];
1090      Double_t exx[50];
1091      Double_t yyy[50];
1092      Double_t eyy[50];
1093      Int_t numpo = 0;
1094      memset(xxx,0,50*sizeof(Double_t));
1095      memset(exx,0,50*sizeof(Double_t));
1096      memset(yyy,0,50*sizeof(Double_t));
1097      memset(eyy,0,50*sizeof(Double_t));
1098      //
1099      // AGH, BUG!
1100      //
1101      Int_t mmin = 0;
1102      Int_t mmax = 0;
1103      if ( cont ){
1104        mmin = N;
1105        mmax = NC+N;
1106      } else {
1107        mmin = 0;
1108        mmax = NC;
1109      };
1110      //
1111      Float_t emax = 0.;
1112      Float_t qtotparz = 0.;
1113      for (Int_t st=mmin;st<mmax+1;st++){
1114        enemip = 0.;
1115        xpos = (st - mmin) * X0pl;
1116        //
1117        if ( xyaverage ){
1118          //
1119          if ( st > mmin && st < mmax ){      
1120            if ( no18x && ( st == 18+1 || st == mask18b+1 )){
1121              if ( !maskYO ){
1122                enemip = 2. * eplane[1][st];
1123              } else {
1124                enemip = eplane[1][st];
1125              };
1126            } else {
1127              enemip = eplane[0][st-1] + eplane[1][st];
1128            };
1129          } else {
1130            if ( st == mmin ){
1131              if ( !maskYE ){
1132                enemip = 2. * eplane[1][st];
1133              } else {
1134                enemip = eplane[1][st];
1135              };
1136            };
1137            if ( st == mmax ){
1138              if ( !maskXE ){
1139                enemip = 2. * eplane[0][st-1];
1140              } else {
1141                enemip = eplane[0][st-1];
1142              };
1143            };
1144          };
1145          //
1146          qtotparz += enemip;
1147          if ( enemip > emax ) emax = enemip;
1148          if ( enemip > 0. ){
1149            xxx[numpo] = xpos;
1150            exx[numpo] = 0.1;
1151            yyy[numpo] = enemip;
1152            eyy[numpo] = sqrt(enemip*3.)+sqrt(5.);
1153            numpo++;
1154            //      th->Fill(xpos,enemip);
1155            if ( debug ) printf(" AVE Filling: st %i xpos %f energy %f \n",st,xpos,enemip);
1156          };
1157        } else {
1158          for (Int_t jj=0; jj<2; jj++){
1159            if ( st > mmin && st < mmax ){
1160              if ( jj == 0 && no18x && ( st == 18+1 || st == mask18b+1 )){
1161                enemip = eplane[1][st];
1162              } else {
1163                if ( jj == 0 ){
1164                  enemip = eplane[jj][st-1];
1165                } else {
1166                  enemip = eplane[jj][st];
1167                };        
1168              };
1169            } else {
1170              if ( st == mmin && jj == 1 ){
1171                enemip = eplane[jj][st];
1172              };
1173              if ( st == mmax && jj == 0){
1174                enemip = eplane[jj][st-1];
1175              };
1176            };
1177            //
1178            qtotparz += enemip;
1179          if ( enemip > emax ) emax = enemip;
1180            if ( enemip > 0. ){
1181              xxx[numpo] = xpos;
1182              exx[numpo] = 0.1;
1183              yyy[numpo] = enemip;
1184              eyy[numpo] = sqrt(enemip*3.)+sqrt(5.);
1185              //      eyy[numpo] = sqrt(enemip)/(st*0.95);
1186              numpo++;
1187              //      th->Fill(xpos,enemip);
1188              if ( debug ) printf(" XY Filling: st %i xpos %f energy %f \n",st,xpos,enemip);
1189            };
1190          };        
1191        };
1192    
1193        //
1194        //    for (Int_t v=1; v>=0;v--)// {
1195    //       //
1196    //       if ( v == 1 ){
1197    //      xpos = (st - N) * X0pl;
1198    //       } else {
1199    //      xpos = (st + 1 - N) * X0pl;
1200    //       };
1201    //       //
1202    //       if ( no18x && st == 18 && v == 0 ){
1203    //      // skip plane 18x
1204    //       } else {
1205    //      if ( v == 1 && st == mask18b ){
1206    //        // emulate plane 18x
1207    //      } else {
1208    //        if ( eplane[v][st] > 0. ){
1209    //          th->Fill(xpos,eplane[v][st]);
1210    //          if ( debug ) printf(" Filling: st %i v %i xpos %f energy %f \n",st,v,xpos,eplane[v][st]);
1211    //        };
1212    //      };
1213    //       };
1214    //       //
1215    //     };
1216      };
1217      //
1218      //  th = new TH2F(thid,thid,int(NC*1.5),-0.2,xmax);
1219      th = new TH2F(thid,thid,1000,-0.2,xmax,1000,0.,emax*1.2);
1220      gh = new TGraphErrors(numpo,xxx,yyy,exx,eyy);
1221      TF1 *lfit  = dynamic_cast<TF1*>(gDirectory->FindObject("lfit"));
1222      if ( lfit ) lfit->Delete();
1223      lfit = new TF1("lfit",ccurve,0.,xmax,3);
1224      if ( debug ) printf("qtot %f qtotparz %f \n",clp->qtot,qtotparz);
1225      E0 = qtotparz;
1226      //  E0 = clp->qtot;
1227      a = 5.;
1228      b = 0.5;
1229      if ( debug ) printf(" STARTING PARAMETERS: E0 %f a %f b %f \n",E0,a,b);
1230      lfit->SetParameters(E0,a,b);
1231      //  lfit->SetParLimits(0,0.,1000.);
1232      //  lfit->SetParLimits(1,-1.,80.);
1233      //  lfit->SetParLimits(2,-1.,10.);
1234      //  TString optio = "ROW"; // "RO"
1235      TString optio = "RO"; // "RO"
1236      if ( debug ){
1237        optio += "NV";
1238        if ( draw ) optio += "V";
1239      } else {
1240        optio += "NQ";
1241        if ( draw ) optio += "Q";
1242      };
1243      //
1244      if ( debug ) printf(" OK, start the fitting procedure...\n");
1245      //
1246      //  fitresult = th->Fit("lfit",optio);
1247      fitresult = gh->Fit("lfit",optio);
1248      //
1249      if ( debug ) printf(" the fit is done! result: %i \n",fitresult);
1250      //
1251      E0 = lfit->GetParameter(0);
1252      a = lfit->GetParameter(1);
1253      b = lfit->GetParameter(2);
1254      errE0 = lfit->GetParError(0);
1255      erra = lfit->GetParError(1);
1256      errb = lfit->GetParError(2);
1257      chi2 = lfit->GetChisquare();
1258      ndf = lfit->GetNDF();
1259      Float_t tmax = 0.;
1260      if ( debug ) printf(" Parameters are retrieved \n");
1261      if ( b != 0 ) tmax = (a - 1.)/b;
1262      //
1263      if ( fitresult != 0 ){
1264        if ( debug ) printf(" The fit failed, no integrals calculation and asymm is set to -1. \n");
1265        asymm = -1.;
1266        defE0 = -1.;
1267      } else {
1268        if ( slmax.MaybeRegexp() ) lmax = this->Evaluate(slmax,tmax);
1269        if ( sumax.MaybeRegexp() ) umax = this->Evaluate(sumax,tmax);
1270        Int_t npp = 1000;
1271        double *xpp=new double[npp];
1272        double *wpp=new double[npp];
1273        lfit->CalcGaussLegendreSamplingPoints(npp,xpp,wpp,1e-12);
1274        defE0 = lfit->IntegralFast(npp,xpp,wpp,lmax,umax);
1275        //
1276        double *xp=new double[npp];
1277        double *wp=new double[npp];
1278        lfit->CalcGaussLegendreSamplingPoints(npp,xp,wp,1e-12);
1279        Float_t imax = lfit->IntegralFast(npp,xp,wp,0.,tmax);
1280        //    Float_t imax = lfit->Integral(0.,tmax);
1281        if ( debug ) printf(" Integral till maximum (%f): %f \n",tmax,imax);
1282        Int_t np = 1000;
1283        double *x=new double[np];
1284        double *w=new double[np];
1285        lfit->CalcGaussLegendreSamplingPoints(np,x,w,1e-12);
1286        Float_t i10max = lfit->IntegralFast(np,x,w,0.,10.*tmax);
1287        delete x;
1288        delete w;
1289        delete xp;
1290        delete wp;
1291        //    Float_t i10max = lfit->Integral(0.,10.*tmax);
1292        if ( debug ) printf(" Integral: %f \n",i10max);
1293        //
1294        if ( i10max != imax ){
1295          asymm = imax / (i10max-imax);
1296        } else {
1297          if ( debug ) printf(" i10max == imax, asymm undefined\n");
1298          asymm = -2.;
1299        };
1300        if ( asymm != asymm ){
1301          if ( debug ) printf(" asymm is nan \n");      
1302          asymm = -3.;
1303        };
1304        //lfit->Integral(0.,tmax)/(lfit->Integral(0.,10.*tmax)-lfit->Integral(0.,tmax));
1305        if ( debug ) printf(" Asymmetry has been calculated \n");
1306      };
1307      //
1308      if ( asymm < 0. || ndf <= 0. || chi2 < 0. || tmax < 0. ){
1309        if ( debug ) printf(" Funny asymm||ndf||chi2||tmax values, fit failed \n");
1310        fitresult = 100;
1311      };
1312      //
1313      if ( draw ){
1314        //
1315        tc->cd();    
1316        //    gStyle->SetOptStat(11111);
1317        tc->SetTitle();
1318        th->SetTitle("");
1319        th->SetName("");
1320        th->SetMarkerStyle(20);
1321        // axis titles
1322        th->SetXTitle("Depth [X0]");
1323        th->SetYTitle("Energy [MIP]");
1324        //    th->DrawCopy("Perror");
1325        th->DrawCopy();
1326        gh->SetMarkerSize(1.);
1327        gh->SetMarkerStyle(21);
1328        //    gh->SetMarkerColor(kRed);
1329        gh->Draw("Psame");
1330        lfit->Draw("same");
1331        tc->Modified();
1332        tc->Update();
1333        //
1334        gStyle->SetLabelSize(0);
1335        gStyle->SetNdivisions(1,"XY");
1336        //
1337      } else {
1338        if ( th ) th->Delete();
1339      };
1340      //
1341      //  delete lfit;
1342      //
1343    };
1344    
1345    void CaloLong::Draw(){
1346      //
1347      Process();
1348      Draw(-1);
1349    };
1350    
1351    void CaloLong::Draw(Int_t view){
1352      //
1353      Int_t minv = 0;
1354      Int_t maxv = 0;
1355      //
1356      if ( view == -1 ){
1357        maxv = -1;
1358      } else {
1359        minv = view;
1360        maxv = view+1;
1361      };
1362      //
1363      Process();
1364      //
1365      gStyle->SetLabelSize(0.04);
1366      gStyle->SetNdivisions(510,"XY");
1367      //
1368      if ( maxv != -1 ){
1369        for (Int_t v=minv; v<maxv;v++){
1370          TString hid = Form("clongv%i",v);
1371          TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
1372          if ( tc ){
1373            //       tc->Clear();
1374          } else {
1375            tc = new TCanvas(hid,hid);
1376          };
1377          //
1378          TString thid = Form("hlongv%i",v);        
1379          TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
1380          if ( th ) th->Delete();
1381          //         th->Clear();
1382          //         th->Reset();
1383          //        } else {
1384          th = new TH1F(thid,thid,22,-0.5,21.5);
1385          //        };
1386          tc->cd();
1387          //
1388          for (Int_t st=0;st<22;st++){
1389            th->Fill(st,eplane[v][st]);
1390          };
1391          th->Draw();
1392          tc->Modified();
1393          tc->Update();
1394        };
1395      } else {
1396        //
1397        TString hid = Form("clongvyvx");    
1398        TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
1399        if ( tc ){
1400        } else {
1401          tc = new TCanvas(hid,hid);
1402        };
1403        //
1404        TString thid = Form("hlongvyvx");  
1405        TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
1406        if ( th ) th->Delete();
1407        th = new TH1F(thid,thid,44,-0.5,43.5);
1408        tc->cd();
1409        Int_t pp=0;
1410        for (Int_t st=0;st<22;st++){
1411          for (Int_t v=1; v>=0;v--){
1412            //
1413            th->Fill(pp,eplane[v][st]);
1414            //
1415            pp++;
1416          };
1417        };
1418        th->Draw();
1419        tc->Modified();
1420        tc->Update();
1421      };
1422      //
1423      gStyle->SetLabelSize(0);
1424      gStyle->SetNdivisions(1,"XY");
1425      //
1426    };
1427    
1428    void CaloLong::Delete(){
1429      Clear();
1430      //delete this;
1431    };
1432    
1433    
1434    
1435    
1436    
1437    
1438    
1439    
1440    
1441    
1442    
1443    
1444    
1445    
1446    
1447    
1448    
1449    
1450    
1451    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.23