/[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.7 by mocchiut, Tue Aug 4 15:05:21 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  /**  /**
# Line 60  CaloLat::CaloLat(){ Line 64  CaloLat::CaloLat(){
64  /**  /**
65   * Default constructor   * Default constructor
66   */   */
67  CaloLong::CaloLong(){  Calo2D::Calo2D(){
68    Clear();    Clear();
69  };  };
70    
# Line 80  CaloLat::CaloLat(PamLevel2 *l2p){   Line 84  CaloLat::CaloLat(PamLevel2 *l2p){  
84    //    //
85  };  };
86    
87  CaloLong::CaloLong(PamLevel2 *l2p){    Calo2D::Calo2D(PamLevel2 *l2p){  
88    //    //
89    Clear();    Clear();
90    //    //
# Line 100  void CaloLat::Clear(){ Line 104  void CaloLat::Clear(){
104    //    //
105  };  };
106    
107  void CaloLong::Clear(){  void Calo2D::Clear(){
108    //    //
109  };  };
110    
# Line 116  void CaloLat::Print(){ Line 120  void CaloLat::Print(){
120    //    //
121  };  };
122    
123  void CaloLong::Print(){  void Calo2D::Print(){
124    //    //
125    Process();    Process();
126    //    //
127    printf("==================== Calorimeter Lateral Profile =======================\n");    printf("==================== Calorimeter 2D Profile =======================\n");
128    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
129  //  printf(" nx [number of X combination]:.. %i\n",nx);  //  printf(" nx [number of X combination]:.. %i\n",nx);
130  //  printf(" ny [number of Y combination]:.. %i\n",ny);  //  printf(" ny [number of Y combination]:.. %i\n",ny);
# Line 134  void CaloLat::Draw(){ Line 138  void CaloLat::Draw(){
138    Draw(-1,-1);    Draw(-1,-1);
139  };  };
140    
141  void CaloLong::Draw(){  void Calo2D::Draw(){
142    //    //
143    Process();    Process();
144    Draw(-1);    Draw(-1);
145  };  };
146    
   
