/[PAMELA software]/calo/flight/CaloBragg/src/CaloBragg.cpp
ViewVC logotype

Diff of /calo/flight/CaloBragg/src/CaloBragg.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.8 by mocchiut, Fri Jan 28 10:40:44 2011 UTC revision 1.12 by mocchiut, Fri Jun 17 11:56:42 2011 UTC
# Line 26  CaloBragg::CaloBragg(PamLevel2 *l2p){   Line 26  CaloBragg::CaloBragg(PamLevel2 *l2p){  
26    debug = false;    debug = false;
27    usetrack = false;    usetrack = false;
28    usepl18x = false;    usepl18x = false;
29      newchi2 = false;
30      usenewBB = false;
31      fzeta = -1.;
32    //    //
33  };  };
34    
35  void CaloBragg::Clear(){  void CaloBragg::Clear(){
36    //    //
37      ndf = 0;
38    tr = 0;    tr = 0;
39    sntr = 0;    sntr = 0;
40    //   qtchi2 = 0.;    //   qtchi2 = 0.;
# Line 42  void CaloBragg::Clear(){ Line 46  void CaloBragg::Clear(){
46    lpetot = 0.;    lpetot = 0.;
47    lppskip = 0.;    lppskip = 0.;
48    memset(calorimetro,0,44*2*sizeof(Float_t));    memset(calorimetro,0,44*2*sizeof(Float_t));
49    memset(spessore,0,3*sizeof(Float_t));    memset(spessore,0,4*sizeof(Float_t));
50    memset(estremi,0,2*2*sizeof(Float_t));    memset(estremi,0,2*2*sizeof(Float_t));
51    Integrale=0.;    Integrale=0.;
52    
# Line 69  void CaloBragg::Print(){ Line 73  void CaloBragg::Print(){
73    printf(" Z from loop %f: \n", lpz);    printf(" Z from loop %f: \n", lpz);
74    printf(" energy from loop %f: \n", lpetot);    printf(" energy from loop %f: \n", lpetot);
75    printf(" plane not used for loop %f: \n", lppskip);    printf(" plane not used for loop %f: \n", lppskip);
76      printf(" ndf: %i \n",ndf);
77    printf("========================================================================\n");    printf("========================================================================\n");
78    //    //
79  };  };
# Line 301  void CaloBragg::Process(Int_t ntr){ Line 306  void CaloBragg::Process(Int_t ntr){
306    };    };
307    //  };    //  };
308    
309      if ( startZero ) {
310        estremi[0][0] = 0.;
311        //    estremi[0][1] = 0.;
312      }
313        
314    /*integrale: energia totale rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/    /*integrale: energia totale rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/
315    for(Int_t pl=0; pl<(2*NPLA); pl++){    for(Int_t pl=0; pl<(2*NPLA); pl++){
# Line 320  void CaloBragg::Process(Int_t ntr){ Line 328  void CaloBragg::Process(Int_t ntr){
328    //  mediatroncata();  // out:  1)chi2, 2)z, 3)Etot, 4)Pskip    //  mediatroncata();  // out:  1)chi2, 2)z, 3)Etot, 4)Pskip
329        
330    /*z ed energia con loop*/    /*z ed energia con loop*/
331      if ( debug ) printf(" call Zdaloop with integrale %f \n",Integrale);
332    Zdaloop(); // out:  1)chi2, 2)z, 3)Etot, 4)Pskip    Zdaloop(); // out:  1)chi2, 2)z, 3)Etot, 4)Pskip
333    
334        
# Line 329  void CaloBragg::Process(Int_t ntr){ Line 338  void CaloBragg::Process(Int_t ntr){
338  };  };
339    
340    
341    Float_t CaloBragg::Integral(){
342      Process();
343    
344      Float_t dEpianiloop[44];
345      Int_t tz1=(Int_t)lpz;
346      Enetrack(&tz1, &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop
347    
348    
349      Float_t integ = 0.;
350      for(Int_t i=0;i<=estremi[1][0];i++){
351        //    integ += dEplan[i];
352        //printf(" step %i integ %f deplan %f \n",i,integ,dEplan[i]);
353        integ += dEpianiloop[i];
354        //    printf(" step %i integ %f deplan %f \n",i,integ,dEpianiloop[i]);
355      }
356      return integ;
357    }
358    
359    Float_t CaloBragg::LastIntegral(){
360      Process();
361    
362      Float_t integ = 0.;
363      for(Int_t i=0;i<=estremi[1][0];i++){
364        integ += dEplan[i];
365        //printf(" step %i integ %f deplan %f \n",i,integ,dEplan[i]);
366      }
367      return integ;
368    }
369    
370    
371  void CaloBragg::Draw(){  void CaloBragg::Draw(){
372    
373    Process();    Process();
374    
375      this->Draw(0.,0.);
376    
377    }
378    
379    void CaloBragg::Draw(Int_t Z, Float_t enetot){
380    
381    //  Float_t dEpianimean[44];    //  Float_t dEpianimean[44];
382    Float_t dEpianiloop[44];    Float_t dEpianiloop[44];
383    Float_t Depth[44];    Float_t Depth[44];
384    //  Int_t tz=(Int_t)qtz;    //  Int_t tz=(Int_t)qtz;
385    Int_t tz1=(Int_t)lpz;    Int_t tz1= Z;
386    //  Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata    Float_t enet = enetot;
387    Enetrack(&tz1, &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop    //  Float_t enet = lpetot;
388      
389      if ( Z > 0. && enetot > 0. ){
390        estremi[0][0] = 0;
391        estremi[1][0] = 43;
392    
393    
394        Float_t    ytgx = 0.;
395        Float_t    ytgy = 0.;
396        
397        //lunghezza effettiva di silicio attraversata (mm)
398        Float_t    SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy));
399        
400        spessore[0] = (SiCross/10.) * rhoSi; //spessore silicio in g/cm2
401        
402        /*tungsteno*/
403        
404        //rapporto tra rilasci energetici nei due materiali
405        Float_t WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad
406        //gcm2W = WCross/10. * rhoW;
407        
408        //       (g/cm2W)/(g/cm2Si)
409        spessore[3] =  (WCross/10.) * rhoW;
410        Float_t a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664);  //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664);
411        spessore[1] =  a;
412        //riscala mip allo spessore attraversato
413        spessore[2] = MIP*(SiCross/ySi);  
414    
415      } else {
416        tz1=(Int_t)lpz;
417        enet = lpetot;
418        //  Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata
419      
420      }
421      Enetrack(&tz1, &enet, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop
422    
423    Float_t sp= spessore[0]*spessore[1];    Float_t sp= spessore[0]*spessore[1];
424    for(Int_t i=0;i<44;i++)Depth[i]=i*sp;    for(Int_t i=0;i<44;i++)Depth[i]=i*sp;
# Line 392  void CaloBragg::Draw(){ Line 471  void CaloBragg::Draw(){
471    //  tc->cd(2);    //  tc->cd(2);
472    tc->cd();    tc->cd();
473    //    //
474    for(Int_t i=0;i<=estremi[1][0];i++)th3->Fill(Depth[i],dEpianiloop[i]);    for(Int_t i=0;i<=estremi[1][0];i++){
475        th3->Fill(Depth[i],dEpianiloop[i]);
476        //    printf(" i %i Depth %f depianiloop %f \n",i,Depth[i],dEpianiloop[i]);
477      }
478    th3->Draw();    th3->Draw();
479    th2->Draw("same");    th2->Draw("same");
480    
# Line 411  void CaloBragg::LoadParam(){ Line 493  void CaloBragg::LoadParam(){
493    
494    //    //
495    elem[0] = 1.00794; //H  1    elem[0] = 1.00794; //H  1
496    elem[1] = 4.0026;  //He 2    elem[1] = 4.002602;  //He 2
497    elem[2] = 6.941;   //Li 3    elem[2] = 6.941;   //Li 3
498    elem[3] = 9.012182;//Be 4    elem[3] = 9.012182;//Be 4
499    elem[4] = 10.811;  //B  5    elem[4] = 10.811;  //B  5
500    elem[5] = 12.0107; //C  6    elem[5] = 12.0107; //C  6
501    elem[6] = 14.00674;//N  7    elem[6] = 14.00674;//N  7
502    elem[7] = 15.9994; //O  8    elem[7] = 15.9994; //O  8
503    elem[8] = 18.9984; //F  9    elem[8] = 18.9984032; //F  9
504    elem[9] = 20.1797; //Ne 10    elem[9] = 20.1797; //Ne 10
505    elem[10] = 22.98977;//Na 11    elem[10] = 22.98977;//Na 11
506    elem[11] = 24.3050; //Mg 12    elem[11] = 24.3050; //Mg 12
# Line 465  void CaloBragg::LoadParam(){ Line 547  void CaloBragg::LoadParam(){
547    pigr = 3.1415;    pigr = 3.1415;
548    Na = 6.02e-23;    Na = 6.02e-23;
549    ZA = 0.49; /*Z/A per Si*/    ZA = 0.49; /*Z/A per Si*/
550    ISi =182e-06; /*MeV*/    //  ISi =182e-06; /*MeV*/
551      ISi = 171e-06; /*MeV*/
552      IW  = 735e-06; /*MeV*/
553      //  ISi =0.0001059994; /*GeV!!*/ no era giusto!!
554    Me = 0.511; /* MeV*/    Me = 0.511; /* MeV*/
555    MassP = 931.27;/*MeV*/    MassP = 931.27;/*MeV*/
556    r2 = 7.95e-26; /*ro*ro in cm */    r2 = 7.95e-26; /*ro*ro in cm */
# Line 500  void CaloBragg::conversione(){ Line 585  void CaloBragg::conversione(){
585    //lunghezza effettiva di silicio attraversata (mm)    //lunghezza effettiva di silicio attraversata (mm)
586    SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy));    SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy));
587        
588    spessore[0] = SiCross/10. * rhoSi; //spessore silicio in g/cm2    spessore[0] = (SiCross/10.) * rhoSi; //spessore silicio in g/cm2
589        
590    /*tungsteno*/    /*tungsteno*/
591    ytgx = yW * L2->GetCaloLevel2()->tanx[0];    ytgx = yW * L2->GetCaloLevel2()->tanx[0];
# Line 510  void CaloBragg::conversione(){ Line 595  void CaloBragg::conversione(){
595    WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad    WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad
596    //gcm2W = WCross/10. * rhoW;    //gcm2W = WCross/10. * rhoW;
597    
   a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664);  //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664);  
   
