/[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.1.1 by mocchiut, Fri Nov 9 09:11:24 2007 UTC revision 1.10 by mocchiut, Wed Aug 12 14:57:00 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      //
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(){
698      //
699    memset(eplane,0, 2*22*sizeof(Float_t));    memset(eplane,0, 2*22*sizeof(Float_t));
700    //    //
701      chi2 = 0.;
702      ndf = 0.;
703      E0 = 0.;
704      defE0 = 0.;
705      a = 0.;
706      b = 0.;
707      errE0 = 0.;
708      erra = 0.;
709      errb = 0.;
710      etmax = 0.;
711      asymm = 0.;
712      fitresult = 0;
713      //
714      X0pl = 0.76;
715      //
716    };
717    
718    void CaloLong::Print(){
719      //
720      Fit();
721      //
722      printf("==================== Calorimeter Longitudinal Profile =======================\n");
723      printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
724      printf(" fitresult:.. %i\n",fitresult);
725      printf(" chi2     :.. %f\n",chi2);
726      printf(" ndf      :.. %f\n",ndf);
727      printf(" nchi2    :.. %f\n",chi2/ndf);
728      printf(" E0       :.. %f\n",E0);
729      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);
736      printf(" b        :.. %f\n",b);
737      printf(" errE0    :.. %f\n",errE0);
738      printf(" erra     :.. %f\n",erra);
739      printf(" errb     :.. %f\n",errb);
740      printf(" asymm    :.. %f\n",asymm);
741      printf(" tmax     :.. %f\n",((a-1.)/b));
742      printf(" etmax    :.. %f\n",etmax);
743      printf(" X0pl     :.. %f\n",X0pl);
744      printf("========================================================================\n");
745      //
746    };
747    
748    void CaloLong::SetNoWpreSampler(Int_t n){
749      //
750      if ( NC+n <= 22 && NC+n >= 0 ){
751        N = n;
752      } else {
753        printf(" ERROR! Calorimeter is made of 22 W planes\n");
754        printf(" you are giving N presampler = %i and N calo = %i \n",n,NC);
755        printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
756        NC = 22;
757        N = 0;
758      };
759    }
760    
761    void CaloLong::SetNoWcalo(Int_t n){
762      if ( N+n <= 22 && N+n >= 0 ){
763        NC = n;
764      } else {
765        printf(" ERROR! Calorimeter is made of 22 W planes\n");
766        printf(" you are giving N W presampler = %i and N W calo = %i \n",N,n);
767        printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
768        NC = 22;
769        N = 0;
770      };
771    }
772    
773    void CaloLong::SplitInto(Int_t NoWpreSampler, Int_t NoWcalo){
774      this->SetNoWpreSampler(0);
775      this->SetNoWcalo(0);
776      if ( NoWpreSampler < NoWcalo ){
777              this->SetNoWpreSampler(NoWpreSampler);
778              this->SetNoWcalo(NoWcalo);    
779      } else {
780              this->SetNoWcalo(NoWcalo);    
781              this->SetNoWpreSampler(NoWpreSampler);
782      };
783    }
784    
785    void CaloLong::Process(){
786      //  
787      if ( !L2 ){
788        printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
789        printf(" ERROR: CaloHough variables not filled \n");
790        return;
791      };
792      //
793      Bool_t newentry = false;
794      //
795      if ( L2->IsORB() ){
796        if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){
797          newentry = true;
798          OBT = L2->GetOrbitalInfo()->OBT;
799          PKT = L2->GetOrbitalInfo()->pkt_num;
800          atime = L2->GetOrbitalInfo()->absTime;
801        };
802      } else {
803        newentry = true;
804      };
805      //
806      if ( !newentry ) return;
807      //
808      if ( debug ) printf(" Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
809      //
810      Clear();
811      //
812      // let's start
813      //
814      if ( cont ){
815        for (Int_t i=0; i<22; i++){
816          if ( i == (18+N) ){
817            mask18b =  18 + N;
818            break;
819          };
820        };
821      };  
822      //  
823      if ( sel ){
824        for (Int_t i=0; i<22; i++){
825          if ( i == (18-N) ){
826            mask18b =  18 - N;
827            break;
828          };
829        };
830      };
831      //
832      //  if ( mask18b == 18 ) mask18b = -1;
833      //
834      Int_t view = 0;
835      Int_t plane = 0;
836      Int_t strip = 0;
837      Float_t mip = 0.;
838      Bool_t gof = true;
839      for (Int_t i=0; i < L2->GetCaloLevel1()->istrip; i++){
840        mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
841        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)
850      //
851      Float_t ytgx = 0;
852      Float_t ytgy = 0;
853      ytgx = 0.76 * clp->tanx[0];
854      ytgy = 0.76 * clp->tany[0];  
855      X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) );
856      //
857      // Find experimental plane of maximum
858      //
859      Int_t pmax = 0;
860      Int_t vmax = 0;
861      Float_t emax = 0.;
862      for (Int_t v=0; v<2; v++){
863       for (Int_t i=0; i<22; i++){
864         if ( eplane[v][i] > emax ){
865           emax = eplane[v][i];
866           vmax = v;
867           pmax = i;
868         };
869       };
870      };
871      //
872      //
873      //
874      if ( vmax == 0 ) pmax++;
875      etmax = pmax * X0pl;
876      //
877      if ( debug ) this->Print();
878      if ( debug ) printf(" exit \n");
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++){    for (Int_t v=0; v<2; v++){
960     for (Int_t i=0; i<22; i++){     for (Int_t i=0; i<22; i++){
961       eplane[v][i] = L2->GetCaloLevel1()->qtotpl(v,i);       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();    if ( debug ) this->Print();
975    if ( debug ) printf(" exit \n");    if ( debug ) printf(" exit \n");
976    //    //
977  };  };
978    
979    Double_t ccurve(Double_t *ti,Double_t *par){
980      //
981      Double_t t = *ti;
982      Double_t cE0 = par[0];
983      Double_t ca = par[1];
984      Double_t cb = par[2];
985      Double_t gammaa = TMath::Gamma(ca);
986      //
987      Double_t value = cE0 * cb * ( pow((cb * t),(ca - 1)) * exp( -cb * t ) ) / gammaa;
988      //
989      return value;
990      //
991    }
992    
993    void CaloLong::Fit(){
994      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){
1034      //
1035      Process();
1036      //
1037      //  
1038      if ( !L2 ){
1039        printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
1040        printf(" ERROR: CaloHough variables not filled \n");
1041        return;
1042      };
1043      //
1044      Bool_t newentry = false;
1045      //
1046      if ( L2->IsORB() ){
1047        if ( L2->GetOrbitalInfo()->pkt_num != fPKT || L2->GetOrbitalInfo()->OBT != fOBT || L2->GetOrbitalInfo()->absTime != fatime ){
1048          newentry = true;
1049          fOBT = L2->GetOrbitalInfo()->OBT;
1050          fPKT = L2->GetOrbitalInfo()->pkt_num;
1051          fatime = L2->GetOrbitalInfo()->absTime;
1052        };
1053      } else {
1054        newentry = true;
1055      };
1056      //
1057      if ( !newentry ) return;
1058      //
1059      if ( debug ) printf(" Start fitting event at OBT %u PKT %u time %u \n",fOBT,fPKT,fatime);
1060      //
1061      if ( draw ){
1062        gStyle->SetLabelSize(0.04);
1063        gStyle->SetNdivisions(510,"XY");
1064      };
1065      //
1066      TString hid = Form("clongfit");      
1067      TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
1068      //  if ( tc ) tc->Delete();
1069      //  if ( tc ) tc->Close();
1070      if ( !tc && draw ){
1071        tc = new TCanvas(hid,hid);
1072      } else {
1073        if ( tc ) tc->cd();
1074      };
1075      //
1076      TString thid = Form("hlongfit");      
1077      TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
1078      if ( th ) th->Delete();
1079      Float_t xpos = 0.;
1080      Float_t enemip = 0.;
1081      Float_t xmax = NC * X0pl + 0.2;
1082      //  th = new TH1F(thid,thid,int(NC*1.5),-0.2,xmax);
1083      th = new TH1F(thid,thid,100,-0.2,xmax);
1084      //
1085      // AGH, BUG!
1086      //
1087      Int_t mmin = 0;
1088      Int_t mmax = 0;
1089      if ( cont ){
1090        mmin = N;
1091        mmax = NC+N;
1092      } else {
1093        mmin = 0;
1094        mmax = NC;
1095      };
1096      //
1097      Float_t qtotparz = 0.;
1098      for (Int_t st=mmin;st<mmax+1;st++){
1099        enemip = 0.;
1100        xpos = (st - mmin) * X0pl;
1101        if ( st > mmin && st < mmax ){      
1102          if ( no18x && ( st == 18+1 || st == mask18b+1 )){
1103            if ( !maskYO ){
1104              enemip = 2. * eplane[1][st];
1105            } else {
1106              enemip = eplane[1][st];
1107            };
1108          } else {
1109            enemip = eplane[0][st-1] + eplane[1][st];
1110          };
1111        } else {
1112          if ( st == mmin ){
1113            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;
1129        if ( enemip > 0. ){
1130          th->Fill(xpos,enemip);
1131          if ( debug ) printf(" Filling: st %i xpos %f energy %f \n",st,xpos,enemip);
1132        };
1133        //
1134        //    for (Int_t v=1; v>=0;v--)// {
1135    //       //
1136    //       if ( v == 1 ){
1137    //      xpos = (st - N) * X0pl;
1138    //       } else {
1139    //      xpos = (st + 1 - N) * X0pl;
1140    //       };
1141    //       //
1142    //       if ( no18x && st == 18 && v == 0 ){
1143    //      // skip plane 18x
1144    //       } else {
1145    //      if ( v == 1 && st == mask18b ){
1146    //        // emulate plane 18x
1147    //      } else {
1148    //        if ( eplane[v][st] > 0. ){
1149    //          th->Fill(xpos,eplane[v][st]);
1150    //          if ( debug ) printf(" Filling: st %i v %i xpos %f energy %f \n",st,v,xpos,eplane[v][st]);
1151    //        };
1152    //      };
1153    //       };
1154    //       //
1155    //     };
1156      };
1157      //
1158      TF1 *lfit = new TF1("lfit",ccurve,0.,xmax,3);
1159      if ( debug ) printf("qtot %f qtotparz %f \n",clp->qtot,qtotparz);
1160      E0 = qtotparz;
1161      //  E0 = clp->qtot;
1162      a = 5.;
1163      b = 0.5;
1164      if ( debug ) printf(" STARTING PARAMETERS: E0 %f a %f b %f \n",E0,a,b);
1165      lfit->SetParameters(E0,a,b);
1166      //  lfit->SetParLimits(0,0.,1000.);
1167      //  lfit->SetParLimits(1,-1.,80.);
1168      //  lfit->SetParLimits(2,-1.,10.);
1169      TString optio;
1170      if ( debug ){
1171        //    optio = "MERBOV";
1172        //    optio = "MEROV";
1173        //    optio = "EROV";
1174        optio = "RNOV";
1175        if ( draw ) optio = "ROV";
1176      } else {
1177        //    optio = "MERNOQ";
1178        //    optio = "ERNOQ";
1179        optio = "RNOQ";
1180        if ( draw ) optio = "ROQ";
1181      };
1182      //
1183      if ( debug ) printf(" OK, start the fitting procedure...\n");
1184      //
1185      fitresult = th->Fit("lfit",optio);
1186      //
1187      if ( debug ) printf(" the fit is done! result: %i \n",fitresult);
1188      //
1189      E0 = lfit->GetParameter(0);
1190      a = lfit->GetParameter(1);
1191      b = lfit->GetParameter(2);
1192      errE0 = lfit->GetParError(0);
1193      erra = lfit->GetParError(1);
1194      errb = lfit->GetParError(2);
1195      chi2 = lfit->GetChisquare();
1196      ndf = lfit->GetNDF();
1197      Float_t tmax = 0.;
1198      if ( debug ) printf(" Parameters are retrieved \n");
1199      if ( b != 0 ) tmax = (a - 1.)/b;
1200      //
1201      if ( fitresult != 0 ){
1202        if ( debug ) printf(" The fit failed, no integrals calculation and asymm is set to -1. \n");
1203        asymm = -1.;
1204        defE0 = -1.;
1205      } else {
1206        if ( slmax.MaybeRegexp() ) lmax = this->Evaluate(slmax,tmax);
1207        if ( sumax.MaybeRegexp() ) umax = this->Evaluate(sumax,tmax);
1208        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];
1215        double *wp=new double[npp];
1216        lfit->CalcGaussLegendreSamplingPoints(npp,xp,wp,1e-12);
1217        Float_t imax = lfit->IntegralFast(npp,xp,wp,0.,tmax);
1218        //    Float_t imax = lfit->Integral(0.,tmax);
1219        if ( debug ) printf(" Integral till maximum (%f): %f \n",tmax,imax);
1220        Int_t np = 1000;
1221        double *x=new double[np];
1222        double *w=new double[np];
1223        lfit->CalcGaussLegendreSamplingPoints(np,x,w,1e-12);
1224        Float_t i10max = lfit->IntegralFast(np,x,w,0.,10.*tmax);
1225        delete x;
1226        delete w;
1227        delete xp;
1228        delete wp;
1229        //    Float_t i10max = lfit->Integral(0.,10.*tmax);
1230        if ( debug ) printf(" Integral: %f \n",i10max);
1231        //
1232        if ( i10max != imax ){
1233          asymm = imax / (i10max-imax);
1234        } else {
1235          if ( debug ) printf(" i10max == imax, asymm undefined\n");
1236          asymm = -2.;
1237        };
1238        if ( asymm != asymm ){
1239          if ( debug ) printf(" asymm is nan \n");      
1240          asymm = -3.;
1241        };
1242        //lfit->Integral(0.,tmax)/(lfit->Integral(0.,10.*tmax)-lfit->Integral(0.,tmax));
1243        if ( debug ) printf(" Asymmetry has been calculated \n");
1244      };
1245      //
1246      if ( asymm < 0. || ndf <= 0. || chi2 < 0. || tmax < 0. ){
1247        if ( debug ) printf(" Funny asymm||ndf||chi2||tmax values, fit failed \n");
1248        fitresult = 100;
1249      };
1250      //
1251      if ( draw ){
1252        //
1253        tc->cd();    
1254        //    gStyle->SetOptStat(11111);
1255        tc->SetTitle();
1256        th->SetTitle("");
1257        th->SetName("");
1258        th->SetMarkerStyle(20);
1259        // axis titles
1260        th->SetXTitle("Depth [X0]");
1261        th->SetYTitle("Energy [MIP]");
1262        th->DrawCopy("Perror");
1263        lfit->Draw("same");
1264        tc->Modified();
1265        tc->Update();
1266        //
1267        gStyle->SetLabelSize(0);
1268        gStyle->SetNdivisions(1,"XY");
1269        //
1270      } else {
1271        if ( th ) th->Delete();
1272      };
1273      //
1274      delete lfit;
1275      //
1276    };
1277    
1278    void CaloLong::Draw(){
1279      //
1280      Process();
1281      Draw(-1);
1282    };
1283    
1284    void CaloLong::Draw(Int_t view){
1285      //
1286      Int_t minv = 0;
1287      Int_t maxv = 0;
1288      //
1289      if ( view == -1 ){
1290        maxv = -1;
1291      } else {
1292        minv = view;
1293        maxv = view+1;
1294      };
1295      //
1296      Process();
1297      //
1298      gStyle->SetLabelSize(0.04);
1299      gStyle->SetNdivisions(510,"XY");
1300      //
1301      if ( maxv != -1 ){
1302        for (Int_t v=minv; v<maxv;v++){
1303          TString hid = Form("clongv%i",v);
1304          TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
1305          if ( tc ){
1306            //       tc->Clear();
1307          } else {
1308            tc = new TCanvas(hid,hid);
1309          };
1310          //
1311          TString thid = Form("hlongv%i",v);        
1312          TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
1313          if ( th ) th->Delete();
1314          //         th->Clear();
1315          //         th->Reset();
1316          //        } else {
1317          th = new TH1F(thid,thid,22,-0.5,21.5);
1318          //        };
1319          tc->cd();
1320          //
1321          for (Int_t st=0;st<22;st++){
1322            th->Fill(st,eplane[v][st]);
1323          };
1324          th->Draw();
1325          tc->Modified();
1326          tc->Update();
1327        };
1328      } else {
1329        //
1330        TString hid = Form("clongvyvx");    
1331        TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
1332        if ( tc ){
1333        } else {
1334          tc = new TCanvas(hid,hid);
1335        };
1336        //
1337        TString thid = Form("hlongvyvx");  
1338        TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
1339        if ( th ) th->Delete();
1340        th = new TH1F(thid,thid,44,-0.5,43.5);
1341        tc->cd();
1342        Int_t pp=0;
1343        for (Int_t st=0;st<22;st++){
1344          for (Int_t v=1; v>=0;v--){
1345            //
1346            th->Fill(pp,eplane[v][st]);
1347            //
1348            pp++;
1349          };
1350        };
1351        th->Draw();
1352        tc->Modified();
1353        tc->Update();
1354      };
1355      //
1356      gStyle->SetLabelSize(0);
1357      gStyle->SetNdivisions(1,"XY");
1358      //
1359    };
1360    
1361    void CaloLong::Delete(){
1362      Clear();
1363      //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.1.1.1  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.23