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