598    //       (g/cm2W)/(g/cm2Si)    //       (g/cm2W)/(g/cm2Si)
599      spessore[3] =  (WCross/10.) * rhoW;
600      a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664);  //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664);
601    spessore[1] =  a;    spessore[1] =  a;
   
602    //riscala mip allo spessore attraversato    //riscala mip allo spessore attraversato
603    spessore[2] = MIP*(SiCross/ySi);      spessore[2] = MIP*(SiCross/ySi);  
   
604  };//end conversione  };//end conversione
605    
606    
607    
608    
609    
610  void CaloBragg::BetheBloch(Float_t *x, Float_t *z, Float_t *Mass, Float_t *gam, Float_t *Bet, Float_t *out){  void CaloBragg::BetheBloch(Float_t *x, Float_t *z, Float_t *Mass, Float_t *gam, Float_t *Bet, Float_t *out, Float_t II){
611    
612    //rilascio energetico con bethe bloch con correzioni    //rilascio energetico con bethe bloch con correzioni
613    //in:    x: g/cm2    //in:    x: g/cm2
# Line 542  void CaloBragg::BetheBloch(Float_t *x, F Line 625  void CaloBragg::BetheBloch(Float_t *x, F
625    Float_t lg =0.;    Float_t lg =0.;
626    Float_t Energia=0.;    Float_t Energia=0.;
627    Float_t C=0.;    Float_t C=0.;
628      Float_t INo = ISi;
629    
630      if ( usenewBB ) INo = II;
631    
632    eta = (*gam)*(*Bet);    eta = (*gam)*(*Bet);
633    
634    //Bet=3/gam;  SQ(*gam) * SQ(*Bet)    //Bet=3/gam;  SQ(*gam) * SQ(*Bet)
635    Wmax = 2.* Me * SQ(eta) / (1. + 2.*(*gam)*Me/(*Mass) + SQ(Me)/SQ(*Mass));    Wmax = 2.* Me * SQ(eta) / (1. + 2.*(*gam)*Me/(*Mass) + SQ(Me)/SQ(*Mass));
636        
637    lg = 2.* Me * SQ(eta) * Wmax / SQ(ISi);    lg = 2.* Me * SQ(eta) * Wmax / SQ(INo);
638    //  Energia = x* 2 * pigr * Na * r2 * Me * rhoSi *ZA*  SQ(z)/SQ(Bet) * lg;    //  Energia = x* 2 * pigr * Na * r2 * Me * rhoSi *ZA*  SQ(z)/SQ(Bet) * lg;
639    C=(0.42237*pow(eta,-2.) + 0.0304*pow(eta,-4.) - 0.00038*pow(eta,-6.))*pow(10.,-6.)* pow(ISi,2.) +    C=(0.42237*pow(eta,-2.) + 0.0304*pow(eta,-4.) - 0.00038*pow(eta,-6.))*pow(10.,-6.)* pow(INo,2.) +
640      (3.858*pow(eta,-2.) - 0.1668*pow(eta,-4.) + 0.00158*pow(eta,-6.))*pow(10.,-9.)*pow(ISi,3.);      (3.858*pow(eta,-2.) - 0.1668*pow(eta,-4.) + 0.00158*pow(eta,-6.))*pow(10.,-9.)*pow(INo,3.);
641        
642    if(eta <= 0.13) C= C * log(eta/0.0653)/log(0.13/0.0653);          if(eta <= 0.13) C= C * log(eta/0.0653)/log(0.13/0.0653);      
643    
# Line 564  void CaloBragg::BetheBloch(Float_t *x, F Line 650  void CaloBragg::BetheBloch(Float_t *x, F
650    
651    
652    
653  void CaloBragg::ELOSS(Float_t *dx, Int_t *Z, Float_t *Etot, Float_t *out){  void CaloBragg::ELOSS(Float_t *dx, Int_t *Z, Float_t *Etot, Float_t *out, Float_t II){
654        
655    /*perdita di energia per ioni pesanti (come da routine geant)*/    /*perdita di energia per ioni pesanti (come da routine geant)*/
656    //  in : dx    => spessore g/cm2    //  in : dx    => spessore g/cm2
# Line 586  void CaloBragg::ELOSS(Float_t *dx, Int_t Line 672  void CaloBragg::ELOSS(Float_t *dx, Int_t
672    
673    Bet = sqrt((SQ(gam) -1.)/SQ(gam));    Bet = sqrt((SQ(gam) -1.)/SQ(gam));
674        
675    v= 121.4139*(Bet/pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet/pow((*Z),(2./3.))));    //  v= 121.4139*(Bet/pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet/pow((*Z),(2./3.))));
676      v= 121.4139*(Bet*pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet*pow((*Z),(2./3.)))); // EMI AAAAGGH!!
677    
678    //carica effettiva    //carica effettiva
679    Q= (*Z)*(1- (1.034 - 0.1777*exp(-0.08114*(*Z)))*exp(-v));    Q= (*Z)*(1- (1.034 - 0.1777*exp(-0.08114*(*Z)))*exp(-v));
680    
681    //perdita energia per un protone    //perdita energia per un protone
682    Float_t protone =1.;    Float_t protone =1.;
683    Float_t Mass=(elem[*Z-1]*MassP);    //  Float_t Mass=(elem[*Z-1]*MassP); //EMI
684    BetheBloch(dx, &protone, &Mass, &gam, &Bet, &dEP);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP);    //  BetheBloch(dx, &protone, &Mass, &gam, &Bet, &dEP);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP);
685    
686      BetheBloch(dx, &protone, &MassP, &gam, &Bet, &dEP, II);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP); //EMI
687        
688    *out= (SQ(Q)*(dEP));//*dx;    *out= (SQ(Q)*(dEP));//*dx;
689    
# Line 627  void CaloBragg::Enetrack(Int_t* Z, Float Line 716  void CaloBragg::Enetrack(Int_t* Z, Float
716    for( Int_t ipla=((int)(*primo)); ipla<= ((int)(*ultimo)); ipla++){    for( Int_t ipla=((int)(*primo)); ipla<= ((int)(*ultimo)); ipla++){
717      dE=0.;      dE=0.;
718      //spessore silicio corretto x inclinazione, z, energia, out:rilascio      //spessore silicio corretto x inclinazione, z, energia, out:rilascio
719      ELOSS(&spessore[0], Z, &Ezero, &dE);//spessore in g/cm2!!      ELOSS(&spessore[0], Z, &Ezero, &dE, ISi);//spessore in g/cm2!!
720      if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop      if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop
721        out[ipla] = Ezero - Massa; //MeV        out[ipla] = Ezero - Massa; //MeV
722        return;        return;
# Line 635  void CaloBragg::Enetrack(Int_t* Z, Float Line 724  void CaloBragg::Enetrack(Int_t* Z, Float
724      }else{      }else{
725        out[ipla] = dE; //MeV        out[ipla] = dE; //MeV
726        Ezero = Ezero - dE;//energia residua        Ezero = Ezero - dE;//energia residua
727          if ( debug ) printf(" zompa %i out %f dE %f ezero %f \n",ipla,out[ipla],dE,Ezero);
728      };      };
729      //se sono su un piano Y (tutti i pari) dopo c'e' il tungsteno      //se sono su un piano Y (tutti i pari) dopo c'e' il tungsteno
730      if(ipla%2 == 0){      if(ipla%2 == 0){
731        /*tungsteno*/        /*tungsteno*/
732        dE=0.;        dE=0.;
733        Float_t sp= spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato  in g/cm2        Float_t sp = 0.;
734        ELOSS(&sp, Z, &Ezero, &dE);        Float_t II = ISi;
735          if ( usenewBB ){
736            sp = spessore[3];        
737            II = IW;
738          } else {
739            sp = spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato  in g/cm2
740          }
741          //      printf(" sp %f II %f \n",sp,II);
742          ELOSS(&sp, Z, &Ezero, &dE,II);
743        if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop        if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop
744          return;          return;
745        }else{        }else{
# Line 671  void CaloBragg::chiquadro(Float_t dE[], Line 769  void CaloBragg::chiquadro(Float_t dE[],
769    Float_t badplane=0.;    Float_t badplane=0.;
770    Float_t badplanetot=0.;    Float_t badplanetot=0.;
771    Float_t w,wi;    Float_t w,wi;
772      //
773    for(Int_t ipla=0; ipla<2*NPLA; ipla++){    if ( newchi2 ){
774      //tutti i piani attraversati dalla traiettoria      ndf = 0;
775      if(calorimetro[ipla][0] != -1.){ //      sum = 0.;
776        w=0.; //normalizzazione;      for( Int_t ipla=((int)(estremi[0][0])); ipla<= ((int)(estremi[1][0])); ipla++){
777        wi=1.;//peso        sum += pow((dE[ipla] - (calorimetro[ipla][1] * spessore[2]))/(0.05*dE[ipla]),2.);
778          //      printf(" quiqui: dE %f calor %f spessore[2] %f \n",dE[ipla],spessore[2]*calorimetro[ipla][1],spessore[2]);
779          ndf++;
780        }
781        ndf -= 2;
782        if ( ndf > 0 ) sum /= (float)ndf;
783        out[0] = sum;
784        out[1] = 0.;
785        out[2] = (int)(estremi[1][0])-ndf;
786        //    printf(" sum %f ndf %i \n ",sum,ndf);
787      } else {
788        for(Int_t ipla=0; ipla<2*NPLA; ipla++){
789          //tutti i piani attraversati dalla traiettoria
790          if(calorimetro[ipla][0] != -1.){ //
791            w=0.; //normalizzazione;
792            wi=1.;//peso
793    
794        //tolgo piani attraversati dalla traccia ma precedenti il piano individuato come ingresso          //tolgo piani attraversati dalla traccia ma precedenti il piano individuato come ingresso
795        if (ipla<estremi[0][0])  wi=0.;          if (ipla<estremi[0][0])  wi=0.;
796                
797        //tolgo piani attraversati da traccia ma successivi all'ultimo se sono diversi da 0          //tolgo piani attraversati da traccia ma successivi all'ultimo se sono diversi da 0
798        //if((ipla>estremi[1][0]) && (calorimetro[ipla][1] >0.) ) wi=0.;          //if((ipla>estremi[1][0]) && (calorimetro[ipla][1] >0.) ) wi=0.;
799        if((ipla>estremi[1][0])) wi=0.;          if((ipla>estremi[1][0])) wi=0.;
800    
801        //normalizzazione          //normalizzazione      
802        if (calorimetro[ipla][1] != 0.)  w=1./(calorimetro[ipla][1]* MIP);    //            if (calorimetro[ipla][1] != 0.)  w=1./(calorimetro[ipla][1]* MIP);    //        
803                
804        //tolgo piani con rilasci inferiori al 30% del precedente          //tolgo piani con rilasci inferiori al 30% del precedente
805        if(calorimetro[ipla][1] < (0.7*PianoPrecedente)){ // cosi' i piani senza rilascio non vengono considerati nel calcolo del chi2          if(calorimetro[ipla][1] < (0.7*PianoPrecedente)){ // cosi' i piani senza rilascio non vengono considerati nel calcolo del chi2
806          wi=0.;            wi=0.;
807          //se sono piani intermedi (non si e' fermta) li considero non buoni            //se sono piani intermedi (non si e' fermta) li considero non buoni
808          if( (ipla <= estremi[1][0]) && (calorimetro[ipla][1] !=0.)){//            if( (ipla <= estremi[1][0]) && (calorimetro[ipla][1] !=0.)){//
809            badplane+=1.;              badplane+=1.;
810            badplanetot+=1.;              badplanetot+=1.;
811          };            };
812        };          };
813        
814        //meno peso ai piani con rilasci maggiori di 1000 MIP          //meno peso ai piani con rilasci maggiori di 1000 MIP
815        //      if(calorimetro[ipla][1] > 1000) wi=0.5;          //      if(calorimetro[ipla][1] > 1000) wi=0.5;
816        if(calorimetro[ipla][1] > 1200.) wi=0.5;          if(calorimetro[ipla][1] > 1200.) wi=0.5;
817                            
818        Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));          Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));
819    
820        sum += SQ(arg); // w*wi*(dEpiani[p][v]-(eplane[p][v]*MIP))));//( dEpiani[p][v] - (eplane[p][v]*MIP));          sum += SQ(arg); // w*wi*(dEpiani[p][v]-(eplane[p][v]*MIP))));//( dEpiani[p][v] - (eplane[p][v]*MIP));
821        if(debug){          if(debug){
822          printf("dedx  calcolata  %f e reale  %f  \n",dE[ipla],(calorimetro[ipla][1] * MIP));            printf("dedx  calcolata  %f e reale  %f  \n",dE[ipla],(calorimetro[ipla][1] * MIP));
823        }          }
824        //se trovo piano non buono (tolto quindi wi=0) non modifico il piano precedente          //se trovo piano non buono (tolto quindi wi=0) non modifico il piano precedente
825        if(wi != 0.){//          if(wi != 0.){//
826          PianoPrecedente= calorimetro[ipla][1];//tengo piano precedente            PianoPrecedente= calorimetro[ipla][1];//tengo piano precedente
827          badplane = 0.;//azzero contatore piani scartati consecutivi            badplane = 0.;//azzero contatore piani scartati consecutivi
828            };
829        };        };
     };  