147  void CaloLat::Draw(Int_t view,Int_t plane){  void CaloLat::Draw(Int_t view,Int_t plane){
148    //    //
149    Int_t minv = 0;    Int_t minv = 0;
# Line 203  void CaloLat::Draw(Int_t view,Int_t plan Line 206  void CaloLat::Draw(Int_t view,Int_t plan
206    //    //
207  };  };
208    
209  void CaloLong::Draw(Int_t view){  void Calo2D::Draw(Int_t plane){
210    //    //
211    Int_t minv = 0;    Int_t minp = 0;
212    Int_t maxv = 0;    Int_t maxp = 0;
213    //    //
214    if ( view == -1 ){    if ( plane == -1 ){
215          maxv = -1;          minp = 0;
216            maxp = 23;
217    } else {    } else {
218          minv = view;          minp = plane;
219          maxv = view+1;          maxp = plane+1;
220    };    };
221    //    //
222    Process();    Process();
# Line 220  void CaloLong::Draw(Int_t view){ Line 224  void CaloLong::Draw(Int_t view){
224    gStyle->SetLabelSize(0.04);    gStyle->SetLabelSize(0.04);
225    gStyle->SetNdivisions(510,"XY");    gStyle->SetNdivisions(510,"XY");
226    //    //
227    if ( maxv != -1 ){    for (Int_t p=minp; p<maxp;p++){
228    for (Int_t v=minv; v<maxv;v++){      TString hid = Form("c2dp%i",p);    
229          TString hid = Form("clongv%i",v);            TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
230          TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));      if ( tc ){
231          if ( tc ){        //         tc->Clear();
232  //       tc->Clear();      } else {
233          } else {        tc = new TCanvas(hid,hid);
234           tc = new TCanvas(hid,hid);      };
235          };      //
236          //      TString thid = Form("h2dp%i",p);    
237          TString thid = Form("hlongv%i",v);            TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
238          TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));      if ( th ) th->Delete();
239          if ( th ) th->Delete();      //   th->Clear();
240  //       th->Clear();      //   th->Reset();
241  //       th->Reset();      //  } else {
242  //      } else {      Int_t minx = smax[p] - 10;
243           th = new TH1F(thid,thid,22,-0.5,21.5);      if ( minx < 0 ) minx = 0;
244  //      };      Int_t maxx = minx + 20;
245          tc->cd();      if ( maxx > 95 ){
246          //        maxx = 95;
247          for (Int_t st=0;st<22;st++){        minx = 75;
248           th->Fill(st,eplane[v][st]);      };
249          };      Int_t miny = smay[p] - 10;
250          th->Draw();      if ( miny < 0 ) miny = 0;
251          tc->Modified();      Int_t maxy = miny + 20;
252          tc->Update();      if ( maxy > 95 ){
253    };        maxy = 95;
254    } else {        miny = 75;
255          //      };
256          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);
257          TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));      //    th = new TH2F(thid,thid,96,-0.5,95.5,96,-0.5,95.5);
258          if ( tc ){      //  };
259          } else {      tc->cd();
260           tc = new TCanvas(hid,hid);      //
261          };      for (Int_t stx=minx;stx<maxx+1;stx++){
262          //        for (Int_t sty=miny;sty<maxy+1;sty++){
263          TString thid = Form("hlongvyvx");                th->Fill(stx,sty,estrip[p][stx][sty]);
264          TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));        };
265          if ( th ) th->Delete();      };
266          th = new TH1F(thid,thid,44,-0.5,43.5);      gStyle->SetPalette(1);
267          tc->cd();      //    tc->SetLogz();    
268          Int_t pp=0;      //    th->Draw("colbox");
269          for (Int_t st=0;st<22;st++){      th->Draw("cont4");
270            for (Int_t v=1; v>=0;v--){      tc->Modified();
271                  //      tc->Update();
                 th->Fill(pp,eplane[v][st]);  
                 //  
                 pp++;  
           };  
         };  
         th->Draw();  
         tc->Modified();  
         tc->Update();  
