/[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.3 by mocchiut, Thu Dec 18 21:05:53 2008 UTC revision 1.8 by mocchiut, Tue Aug 11 14:23:09 2009 UTC
# Line 451  CaloLong::CaloLong(PamLevel2 *l2p){   Line 451  CaloLong::CaloLong(PamLevel2 *l2p){  
451    PKT = 0;    PKT = 0;
452    atime = 0;    atime = 0;
453    //    //
454      clp = L2->GetCaloLevel2();
455      //
456    sel = true;    sel = true;
457    cont = false;    cont = false;
458    N = 0;    N = 0;
# Line 459  CaloLong::CaloLong(PamLevel2 *l2p){   Line 461  CaloLong::CaloLong(PamLevel2 *l2p){  
461    //    //
462    no18x = true;    no18x = true;
463    debug = false;    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(){  void CaloLong::Clear(){
492    //    //
493    memset(eplane,0, 2*22*sizeof(Float_t));    memset(eplane,0, 2*22*sizeof(Float_t));
# Line 484  void CaloLong::Clear(){ Line 510  void CaloLong::Clear(){
510    
511  void CaloLong::Print(){  void CaloLong::Print(){
512    //    //
513    Process();    Fit();
514    //    //
515    printf("==================== Calorimeter Longitudinal Profile =======================\n");    printf("==================== Calorimeter Longitudinal Profile =======================\n");
516    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
# Line 597  void CaloLong::Process(){ Line 623  void CaloLong::Process(){
623    Int_t plane = 0;    Int_t plane = 0;
624    Int_t strip = 0;    Int_t strip = 0;
625    Float_t mip = 0.;    Float_t mip = 0.;
626      Bool_t gof = true;
627    for (Int_t i=0; i < L2->GetCaloLevel1()->istrip; i++){    for (Int_t i=0; i < L2->GetCaloLevel1()->istrip; i++){
628      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
629      eplane[view][plane] += mip;      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 (steal from Daniele's code)    // inclination factor (stolen from Daniele's code)
638    //    //
639    Float_t ytgx = 0;    Float_t ytgx = 0;
640    Float_t ytgy = 0;    Float_t ytgy = 0;
641    ytgx = 0.76 * L2->GetCaloLevel2()->tanx[0];    ytgx = 0.76 * clp->tanx[0];
642    ytgy = 0.76 * L2->GetCaloLevel2()->tany[0];      ytgy = 0.76 * clp->tany[0];  
643    X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) );    X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) );
644    //    //
645    // Find experimental plane of maximum    // Find experimental plane of maximum
# Line 635  void CaloLong::Process(){ Line 667  void CaloLong::Process(){
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){  Double_t ccurve(Double_t *ti,Double_t *par){
768    //    //
# Line 706  void CaloLong::Fit(Bool_t draw){ Line 834  void CaloLong::Fit(Bool_t draw){
834    //  th = new TH1F(thid,thid,int(NC*1.5),-0.2,xmax);    //  th = new TH1F(thid,thid,int(NC*1.5),-0.2,xmax);
835    th = new TH1F(thid,thid,100,-0.2,xmax);    th = new TH1F(thid,thid,100,-0.2,xmax);
836    //    //
837    for (Int_t st=N;st<(N+NC);st++){    // 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.;      enemip = 0.;
852      xpos = (st - N) * X0pl;      xpos = (st - mmin) * X0pl;
853      if ( st > N && st < (N+NC-1) ){            if ( st > mmin && st < mmax ){      
854        if ( no18x && ( st == 18+1 || st == mask18b+1 )){        if ( no18x && ( st == 18+1 || st == mask18b+1 )){
855          enemip = 2. * eplane[1][st];          if ( !maskYO ){
856              enemip = 2. * eplane[1][st];
857            } else {
858              enemip = eplane[1][st];
859            };
860        } else {        } else {
861          enemip = eplane[0][st-1] + eplane[1][st];          enemip = eplane[0][st-1] + eplane[1][st];
862        };        };
863      } else {      } else {
864        if ( st == N ) enemip = 2. * eplane[1][st];        if ( st == mmin ){
865        if ( st == (N+NC-1) ) enemip = 2. * eplane[0][st];          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. ){      if ( enemip > 0. ){
882        th->Fill(xpos,enemip);        th->Fill(xpos,enemip);
883        if ( debug ) printf(" Filling: st %i xpos %f energy %f \n",st,xpos,enemip);        if ( debug ) printf(" Filling: st %i xpos %f energy %f \n",st,xpos,enemip);
# Line 750  void CaloLong::Fit(Bool_t draw){ Line 908  void CaloLong::Fit(Bool_t draw){
908    };    };
909    //    //
910    TF1 *lfit = new TF1("lfit",ccurve,0.,xmax,3);    TF1 *lfit = new TF1("lfit",ccurve,0.,xmax,3);
911    E0 = L2->GetCaloLevel2()->qtot;    if ( debug ) printf("qtot %f qtotparz %f \n",clp->qtot,qtotparz);
912      E0 = qtotparz;
913      //  E0 = clp->qtot;
914    a = 5.;    a = 5.;
915    b = 0.5;    b = 0.5;
916    if ( debug ) printf(" STARTING PARAMETERS: E0 %f a %f b %f \n",E0,a,b);    if ( debug ) printf(" STARTING PARAMETERS: E0 %f a %f b %f \n",E0,a,b);
# Line 819  void CaloLong::Fit(Bool_t draw){ Line 979  void CaloLong::Fit(Bool_t draw){
979        if ( debug ) printf(" i10max == imax, asymm undefined\n");        if ( debug ) printf(" i10max == imax, asymm undefined\n");
980        asymm = -2.;        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));      //lfit->Integral(0.,tmax)/(lfit->Integral(0.,10.*tmax)-lfit->Integral(0.,tmax));
987      if ( debug ) printf(" Asymmetry has been calculated \n");      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 ){    if ( draw ){
996      //      //
997      tc->cd();          tc->cd();    
# Line 842  void CaloLong::Fit(Bool_t draw){ Line 1011  void CaloLong::Fit(Bool_t draw){
1011      gStyle->SetLabelSize(0);      gStyle->SetLabelSize(0);
1012      gStyle->SetNdivisions(1,"XY");      gStyle->SetNdivisions(1,"XY");
1013      //      //
1014      } else {
1015        if ( th ) th->Delete();
1016    };    };
1017    //    //
1018    delete lfit;    delete lfit;
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.23