830                        
831      //da Emi        //da Emi
832      if(badplane > 2){              if(badplane > 2){      
833        //      printf(" AAAAAAAAAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG\n");          //      printf(" AAAAAAAAAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG\n");
834        out[1] =79.;          out[1] =79.;
835        break;          break;
836      };        };
     
   };//fine loop piani  
   //chi2,frammentato,pskip    
   out[0]=sum;  
   out[2]=badplanetot;  
837        
838        };//fine loop piani
839        //chi2,frammentato,pskip  
840        out[0]=sum;
841        out[2]=badplanetot;
842      }  
843  };//end chiquadro  };//end chiquadro
844    
845    
# Line 743  void CaloBragg::loopze( Float_t step, Fl Line 856  void CaloBragg::loopze( Float_t step, Fl
856    //    //
857    
858    
   Float_t dEplan[2*NPLA];//energia rilasciata calcolata  
859    memset(dEplan,0,2*NPLA*sizeof(Float_t));    memset(dEplan,0,2*NPLA*sizeof(Float_t));
860    
861    Int_t Z = 0;// z iniziale    Int_t Z = 0;// z iniziale
# Line 755  void CaloBragg::loopze( Float_t step, Fl Line 867  void CaloBragg::loopze( Float_t step, Fl
867    Float_t energia =0.;//energia del loop    Float_t energia =0.;//energia del loop
868    
869    Float_t chi2[3] = {0,0,0};//out dal calcolo chi2: chi2, piani consecutivi saltati, piani totali saltati    Float_t chi2[3] = {0,0,0};//out dal calcolo chi2: chi2, piani consecutivi saltati, piani totali saltati
870      
871      Int_t zmin = (int)Zstart;  
872    Int_t max=32;//max z di cui so la massa :P    Int_t max=32;//max z di cui so la massa :P
873    if((Zlimite)<=31) max=(int)(Zlimite) + 1;    if((Zlimite)<=31) max=(int)(Zlimite) + 1;
874    
875      if ( fzeta > 0. ){
876        zmin = fzeta;
877        max = fzeta+1;
878      }
879        
880    Int_t colmax=32;    Int_t colmax=32;
881    Int_t rowmax=3000;    Int_t rowmax=3000;
# Line 769  void CaloBragg::loopze( Float_t step, Fl Line 887  void CaloBragg::loopze( Float_t step, Fl
887    Int_t imax = nostep/2;    Int_t imax = nostep/2;
888    
889    //loop elementi      //loop elementi  
890    for(Int_t inucl=(int)(Zstart); inucl<max; inucl++){    for(Int_t inucl=zmin; inucl<max; inucl++){
891            
892      Z= inucl;      Z= inucl;
893        
# Line 784  void CaloBragg::loopze( Float_t step, Fl Line 902  void CaloBragg::loopze( Float_t step, Fl
902        energia=  Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica..        energia=  Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica..
903    
904        Enetrack(&Z, &energia, &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani        Enetrack(&Z, &energia, &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani
         
905        //calcolo chi2        //calcolo chi2
906        chiquadro(dEplan,chi2);        chiquadro(dEplan,chi2);
907    
908          //      printf(" last deplan from: Z = %i iene %i energia %f chi2 %f \n",inucl,iene,energia,chi2[0]);
909                
910        if( (chi2[1] != 79.) ){//salto quelli che frammentano        if( (chi2[1] != 79.) ){//salto quelli che frammentano
911          matrixchi2[inucl][iene2][0]=chi2[0];//valore chi2 per questo z a questa energia          matrixchi2[inucl][iene2][0]=chi2[0];//valore chi2 per questo z a questa energia
# Line 805  void CaloBragg::loopze( Float_t step, Fl Line 924  void CaloBragg::loopze( Float_t step, Fl
924    
925    
926    //Emi    //Emi
927    for (Int_t nu=(int)(Zstart); nu<max; nu++){    for (Int_t nu=zmin; nu<max; nu++){
928      for (Int_t en=0; en<nostep; en++){      for (Int_t en=0; en<nostep; en++){
929        if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] >0.)){        if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] >0.)){
930          bestchi2[0]= matrixchi2[nu][en][0];// chi2          bestchi2[0]= matrixchi2[nu][en][0];// chi2
# Line 949  void CaloBragg::Zdaloop(){ Line 1068  void CaloBragg::Zdaloop(){
1068    //------------primo loop   ----------------------    //------------primo loop   ----------------------
1069    //     energia   ezero, zstart  zstop    //     energia   ezero, zstart  zstop
1070    //  loopze(Integrale,zero,zmin,zmax);    //  loopze(Integrale,zero,zmin,zmax);
1071    loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,50);  
1072      //->  loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,50);
1073      loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,200);
1074    
1075    //  loopze(Integrale*2.,Integrale/100.,zmin,zmax);    //  loopze(Integrale*2.,Integrale/100.,zmin,zmax);
1076    //  printf(" Integrale %f , outene %f \n",Integrale,bestchi2[2]);    if ( debug ) printf(" Integrale %f , outene %f \n",Integrale,bestchi2[2]);
1077        
1078    //------------secondo loop  ----------------------    //------------secondo loop  ----------------------
1079    for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i];    for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i];
# Line 966  void CaloBragg::Zdaloop(){ Line 1088  void CaloBragg::Zdaloop(){
1088    zmin=bestchitemp[1]-1;    zmin=bestchitemp[1]-1;
1089    zmax=bestchitemp[1]+1;    zmax=bestchitemp[1]+1;
1090    //  loopze(step,zero,zmin,zmax); //    //  loopze(step,zero,zmin,zmax); //
1091    loopze(step,step/2.,zmin,zmax,200); //  
1092      //->  loopze(step,step/2.,zmin,zmax,200); //
1093      loopze(step,step/2.,zmin,zmax,500); //
1094    
1095      if ( debug ) printf(" Integrale2 %f , outene %f step %f \n",Integrale,bestchi2[2],step);
1096    
1097        
1098    //chi2,z,Etot,Pskip    //chi2,z,Etot,Pskip
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.23