272    };    };
273    //    //
274    gStyle->SetLabelSize(0);    gStyle->SetLabelSize(0);
# Line 285  void CaloLat::Delete(){ Line 281  void CaloLat::Delete(){
281    //delete this;    //delete this;
282  };  };
283    
284  void CaloLong::Delete(){  void Calo2D::Delete(){
285    Clear();    Clear();
286    //delete this;    //delete this;
287  };  };
288    
289    
290  void CaloLat::Process(){  void CaloLat::Process(){
291    //      //  
292    if ( !L2 ){    if ( !L2 ){
# Line 335  void CaloLat::Process(){ Line 332  void CaloLat::Process(){
332    //    //
333  };  };
334    
335  void CaloLong::Process(){  void Calo2D::Process(){
336    //      //  
337    if ( !L2 ){    if ( !L2 ){
338      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 361  void CaloLong::Process(){
361    //    //
362    // let's start    // let's start
363    //    //
364      Float_t es[2][22][96];
365      memset(es,0, 4224*sizeof(Float_t));
366      memset(estrip,0, 4224*sizeof(Float_t));
367      Float_t mip1 = 0.;
368      Int_t view1 = 0;
369      Int_t plane1 = 0;
370      Int_t strip1 = 0;
371      //
372      for (Int_t i=0; i<L2->GetCaloLevel1()->istrip ; i++){
373        mip1 = L2->GetCaloLevel1()->DecodeEstrip(i,view1,plane1,strip1);
374        es[view1][plane1][strip1] = mip1;
375      };
376      //
377      Int_t plane2 = 0;
378      Float_t emax[23];
379      memset(emax,0,sizeof(Float_t)*23);
380      memset(smax,0,sizeof(Int_t)*23);
381      memset(smay,0,sizeof(Int_t)*23);
382      //
383      //
384      for (Int_t p=0; p < 23 ; p++){
385        //
386        plane1 = p-1;
387        plane2 = p;
388        //
389        if ( p == 0 ){
390          plane1 = -1;
391          plane2 = 0;
392        };
393        if ( p == 22 ){
394          plane1 = 21;
395          plane2 = -1;
396        };
397        //
398        for (Int_t s=0; s < 96 ; s++){ // x
399          for (Int_t ss=0; ss < 96 ; ss++){ // y
400            if ( p == 0 ){
401              estrip[p][s][ss] += es[1][plane2][ss];
402              if ( (es[1][plane2][ss]) > emax[p] ){
403                smax[p] = 45;
404                smay[p] = ss;
405                emax[p] = es[1][plane2][ss] ;
406              };
407            };
408            if ( p > 0 && p < 22 ){
409              estrip[p][s][ss] += es[0][plane1][s] + es[1][plane2][ss];
410              if ( (es[0][plane1][s] + es[1][plane2][ss]) > emax[p] ){
411                smax[p] = s;
412                smay[p] = ss;
413                emax[p] = es[0][plane1][s] + es[1][plane2][ss] ;
414              };
415            };      
416            if ( p == 22 ){
417              estrip[p][s][ss] += es[0][plane1][s];
418              if ( (es[1][plane2][s]) > emax[p] ){
419                smax[p] = s;
420                smay[p] = 45;
421                emax[p] = es[1][plane2][s] ;
422              };
423            };
424          };
425        };
426        //
427      };
428      //
429      if ( debug ) this->Print();
430      if ( debug ) printf(" exit \n");
431      //
432    };
433    
434    
435    /**
436     * Default constructor
437     */
438    CaloLong::CaloLong(){
439      Clear();
440    };
441    
442    CaloLong::CaloLong(PamLevel2 *l2p){  
443      //
444      Clear();
445      //
446      L2 = l2p;
447      //
448      if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
449      //
450      OBT = 0;
451      PKT = 0;
452      atime = 0;
453      //
454      sel = true;
455      cont = false;
456      N = 0;
457      NC = 22;
458      mask18b = -1;
459      //
460      no18x = true;
461      debug = false;
462      maskXE = false;
463      maskXO = false;
464      maskYE = false;
465      maskYO = false;
466      //
467    };
468    
469    void CaloLong::MaskSection(TString sec){
470      sec.ToUpper();
471      if ( sec.Contains("XO") ) maskXO = true;
472      if ( sec.Contains("YO") ) maskYO = true;
473      if ( sec.Contains("XE") ) maskXE = true;
474      if ( sec.Contains("YE") ) maskYE = true;
475    }
476    
477    void CaloLong::Clear(){
478      //
479    memset(eplane,0, 2*22*sizeof(Float_t));    memset(eplane,0, 2*22*sizeof(Float_t));
480    //    //
481      chi2 = 0.;
482      ndf = 0.;
483      E0 = 0.;
484      a = 0.;
485      b = 0.;
486      errE0 = 0.;
487      erra = 0.;
488      errb = 0.;
489      etmax = 0.;
490      asymm = 0.;
491      fitresult = 0;
492      //
493      X0pl = 0.76;
494      //
495    };
496    
497    void CaloLong::Print(){
498      //
499      Fit();
500      //
501      printf("==================== Calorimeter Longitudinal Profile =======================\n");
502      printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
503      printf(" fitresult:.. %i\n",fitresult);
504      printf(" chi2     :.. %f\n",chi2);
505      printf(" ndf      :.. %f\n",ndf);
506      printf(" nchi2    :.. %f\n",chi2/ndf);
507      printf(" E0       :.. %f\n",E0);
508      printf(" E0/260.  :.. %f\n",E0/260.);
509      printf(" a        :.. %f\n",a);
510      printf(" b        :.. %f\n",b);
511      printf(" errE0    :.. %f\n",errE0);
512      printf(" erra     :.. %f\n",erra);
513      printf(" errb     :.. %f\n",errb);
514      printf(" asymm    :.. %f\n",asymm);
515      printf(" tmax     :.. %f\n",((a-1.)/b));
516      printf(" etmax    :.. %f\n",etmax);
517      printf(" X0pl     :.. %f\n",X0pl);
518      printf("========================================================================\n");
519      //
520    };
521    
522    void CaloLong::SetNoWpreSampler(Int_t n){
523      //
524      if ( NC+n <= 22 && NC+n >= 0 ){
525        N = n;
526      } else {
527        printf(" ERROR! Calorimeter is made of 22 W planes\n");
528        printf(" you are giving N presampler = %i and N calo = %i \n",n,NC);
529        printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
530        NC = 22;
531        N = 0;
532      };
533    }
534    
535    void CaloLong::SetNoWcalo(Int_t n){
536      if ( N+n <= 22 && N+n >= 0 ){
537        NC = n;
538      } else {
539        printf(" ERROR! Calorimeter is made of 22 W planes\n");
540        printf(" you are giving N W presampler = %i and N W calo = %i \n",N,n);
541        printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
542        NC = 22;
543        N = 0;
544      };
545    }
546    
547    void CaloLong::SplitInto(Int_t NoWpreSampler, Int_t NoWcalo){
548      this->SetNoWpreSampler(0);
549      this->SetNoWcalo(0);
550      if ( NoWpreSampler < NoWcalo ){
551              this->SetNoWpreSampler(NoWpreSampler);
552              this->SetNoWcalo(NoWcalo);    
553      } else {
554              this->SetNoWcalo(NoWcalo);    
555              this->SetNoWpreSampler(NoWpreSampler);
556      };
557    }
558    
559    void CaloLong::Process(){
560      //  
561      if ( !L2 ){
562        printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
563        printf(" ERROR: CaloHough variables not filled \n");
564        return;
565      };
566      //
567      Bool_t newentry = false;
568      //
569      if ( L2->IsORB() ){
570        if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){
571          newentry = true;
572          OBT = L2->GetOrbitalInfo()->OBT;
573          PKT = L2->GetOrbitalInfo()->pkt_num;
574          atime = L2->GetOrbitalInfo()->absTime;
575        };
576      } else {
577        newentry = true;
578      };
579      //
580      if ( !newentry ) return;
581      //
582      if ( debug ) printf(" Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
583      //
584      Clear();
585      //
586      // let's start
587      //
588      if ( cont ){
589        for (Int_t i=0; i<22; i++){
590          if ( i == (18+N) ){
591            mask18b =  18 + N;
592            break;
593          };
594        };
595      };  
596      //  
597      if ( sel ){
598        for (Int_t i=0; i<22; i++){
599          if ( i == (18-N) ){
600            mask18b =  18 - N;
601            break;
602          };
603        };
604      };
605      //
606      //  if ( mask18b == 18 ) mask18b = -1;
607      //
608      Int_t view = 0;
609      Int_t plane = 0;
610      Int_t strip = 0;
611      Float_t mip = 0.;
612      Bool_t gof = true;
613      for (Int_t i=0; i < L2->GetCaloLevel1()->istrip; i++){
614        mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
615        gof = true;
616        if ( maskXE && (plane%2)==0 && view==1 ) gof = false;
617        if ( maskXO && (plane%2)!=0 && view==1 ) gof = false;
618        if ( maskYE && (plane%2)!=0 && view==0 ) gof = false;
619        if ( maskYO && (plane%2)==0 && view==0 ) gof = false;
620        if ( gof ) eplane[view][plane] += mip;
621      };
622      //
623      // inclination factor (stolen from Daniele's code)
624      //
625      Float_t ytgx = 0;
626      Float_t ytgy = 0;
627      ytgx = 0.76 * L2->GetCaloLevel2()->tanx[0];
628      ytgy = 0.76 * L2->GetCaloLevel2()->tany[0];  
629      X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) );
630      //
631      // Find experimental plane of maximum
632      //
633      Int_t pmax = 0;
634      Int_t vmax = 0;
635      Float_t emax = 0.;
636    for (Int_t v=0; v<2; v++){    for (Int_t v=0; v<2; v++){
637     for (Int_t i=0; i<22; i++){     for (Int_t i=0; i<22; i++){
638       eplane[v][i] = L2->GetCaloLevel1()->qtotpl(v,i);       if ( eplane[v][i] > emax ){
639           emax = eplane[v][i];
640           vmax = v;
641           pmax = i;
642         };
643     };     };
644    };    };
645    //    //
646      //
647      //
648      if ( vmax == 0 ) pmax++;
649      etmax = pmax * X0pl;
650      //
651    if ( debug ) this->Print();    if ( debug ) this->Print();
652    if ( debug ) printf(" exit \n");    if ( debug ) printf(" exit \n");
653    //    //
654  };  };
655    
656    
657    Double_t ccurve(Double_t *ti,Double_t *par){
658      //
659      Double_t t = *ti;
660      Double_t cE0 = par[0];
661      Double_t ca = par[1];
662      Double_t cb = par[2];
663      Double_t gammaa = TMath::Gamma(ca);
664      //
665      Double_t value = cE0 * cb * ( pow((cb * t),(ca - 1)) * exp( -cb * t ) ) / gammaa;
666      //
667      return value;
668      //
669    }
670    
671    void CaloLong::Fit(){
672      this->Fit(false);
673    };
674    
675    void CaloLong::Fit(Bool_t draw){
676      //
677      Process();
678      //
679      //  
680      if ( !L2 ){
681        printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
682        printf(" ERROR: CaloHough variables not filled \n");
683        return;
684      };
685      //
686      Bool_t newentry = false;
687      //
688      if ( L2->IsORB() ){
689        if ( L2->GetOrbitalInfo()->pkt_num != fPKT || L2->GetOrbitalInfo()->OBT != fOBT || L2->GetOrbitalInfo()->absTime != fatime ){
690          newentry = true;
691          fOBT = L2->GetOrbitalInfo()->OBT;
692          fPKT = L2->GetOrbitalInfo()->pkt_num;
693          fatime = L2->GetOrbitalInfo()->absTime;
694        };
695      } else {
696        newentry = true;
697      };
698      //
699      if ( !newentry ) return;
700      //
701      if ( debug ) printf(" Start fitting event at OBT %u PKT %u time %u \n",fOBT,fPKT,fatime);
702      //
703      if ( draw ){
704        gStyle->SetLabelSize(0.04);
705        gStyle->SetNdivisions(510,"XY");
706      };
707      //
708      TString hid = Form("clongfit");      
709      TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
710      //  if ( tc ) tc->Delete();
711      //  if ( tc ) tc->Close();
712      if ( !tc && draw ){
713        tc = new TCanvas(hid,hid);
714      } else {
715        if ( tc ) tc->cd();
716      };
717      //
718      TString thid = Form("hlongfit");      
719      TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
720      if ( th ) th->Delete();
721      Float_t xpos = 0.;
722      Float_t enemip = 0.;
723      Float_t xmax = NC * X0pl + 0.2;
724      //  th = new TH1F(thid,thid,int(NC*1.5),-0.2,xmax);
725      th = new TH1F(thid,thid,100,-0.2,xmax);
726      //
727      // AGH, BUG!
728      //
729      Int_t mmin = 0;
730      Int_t mmax = 0;
731      if ( cont ){
732        mmin = N;
733        mmax = NC+N;
734      } else {
735        mmin = 0;
736        mmax = NC;
737      };
738      //
739      Float_t qtotparz = 0.;
740      for (Int_t st=mmin;st<mmax+1;st++){
741        enemip = 0.;
742        xpos = (st - mmin) * X0pl;
743        if ( st > mmin && st < mmax ){      
744          if ( no18x && ( st == 18+1 || st == mask18b+1 )){
745            if ( !maskYO ){
746              enemip = 2. * eplane[1][st];
747            } else {
748              enemip = eplane[1][st];
749            };
750          } else {
751            enemip = eplane[0][st-1] + eplane[1][st];
752          };
753        } else {
754          if ( st == mmin ){
755            if ( !maskYE ){
756              enemip = 2. * eplane[1][st];
757            } else {
758              enemip = eplane[1][st];
759            };
760          };
761          if ( st == mmax ){
762            if ( !maskXE ){
763              enemip = 2. * eplane[0][st-1];
764            } else {
765              enemip = eplane[0][st-1];
766            };
767          };
768        };
769        //
770        qtotparz += enemip;
771        if ( enemip > 0. ){
772          th->Fill(xpos,enemip);
773          if ( debug ) printf(" Filling: st %i xpos %f energy %f \n",st,xpos,enemip);
774        };
775        //
776        //    for (Int_t v=1; v>=0;v--)// {
777    //       //
778    //       if ( v == 1 ){
779    //      xpos = (st - N) * X0pl;
780    //       } else {
781    //      xpos = (st + 1 - N) * X0pl;
782    //       };
783    //       //
784    //       if ( no18x && st == 18 && v == 0 ){
785    //      // skip plane 18x
786    //       } else {
787    //      if ( v == 1 && st == mask18b ){
788    //        // emulate plane 18x
789    //      } else {
790    //        if ( eplane[v][st] > 0. ){
791    //          th->Fill(xpos,eplane[v][st]);
792    //          if ( debug ) printf(" Filling: st %i v %i xpos %f energy %f \n",st,v,xpos,eplane[v][st]);
793    //        };
794    //      };
795    //       };
796    //       //
797    //     };
798      };
799      //
800      TF1 *lfit = new TF1("lfit",ccurve,0.,xmax,3);
801      if ( debug ) printf("qtot %f qtotparz %f \n",L2->GetCaloLevel2()->qtot,qtotparz);
802      E0 = qtotparz;
803      //  E0 = L2->GetCaloLevel2()->qtot;
804      a = 5.;
805      b = 0.5;
806      if ( debug ) printf(" STARTING PARAMETERS: E0 %f a %f b %f \n",E0,a,b);
807      lfit->SetParameters(E0,a,b);
808      //  lfit->SetParLimits(0,0.,1000.);
809      //  lfit->SetParLimits(1,-1.,80.);
810      //  lfit->SetParLimits(2,-1.,10.);
811      TString optio;
812      if ( debug ){
813        //    optio = "MERBOV";
814        //    optio = "MEROV";
815        //    optio = "EROV";
816        optio = "RNOV";
817        if ( draw ) optio = "ROV";
818      } else {
819        //    optio = "MERNOQ";
820        //    optio = "ERNOQ";
821        optio = "RNOQ";
822        if ( draw ) optio = "ROQ";
823      };
824      //
825      if ( debug ) printf(" OK, start the fitting procedure...\n");
826      //
827      fitresult = th->Fit("lfit",optio);
828      //
829      if ( debug ) printf(" the fit is done! result: %i \n",fitresult);
830      //
831      E0 = lfit->GetParameter(0);
832      a = lfit->GetParameter(1);
833      b = lfit->GetParameter(2);
834      errE0 = lfit->GetParError(0);
835      erra = lfit->GetParError(1);
836      errb = lfit->GetParError(2);
837      chi2 = lfit->GetChisquare();
838      ndf = lfit->GetNDF();
839      Float_t tmax = 0.;
840      if ( debug ) printf(" Parameters are retrieved \n");
841      if ( b != 0 ) tmax = (a - 1.)/b;
842      //
843      if ( fitresult != 0 ){
844        if ( debug ) printf(" The fit failed, no integrals calculation and asymm is set to -1. \n");
845        asymm = -1.;
846      } else {
847        Int_t npp = 1000;
848        double *xp=new double[npp];
849        double *wp=new double[npp];
850        lfit->CalcGaussLegendreSamplingPoints(npp,xp,wp,1e-12);
851        Float_t imax = lfit->IntegralFast(npp,xp,wp,0.,tmax);
852        //    Float_t imax = lfit->Integral(0.,tmax);
853        if ( debug ) printf(" Integral till maximum (%f): %f \n",tmax,imax);
854        Int_t np = 1000;
855        double *x=new double[np];
856        double *w=new double[np];
857        lfit->CalcGaussLegendreSamplingPoints(np,x,w,1e-12);
858        Float_t i10max = lfit->IntegralFast(np,x,w,0.,10.*tmax);
859        delete x;
860        delete w;
861        delete xp;
862        delete wp;
863        //    Float_t i10max = lfit->Integral(0.,10.*tmax);
864        if ( debug ) printf(" Integral: %f \n",i10max);
865        //
866        if ( i10max != imax ){
867          asymm = imax / (i10max-imax);
868        } else {
869          if ( debug ) printf(" i10max == imax, asymm undefined\n");
870          asymm = -2.;
871        };
872        if ( asymm != asymm ){
873          if ( debug ) printf(" asymm is nan \n");      
874          asymm = -3.;
875        };
876        //lfit->Integral(0.,tmax)/(lfit->Integral(0.,10.*tmax)-lfit->Integral(0.,tmax));
877        if ( debug ) printf(" Asymmetry has been calculated \n");
878      };
879      //
880      if ( asymm < 0. || ndf <= 0. || chi2 < 0. || tmax < 0. ){
881        if ( debug ) printf(" Funny asymm||ndf||chi2||tmax values, fit failed \n");
882        fitresult = 100;
883      };
884      //
885      if ( draw ){
886        //
887        tc->cd();    
888        //    gStyle->SetOptStat(11111);
889        tc->SetTitle();
890        th->SetTitle("");
891        th->SetName("");
892        th->SetMarkerStyle(20);
893        // axis titles
894        th->SetXTitle("Depth [X0]");
895        th->SetYTitle("Energy [MIP]");
896        th->DrawCopy("Perror");
897        lfit->Draw("same");
898        tc->Modified();
899        tc->Update();
900        //
901        gStyle->SetLabelSize(0);
902        gStyle->SetNdivisions(1,"XY");
903        //
904      } else {
905        if ( th ) th->Delete();
906      };
907      //
908      delete lfit;
909      //
910    };
911    
912    void CaloLong::Draw(){
913      //
914      Process();
915      Draw(-1);
916    };
917    
918    void CaloLong::Draw(Int_t view){
919      //
920      Int_t minv = 0;
921      Int_t maxv = 0;
922      //
923      if ( view == -1 ){
924        maxv = -1;
925      } else {
926        minv = view;
927        maxv = view+1;
928      };
929      //
930      Process();
931      //
932      gStyle->SetLabelSize(0.04);
933      gStyle->SetNdivisions(510,"XY");
934      //
935      if ( maxv != -1 ){
936        for (Int_t v=minv; v<maxv;v++){
937          TString hid = Form("clongv%i",v);
938          TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
939          if ( tc ){
940            //       tc->Clear();
941          } else {
942            tc = new TCanvas(hid,hid);
943          };
944          //
945          TString thid = Form("hlongv%i",v);        
946          TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
947          if ( th ) th->Delete();
948          //         th->Clear();
949          //         th->Reset();
950          //        } else {
951          th = new TH1F(thid,thid,22,-0.5,21.5);
952          //        };
953          tc->cd();
954          //
955          for (Int_t st=0;st<22;st++){
956            th->Fill(st,eplane[v][st]);
957          };
958          th->Draw();
959          tc->Modified();
960          tc->Update();
961        };
962      } else {
963        //
964        TString hid = Form("clongvyvx");    
965        TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
966        if ( tc ){
967        } else {
968          tc = new TCanvas(hid,hid);
969        };
970        //
971        TString thid = Form("hlongvyvx");  
972        TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
973        if ( th ) th->Delete();
974        th = new TH1F(thid,thid,44,-0.5,43.5);
975        tc->cd();
976        Int_t pp=0;
977        for (Int_t st=0;st<22;st++){
978          for (Int_t v=1; v>=0;v--){
979            //
980            th->Fill(pp,eplane[v][st]);
981            //
982            pp++;
983          };
984        };
985        th->Draw();
986        tc->Modified();
987        tc->Update();
988      };
989      //
990      gStyle->SetLabelSize(0);
991      gStyle->SetNdivisions(1,"XY");
992      //
993    };
994    
995    void CaloLong::Delete(){
996      Clear();
997      //delete this;
998    };
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.23