/[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.6 by mocchiut, Wed Jun 24 14:12:55 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            enemip = 2. * eplane[1][st];
746          } else {
747            enemip = eplane[0][st-1] + eplane[1][st];
748          };
749        } else {
750          if ( st == mmin ) enemip = 2. * eplane[1][st];
751          if ( st == mmax ) enemip = 2. * eplane[0][st-1];
752        };
753        //
754        qtotparz += enemip;
755        if ( enemip > 0. ){
756          th->Fill(xpos,enemip);
757          if ( debug ) printf(" Filling: st %i xpos %f energy %f \n",st,xpos,enemip);
758        };
759        //
760        //    for (Int_t v=1; v>=0;v--)// {
761    //       //
762    //       if ( v == 1 ){
763    //      xpos = (st - N) * X0pl;
764    //       } else {
765    //      xpos = (st + 1 - N) * X0pl;
766    //       };
767    //       //
768    //       if ( no18x && st == 18 && v == 0 ){
769    //      // skip plane 18x
770    //       } else {
771    //      if ( v == 1 && st == mask18b ){
772    //        // emulate plane 18x
773    //      } else {
774    //        if ( eplane[v][st] > 0. ){
775    //          th->Fill(xpos,eplane[v][st]);
776    //          if ( debug ) printf(" Filling: st %i v %i xpos %f energy %f \n",st,v,xpos,eplane[v][st]);
777    //        };
778    //      };
779    //       };
780    //       //
781    //     };
782      };
783      //
784      TF1 *lfit = new TF1("lfit",ccurve,0.,xmax,3);
785      if ( debug ) printf("qtot %f qtotparz %f \n",L2->GetCaloLevel2()->qtot,qtotparz);
786      E0 = qtotparz;
787      //  E0 = L2->GetCaloLevel2()->qtot;
788      a = 5.;
789      b = 0.5;
790      if ( debug ) printf(" STARTING PARAMETERS: E0 %f a %f b %f \n",E0,a,b);
791      lfit->SetParameters(E0,a,b);
792      //  lfit->SetParLimits(0,0.,1000.);
793      //  lfit->SetParLimits(1,-1.,80.);
794      //  lfit->SetParLimits(2,-1.,10.);
795      TString optio;
796      if ( debug ){
797        //    optio = "MERBOV";
798        //    optio = "MEROV";
799        //    optio = "EROV";
800        optio = "RNOV";
801        if ( draw ) optio = "ROV";
802      } else {
803        //    optio = "MERNOQ";
804        //    optio = "ERNOQ";
805        optio = "RNOQ";
806        if ( draw ) optio = "ROQ";
807      };
808      //
809      if ( debug ) printf(" OK, start the fitting procedure...\n");
810      //
811      fitresult = th->Fit("lfit",optio);
812      //
813      if ( debug ) printf(" the fit is done! result: %i \n",fitresult);
814      //
815      E0 = lfit->GetParameter(0);
816      a = lfit->GetParameter(1);
817      b = lfit->GetParameter(2);
818      errE0 = lfit->GetParError(0);
819      erra = lfit->GetParError(1);
820      errb = lfit->GetParError(2);
821      chi2 = lfit->GetChisquare();
822      ndf = lfit->GetNDF();
823      Float_t tmax = 0.;
824      if ( debug ) printf(" Parameters are retrieved \n");
825      if ( b != 0 ) tmax = (a - 1.)/b;
826      //
827      if ( fitresult != 0 ){
828        if ( debug ) printf(" The fit failed, no integrals calculation and asymm is set to -1. \n");
829        asymm = -1.;
830      } else {
831        Int_t npp = 1000;
832        double *xp=new double[npp];
833        double *wp=new double[npp];
834        lfit->CalcGaussLegendreSamplingPoints(npp,xp,wp,1e-12);
835        Float_t imax = lfit->IntegralFast(npp,xp,wp,0.,tmax);
836        //    Float_t imax = lfit->Integral(0.,tmax);
837        if ( debug ) printf(" Integral till maximum (%f): %f \n",tmax,imax);
838        Int_t np = 1000;
839        double *x=new double[np];
840        double *w=new double[np];
841        lfit->CalcGaussLegendreSamplingPoints(np,x,w,1e-12);
842        Float_t i10max = lfit->IntegralFast(np,x,w,0.,10.*tmax);
843        delete x;
844        delete w;
845        delete xp;
846        delete wp;
847        //    Float_t i10max = lfit->Integral(0.,10.*tmax);
848        if ( debug ) printf(" Integral: %f \n",i10max);
849        //
850        if ( i10max != imax ){
851          asymm = imax / (i10max-imax);
852        } else {
853          if ( debug ) printf(" i10max == imax, asymm undefined\n");
854          asymm = -2.;
855        };
856        if ( asymm != asymm ){
857          if ( debug ) printf(" asymm is nan \n");      
858          asymm = -3.;
859        };
860        //lfit->Integral(0.,tmax)/(lfit->Integral(0.,10.*tmax)-lfit->Integral(0.,tmax));
861        if ( debug ) printf(" Asymmetry has been calculated \n");
862      };
863      //
864      if ( asymm < 0. || ndf <= 0. || chi2 < 0. || tmax < 0. ){
865        if ( debug ) printf(" Funny asymm||ndf||chi2||tmax values, fit failed \n");
866        fitresult = 100;
867      };
868      //
869      if ( draw ){
870        //
871        tc->cd();    
872        //    gStyle->SetOptStat(11111);
873        tc->SetTitle();
874        th->SetTitle("");
875        th->SetName("");
876        th->SetMarkerStyle(20);
877        // axis titles
878        th->SetXTitle("Depth [X0]");
879        th->SetYTitle("Energy [MIP]");
880        th->DrawCopy("Perror");
881        lfit->Draw("same");
882        tc->Modified();
883        tc->Update();
884        //
885        gStyle->SetLabelSize(0);
886        gStyle->SetNdivisions(1,"XY");
887        //
888      } else {
889        if ( th ) th->Delete();
890      };
891      //
892      delete lfit;
893      //
894    };
895    
896    void CaloLong::Draw(){
897      //
898      Process();
899      Draw(-1);
900    };
901    
902    void CaloLong::Draw(Int_t view){
903      //
904      Int_t minv = 0;
905      Int_t maxv = 0;
906      //
907      if ( view == -1 ){
908        maxv = -1;
909      } else {
910        minv = view;
911        maxv = view+1;
912      };
913      //
914      Process();
915      //
916      gStyle->SetLabelSize(0.04);
917      gStyle->SetNdivisions(510,"XY");
918      //
919      if ( maxv != -1 ){
920        for (Int_t v=minv; v<maxv;v++){
921          TString hid = Form("clongv%i",v);
922          TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
923          if ( tc ){
924            //       tc->Clear();
925          } else {
926            tc = new TCanvas(hid,hid);
927          };
928          //
929          TString thid = Form("hlongv%i",v);        
930          TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
931          if ( th ) th->Delete();
932          //         th->Clear();
933          //         th->Reset();
934          //        } else {
935          th = new TH1F(thid,thid,22,-0.5,21.5);
936          //        };
937          tc->cd();
938          //
939          for (Int_t st=0;st<22;st++){
940            th->Fill(st,eplane[v][st]);
941          };
942          th->Draw();
943          tc->Modified();
944          tc->Update();
945        };
946      } else {
947        //
948        TString hid = Form("clongvyvx");    
949        TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
950        if ( tc ){
951        } else {
952          tc = new TCanvas(hid,hid);
953        };
954        //
955        TString thid = Form("hlongvyvx");  
956        TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
957        if ( th ) th->Delete();
958        th = new TH1F(thid,thid,44,-0.5,43.5);
959        tc->cd();
960        Int_t pp=0;
961        for (Int_t st=0;st<22;st++){
962          for (Int_t v=1; v>=0;v--){
963            //
964            th->Fill(pp,eplane[v][st]);
965            //
966            pp++;
967          };
968        };
969        th->Draw();
970        tc->Modified();
971        tc->Update();
972      };
973      //
974      gStyle->SetLabelSize(0);
975      gStyle->SetNdivisions(1,"XY");
976      //
977    };
978    
979    void CaloLong::Delete(){
980      Clear();
981      //delete this;
982    };
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.23