/[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.4 by mocchiut, Wed Jun 10 13:00:23 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      //
463    };
464    
465    void CaloLong::Clear(){
466      //
467    memset(eplane,0, 2*22*sizeof(Float_t));    memset(eplane,0, 2*22*sizeof(Float_t));
468    //    //
469      chi2 = 0.;
470      ndf = 0.;
471      E0 = 0.;
472      a = 0.;
473      b = 0.;
474      errE0 = 0.;
475      erra = 0.;
476      errb = 0.;
477      etmax = 0.;
478      asymm = 0.;
479      fitresult = 0;
480      //
481      X0pl = 0.76;
482      //
483    };
484    
485    void CaloLong::Print(){
486      //
487      Process();
488      //
489      printf("==================== Calorimeter Longitudinal Profile =======================\n");
490      printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
491      printf(" fitresult:.. %i\n",fitresult);
492      printf(" chi2     :.. %f\n",chi2);
493      printf(" ndf      :.. %f\n",ndf);
494      printf(" nchi2    :.. %f\n",chi2/ndf);
495      printf(" E0       :.. %f\n",E0);
496      printf(" E0/260.  :.. %f\n",E0/260.);
497      printf(" a        :.. %f\n",a);
498      printf(" b        :.. %f\n",b);
499      printf(" errE0    :.. %f\n",errE0);
500      printf(" erra     :.. %f\n",erra);
501      printf(" errb     :.. %f\n",errb);
502      printf(" asymm    :.. %f\n",asymm);
503      printf(" tmax     :.. %f\n",((a-1.)/b));
504      printf(" etmax    :.. %f\n",etmax);
505      printf(" X0pl     :.. %f\n",X0pl);
506      printf("========================================================================\n");
507      //
508    };
509    
510    void CaloLong::SetNoWpreSampler(Int_t n){
511      //
512      if ( NC+n <= 22 && NC+n >= 0 ){
513        N = n;
514      } else {
515        printf(" ERROR! Calorimeter is made of 22 W planes\n");
516        printf(" you are giving N presampler = %i and N calo = %i \n",n,NC);
517        printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
518        NC = 22;
519        N = 0;
520      };
521    }
522    
523    void CaloLong::SetNoWcalo(Int_t n){
524      if ( N+n <= 22 && N+n >= 0 ){
525        NC = n;
526      } else {
527        printf(" ERROR! Calorimeter is made of 22 W planes\n");
528        printf(" you are giving N W presampler = %i and N W calo = %i \n",N,n);
529        printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
530        NC = 22;
531        N = 0;
532      };
533    }
534    
535    void CaloLong::SplitInto(Int_t NoWpreSampler, Int_t NoWcalo){
536      this->SetNoWpreSampler(0);
537      this->SetNoWcalo(0);
538      if ( NoWpreSampler < NoWcalo ){
539              this->SetNoWpreSampler(NoWpreSampler);
540              this->SetNoWcalo(NoWcalo);    
541      } else {
542              this->SetNoWcalo(NoWcalo);    
543              this->SetNoWpreSampler(NoWpreSampler);
544      };
545    }
546    
547    void CaloLong::Process(){
548      //  
549      if ( !L2 ){
550        printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
551        printf(" ERROR: CaloHough variables not filled \n");
552        return;
553      };
554      //
555      Bool_t newentry = false;
556      //
557      if ( L2->IsORB() ){
558        if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){
559          newentry = true;
560          OBT = L2->GetOrbitalInfo()->OBT;
561          PKT = L2->GetOrbitalInfo()->pkt_num;
562          atime = L2->GetOrbitalInfo()->absTime;
563        };
564      } else {
565        newentry = true;
566      };
567      //
568      if ( !newentry ) return;
569      //
570      if ( debug ) printf(" Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
571      //
572      Clear();
573      //
574      // let's start
575      //
576      if ( cont ){
577        for (Int_t i=0; i<22; i++){
578          if ( i == (18+N) ){
579            mask18b =  18 + N;
580            break;
581          };
582        };
583      };  
584      //  
585      if ( sel ){
586        for (Int_t i=0; i<22; i++){
587          if ( i == (18-N) ){
588            mask18b =  18 - N;
589            break;
590          };
591        };
592      };
593      //
594      //  if ( mask18b == 18 ) mask18b = -1;
595      //
596      Int_t view = 0;
597      Int_t plane = 0;
598      Int_t strip = 0;
599      Float_t mip = 0.;
600      for (Int_t i=0; i < L2->GetCaloLevel1()->istrip; i++){
601        mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
602        eplane[view][plane] += mip;
603      };
604      //
605      // inclination factor (stolen from Daniele's code)
606      //
607      Float_t ytgx = 0;
608      Float_t ytgy = 0;
609      ytgx = 0.76 * L2->GetCaloLevel2()->tanx[0];
610      ytgy = 0.76 * L2->GetCaloLevel2()->tany[0];  
611      X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) );
612      //
613      // Find experimental plane of maximum
614      //
615      Int_t pmax = 0;
616      Int_t vmax = 0;
617      Float_t emax = 0.;
618    for (Int_t v=0; v<2; v++){    for (Int_t v=0; v<2; v++){
619     for (Int_t i=0; i<22; i++){     for (Int_t i=0; i<22; i++){
620       eplane[v][i] = L2->GetCaloLevel1()->qtotpl(v,i);       if ( eplane[v][i] > emax ){
621           emax = eplane[v][i];
622           vmax = v;
623           pmax = i;
624         };
625     };     };
626    };    };
627    //    //
628      //
629      //
630      if ( vmax == 0 ) pmax++;
631      etmax = pmax * X0pl;
632      //
633    if ( debug ) this->Print();    if ( debug ) this->Print();
634    if ( debug ) printf(" exit \n");    if ( debug ) printf(" exit \n");
635    //    //
636  };  };
637    
638    
639    Double_t ccurve(Double_t *ti,Double_t *par){
640      //
641      Double_t t = *ti;
642      Double_t cE0 = par[0];
643      Double_t ca = par[1];
644      Double_t cb = par[2];
645      Double_t gammaa = TMath::Gamma(ca);
646      //
647      Double_t value = cE0 * cb * ( pow((cb * t),(ca - 1)) * exp( -cb * t ) ) / gammaa;
648      //
649      return value;
650      //
651    }
652    
653    void CaloLong::Fit(){
654      this->Fit(false);
655    };
656    
657    void CaloLong::Fit(Bool_t draw){
658      //
659      Process();
660      //
661      //  
662      if ( !L2 ){
663        printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
664        printf(" ERROR: CaloHough variables not filled \n");
665        return;
666      };
667      //
668      Bool_t newentry = false;
669      //
670      if ( L2->IsORB() ){
671        if ( L2->GetOrbitalInfo()->pkt_num != fPKT || L2->GetOrbitalInfo()->OBT != fOBT || L2->GetOrbitalInfo()->absTime != fatime ){
672          newentry = true;
673          fOBT = L2->GetOrbitalInfo()->OBT;
674          fPKT = L2->GetOrbitalInfo()->pkt_num;
675          fatime = L2->GetOrbitalInfo()->absTime;
676        };
677      } else {
678        newentry = true;
679      };
680      //
681      if ( !newentry ) return;
682      //
683      if ( debug ) printf(" Start fitting event at OBT %u PKT %u time %u \n",fOBT,fPKT,fatime);
684      //
685      if ( draw ){
686        gStyle->SetLabelSize(0.04);
687        gStyle->SetNdivisions(510,"XY");
688      };
689      //
690      TString hid = Form("clongfit");      
691      TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
692      //  if ( tc ) tc->Delete();
693      //  if ( tc ) tc->Close();
694      if ( !tc && draw ){
695        tc = new TCanvas(hid,hid);
696      } else {
697        if ( tc ) tc->cd();
698      };
699      //
700      TString thid = Form("hlongfit");      
701      TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
702      if ( th ) th->Delete();
703      Float_t xpos = 0.;
704      Float_t enemip = 0.;
705      Float_t xmax = NC * X0pl + 0.2;
706      //  th = new TH1F(thid,thid,int(NC*1.5),-0.2,xmax);
707      th = new TH1F(thid,thid,100,-0.2,xmax);
708      //
709      // AGH, BUG!
710      //
711      Int_t mmin = 0;
712      Int_t mmax = 0;
713      if ( cont ){
714        mmin = N;
715        mmax = NC+N;
716      } else {
717        mmin = 0;
718        mmax = NC;
719      };
720      //
721      Float_t qtotparz = 0.;
722      for (Int_t st=mmin;st<mmax+1;st++){
723        enemip = 0.;
724        xpos = (st - mmin) * X0pl;
725        if ( st > mmin && st < mmax ){      
726          if ( no18x && ( st == 18+1 || st == mask18b+1 )){
727            enemip = 2. * eplane[1][st];
728          } else {
729            enemip = eplane[0][st-1] + eplane[1][st];
730          };
731        } else {
732          if ( st == mmin ) enemip = 2. * eplane[1][st];
733          if ( st == mmax ) enemip = 2. * eplane[0][st-1];
734        };
735        //
736        qtotparz += enemip;
737        if ( enemip > 0. ){
738          th->Fill(xpos,enemip);
739          if ( debug ) printf(" Filling: st %i xpos %f energy %f \n",st,xpos,enemip);
740        };
741        //
742        //    for (Int_t v=1; v>=0;v--)// {
743    //       //
744    //       if ( v == 1 ){
745    //      xpos = (st - N) * X0pl;
746    //       } else {
747    //      xpos = (st + 1 - N) * X0pl;
748    //       };
749    //       //
750    //       if ( no18x && st == 18 && v == 0 ){
751    //      // skip plane 18x
752    //       } else {
753    //      if ( v == 1 && st == mask18b ){
754    //        // emulate plane 18x
755    //      } else {
756    //        if ( eplane[v][st] > 0. ){
757    //          th->Fill(xpos,eplane[v][st]);
758    //          if ( debug ) printf(" Filling: st %i v %i xpos %f energy %f \n",st,v,xpos,eplane[v][st]);
759    //        };
760    //      };
761    //       };
762    //       //
763    //     };
764      };
765      //
766      TF1 *lfit = new TF1("lfit",ccurve,0.,xmax,3);
767      if ( debug ) printf("qtot %f qtotparz %f \n",L2->GetCaloLevel2()->qtot,qtotparz);
768      E0 = qtotparz;
769      //  E0 = L2->GetCaloLevel2()->qtot;
770      a = 5.;
771      b = 0.5;
772      if ( debug ) printf(" STARTING PARAMETERS: E0 %f a %f b %f \n",E0,a,b);
773      lfit->SetParameters(E0,a,b);
774      //  lfit->SetParLimits(0,0.,1000.);
775      //  lfit->SetParLimits(1,-1.,80.);
776      //  lfit->SetParLimits(2,-1.,10.);
777      TString optio;
778      if ( debug ){
779        //    optio = "MERBOV";
780        //    optio = "MEROV";
781        //    optio = "EROV";
782        optio = "RNOV";
783        if ( draw ) optio = "ROV";
784      } else {
785        //    optio = "MERNOQ";
786        //    optio = "ERNOQ";
787        optio = "RNOQ";
788        if ( draw ) optio = "ROQ";
789      };
790      //
791      if ( debug ) printf(" OK, start the fitting procedure...\n");
792      //
793      fitresult = th->Fit("lfit",optio);
794      //
795      if ( debug ) printf(" the fit is done! result: %i \n",fitresult);
796      //
797      E0 = lfit->GetParameter(0);
798      a = lfit->GetParameter(1);
799      b = lfit->GetParameter(2);
800      errE0 = lfit->GetParError(0);
801      erra = lfit->GetParError(1);
802      errb = lfit->GetParError(2);
803      chi2 = lfit->GetChisquare();
804      ndf = lfit->GetNDF();
805      Float_t tmax = 0.;
806      if ( debug ) printf(" Parameters are retrieved \n");
807      if ( b != 0 ) tmax = (a - 1.)/b;
808      //
809      if ( fitresult != 0 ){
810        if ( debug ) printf(" The fit failed, no integrals calculation and asymm is set to -1. \n");
811        asymm = -1.;
812      } else {
813        Int_t npp = 1000;
814        double *xp=new double[npp];
815        double *wp=new double[npp];
816        lfit->CalcGaussLegendreSamplingPoints(npp,xp,wp,1e-12);
817        Float_t imax = lfit->IntegralFast(npp,xp,wp,0.,tmax);
818        //    Float_t imax = lfit->Integral(0.,tmax);
819        if ( debug ) printf(" Integral till maximum (%f): %f \n",tmax,imax);
820        Int_t np = 1000;
821        double *x=new double[np];
822        double *w=new double[np];
823        lfit->CalcGaussLegendreSamplingPoints(np,x,w,1e-12);
824        Float_t i10max = lfit->IntegralFast(np,x,w,0.,10.*tmax);
825        delete x;
826        delete w;
827        delete xp;
828        delete wp;
829        //    Float_t i10max = lfit->Integral(0.,10.*tmax);
830        if ( debug ) printf(" Integral: %f \n",i10max);
831        //
832        if ( i10max != imax ){
833          asymm = imax / (i10max-imax);
834        } else {
835          if ( debug ) printf(" i10max == imax, asymm undefined\n");
836          asymm = -2.;
837        };
838        if ( asymm != asymm ){
839          if ( debug ) printf(" asymm is nan \n");      
840          asymm = -3.;
841        };
842        //lfit->Integral(0.,tmax)/(lfit->Integral(0.,10.*tmax)-lfit->Integral(0.,tmax));
843        if ( debug ) printf(" Asymmetry has been calculated \n");
844      };
845      //
846      if ( asymm < 0. || ndf <= 0. || chi2 < 0. || tmax < 0. ){
847        if ( debug ) printf(" Funny asymm||ndf||chi2||tmax values, fit failed \n");
848        fitresult = 100;
849      };
850      //
851      if ( draw ){
852        //
853        tc->cd();    
854        //    gStyle->SetOptStat(11111);
855        tc->SetTitle();
856        th->SetTitle("");
857        th->SetName("");
858        th->SetMarkerStyle(20);
859        // axis titles
860        th->SetXTitle("Depth [X0]");
861        th->SetYTitle("Energy [MIP]");
862        th->DrawCopy("Perror");
863        lfit->Draw("same");
864        tc->Modified();
865        tc->Update();
866        //
867        gStyle->SetLabelSize(0);
868        gStyle->SetNdivisions(1,"XY");
869        //
870      } else {
871        if ( th ) th->Delete();
872      };
873      //
874      delete lfit;
875      //
876    };
877    
878    void CaloLong::Draw(){
879      //
880      Process();
881      Draw(-1);
882    };
883    
884    void CaloLong::Draw(Int_t view){
885      //
886      Int_t minv = 0;
887      Int_t maxv = 0;
888      //
889      if ( view == -1 ){
890        maxv = -1;
891      } else {
892        minv = view;
893        maxv = view+1;
894      };
895      //
896      Process();
897      //
898      gStyle->SetLabelSize(0.04);
899      gStyle->SetNdivisions(510,"XY");
900      //
901      if ( maxv != -1 ){
902        for (Int_t v=minv; v<maxv;v++){
903          TString hid = Form("clongv%i",v);
904          TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
905          if ( tc ){
906            //       tc->Clear();
907          } else {
908            tc = new TCanvas(hid,hid);
909          };
910          //
911          TString thid = Form("hlongv%i",v);        
912          TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
913          if ( th ) th->Delete();
914          //         th->Clear();
915          //         th->Reset();
916          //        } else {
917          th = new TH1F(thid,thid,22,-0.5,21.5);
918          //        };
919          tc->cd();
920          //
921          for (Int_t st=0;st<22;st++){
922            th->Fill(st,eplane[v][st]);
923          };
924          th->Draw();
925          tc->Modified();
926          tc->Update();
927        };
928      } else {
929        //
930        TString hid = Form("clongvyvx");    
931        TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
932        if ( tc ){
933        } else {
934          tc = new TCanvas(hid,hid);
935        };
936        //
937        TString thid = Form("hlongvyvx");  
938        TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
939        if ( th ) th->Delete();
940        th = new TH1F(thid,thid,44,-0.5,43.5);
941        tc->cd();
942        Int_t pp=0;
943        for (Int_t st=0;st<22;st++){
944          for (Int_t v=1; v>=0;v--){
945            //
946            th->Fill(pp,eplane[v][st]);
947            //
948            pp++;
949          };
950        };
951        th->Draw();
952        tc->Modified();
953        tc->Update();
954      };
955      //
956      gStyle->SetLabelSize(0);
957      gStyle->SetNdivisions(1,"XY");
958      //
959    };
960    
961    void CaloLong::Delete(){
962      Clear();
963      //delete this;
964    };
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.23