/[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.3 by mocchiut, Mon Sep 22 20:01:19 2008 UTC revision 1.7 by mocchiut, Mon Dec 14 14:37:58 2009 UTC
# Line 3  Line 3 
3    
4  ClassImp(CaloBragg);  ClassImp(CaloBragg);
5  //--------------------------------------  //--------------------------------------
6  /**  /*
7   * Default constructor   * Default constructor
8   */   */
9  CaloBragg::CaloBragg(){  CaloBragg::CaloBragg(){
# Line 25  CaloBragg::CaloBragg(PamLevel2 *l2p){   Line 25  CaloBragg::CaloBragg(PamLevel2 *l2p){  
25    //    //
26    debug = false;    debug = false;
27    usetrack = false;    usetrack = false;
28      usepl18x = false;
29    //    //
30  };  };
31    
# Line 32  void CaloBragg::Clear(){ Line 33  void CaloBragg::Clear(){
33    //    //
34    tr = 0;    tr = 0;
35    sntr = 0;    sntr = 0;
36    qtchi2 = 0.;  //   qtchi2 = 0.;
37    qtz = 0.;  //   qtz = 0.;
38    qtetot = 0.;  //   qtetot = 0.;
39    qtpskip = 0.;  //   qtpskip = 0.;
40    lpchi2 = 0.;    lpchi2 = 0.;
41    lpz = 0.;    lpz = 0.;
42    lpetot = 0.;    lpetot = 0.;
# Line 44  void CaloBragg::Clear(){ Line 45  void CaloBragg::Clear(){
45    memset(spessore,0,3*sizeof(Float_t));    memset(spessore,0,3*sizeof(Float_t));
46    memset(estremi,0,2*2*sizeof(Float_t));    memset(estremi,0,2*2*sizeof(Float_t));
47    Integrale=0.;    Integrale=0.;
48    
49      for(Int_t l=0;l<44;l++){
50        calorimetro[l][0]=-1.;
51      }
52    
   //  
53  };  };
54    
55  void CaloBragg::Print(){  void CaloBragg::Print(){
# Line 55  void CaloBragg::Print(){ Line 59  void CaloBragg::Print(){
59    //    //
60    printf("========================================================================\n");    printf("========================================================================\n");
61    printf(" OBT: %u PKT: %u ATIME: %u Track %i Use track %i \n",OBT,PKT,atime,tr,usetrack);    printf(" OBT: %u PKT: %u ATIME: %u Track %i Use track %i \n",OBT,PKT,atime,tr,usetrack);
62    printf(" chi 2 from truncated mean: %f \n", qtchi2);    printf(" first plane: %f \n", estremi[0][0]);
63    printf(" Z from truncated mean %f: \n", qtz);    printf(" last plane: %f \n", estremi[1][0]);
64    printf(" energy from truncated mean %f: \n", qtetot);  //   printf(" chi 2 from truncated mean: %f \n", qtchi2);
65    printf(" plane not used for truncated mean %f: \n", qtpskip);  //   printf(" Z from truncated mean %f: \n", qtz);
66    //   printf(" energy from truncated mean %f: \n", qtetot);
67    //   printf(" plane not used for truncated mean %f: \n", qtpskip);
68    printf(" chi 2 from loop %f:  \n", lpchi2);    printf(" chi 2 from loop %f:  \n", lpchi2);
69    printf(" Z from loop %f: \n", lpz);    printf(" Z from loop %f: \n", lpz);
70    printf(" energy from loop %f: \n", lpetot);    printf(" energy from loop %f: \n", lpetot);
# Line 119  void CaloBragg::Process(Int_t ntr){ Line 125  void CaloBragg::Process(Int_t ntr){
125    for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){    for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
126      //      //
127      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
128        //
129        if ( !usepl18x && view==0 && plane==18 ) mip = 0.;
130        //
131      epiano[plane][view]+=mip;      epiano[plane][view]+=mip;
132      //      //
133      //      //
# Line 162  void CaloBragg::Process(Int_t ntr){ Line 171  void CaloBragg::Process(Int_t ntr){
171    
172    /*per ogni evento calcolo la conversione mip e w attraversato in equivalente Si*/    /*per ogni evento calcolo la conversione mip e w attraversato in equivalente Si*/
173    conversione(); // out: 1) g/cm2 Si , 2) spessoreW equivalente in Si, 3)Mip corretta per inclinazione    conversione(); // out: 1) g/cm2 Si , 2) spessoreW equivalente in Si, 3)Mip corretta per inclinazione
   //cout<<"spessore= "<<spessore[0]<<" gcm2si; "<<spessore[1]<<"W in Si; "<<spessore[2]<<"mip "<<endl;  
174    
175      /*settaggio della soglia per il loop sulla determinazione del piano di partenza */
176      Float_t ordplane[44];//mi serve per la media troncata
177      memset(ordplane,0,44*sizeof(Float_t));
178    
179      for(Int_t ipla=0; ipla< 2*NPLA; ipla++)  ordplane[ipla]=calorimetro[ipla][1]; //energia del piano
180    
181      
182      //ordino tutte le energie dei piani in ordine crescente
183    
184        Long64_t work[200];
185        Int_t ind = 0;
186        //Int_t l = 0;
187        Int_t RN = 0;
188        Float_t sum4 = 0.;
189        Float_t qm = 0.;
190        while ( RN < 4 && ind < 44 ){
191          qm = TMath::KOrdStat((Long64_t)44,ordplane,(Long64_t)ind,work);
192          if (qm >= 0.7 ){  
193            if ( RN < 4 ){
194              sum4 += qm;
195              RN++;
196            };
197          };
198          ind++;
199        };
200        //
201        //sum4 /= (Float_t)RN;
202        Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2])));
203        if(Zmean ==0.) Zmean=1.;
204        if ( Zmean < 1. ) Zmean = 1.;
205    
206    
207    /*trova primo e ultimo piano attraversati*/    /*trova primo e ultimo piano attraversati*/
208    Int_t p = 0;//contatore piani    Int_t p = 0;//contatore piani
209    //per il primo parte da 0 e va in giù    //per il primo parte da 0 e va in giu'
210    while( estremi[0][1] == 0  &&  p<(2*NPLA) ){    while( estremi[0][1] <= 0.  &&  p<(2*NPLA) ){ // era ==0 ma ricorda i problemi con Float == !!!!!
     //cout<<"piano "<<p <<"; strip "<<calorimetro[p][0]<<"; en "<<calorimetro[p][1]<<endl;  
211      //    if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){      //    if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){
212      if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){ //0.7 soglia minima      //   if( (calorimetro[p][0] >0) && (calorimetro[p][1]*MIP >0.3)){ //0.7 mip = 70MeV soglia minima
213        if( (calorimetro[p][0] >0) && (calorimetro[p][1]*MIP >Zmean*0.7)){ // 70% della MIP
214        estremi[0][0]=p;        estremi[0][0]=p;
215        estremi[0][1]=calorimetro[p][1] *MIP; //energia in MeV        estremi[0][1]=calorimetro[p][1] *MIP; //energia in MeV
216      };      };
# Line 180  void CaloBragg::Process(Int_t ntr){ Line 218  void CaloBragg::Process(Int_t ntr){
218    };    };
219   //ultimo parte da 44 e sale   //ultimo parte da 44 e sale
220    p=43;    p=43;
221    while( (estremi[1][1] == 0)  &&  (p>(int)estremi[0][0]) ){    while( (estremi[1][1] <= 0.)  &&  (p>(int)estremi[0][0]) ){
     //cout<<"piano "<<p <<"; strip "<<calorimetro[p][0]<<"; en "<<calorimetro[p][1]<<endl;  
   //    if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){  
222      if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){      if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){
223        estremi[1][0]=p;//era p        estremi[1][0]=p;//
224        estremi[1][1]=calorimetro[p][1] *MIP;//energia in MeV        estremi[1][1]=calorimetro[p][1] *MIP;//energia in MeV
225      };      };
226      p = p-1;      p = p-1;
227    };    };
228    //    //
   
229        
230    /*integrale: energia totale rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/    /*integrale: energia totale rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/
231    for(Int_t pl=0; pl<(2*NPLA); pl++){    for(Int_t pl=0; pl<(2*NPLA); pl++){
# Line 199  void CaloBragg::Process(Int_t ntr){ Line 234  void CaloBragg::Process(Int_t ntr){
234      // se non e'il 1o dopo l'Y (tutti i pari) c'e' il W                // se non e'il 1o dopo l'Y (tutti i pari) c'e' il W          
235      if(pl%2!=0){                                                              //equival W in Si      if(pl%2!=0){                                                              //equival W in Si
236        Integrale+= 0.5*((calorimetro[pl-1][1] * MIP)+(calorimetro[pl][1] * MIP))*(spessore[1]);        Integrale+= 0.5*((calorimetro[pl-1][1] * MIP)+(calorimetro[pl][1] * MIP))*(spessore[1]);
       //      cout<<" W "<<0.5*((calorimetro[pl-1][1] * MIP)+(calorimetro[pl][1] * MIP))*(spessore[1])<<endl;  
237      };      };
238    };    };
239        Integrale=24000;//Integrale*1000;
240        
241    /*z ed energia con media troncata*/    /*z ed energia con media troncata*/
242  //   cout<<"Media troncata"<<endl;    //  mediatroncata();  // out:  1)chi2, 2)z, 3)Etot, 4)Pskip
   mediatroncata();  // out:  1)chi2, 2)z, 3)Etot, 4)Pskip  
243        
244    /*z ed energia con loop*/    /*z ed energia con loop*/
245  //    cout<<"Zdaloop"<<endl;      Zdaloop(); // out:  1)chi2, 2)z, 3)Etot, 4)Pskip
   //  Zdaloop(); // out:  1)chi2, 2)z, 3)Etot, 4)Pskip  
   
   
   /*energia rilasciata da z migliore*/  
  //  Float_t dEpianimean[2*NPLA];  
   //  Int_t zet=(int)bestchi2mean[1];  
   //Enetrack(&zet, &bestchi2mean[2], &estremi[0][0], dEpianimean);//calcola rilascio energetico sui piani  
   
   // Float_t dEpianiloop[2*NPLA];  
   //zet=(int)bestchi2loop[1];  
   //Enetrack(&zet, &bestchi2loop[2], &estremi[0][0], dEpianiloop);//calcola rilascio energetico sui piani  
246    
247        
248    if ( debug ) this->Print();    if ( debug ) this->Print();
# Line 233  void CaloBragg::Draw(){ Line 255  void CaloBragg::Draw(){
255    
256    Process();    Process();
257    
258   Float_t dEpianimean[44];  //  Float_t dEpianimean[44];
259   Float_t dEpianiloop[44];   Float_t dEpianiloop[44];
260   Float_t Depth[44];   Float_t Depth[44];
261   Int_t tz=(Int_t)qtz;  //  Int_t tz=(Int_t)qtz;
262   Int_t tz1=(Int_t)lpz;   Int_t tz1=(Int_t)lpz;
263   Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata  //  Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata
264   Enetrack(&tz1, &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop   Enetrack(&tz1, &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop
265    
266   Float_t sp= spessore[0]*spessore[1];   Float_t sp= spessore[0]*spessore[1];
# Line 247  void CaloBragg::Draw(){ Line 269  void CaloBragg::Draw(){
269    gStyle->SetLabelSize(0.04);    gStyle->SetLabelSize(0.04);
270    gStyle->SetNdivisions(510,"XY");    gStyle->SetNdivisions(510,"XY");
271    //    //
   
272          TString hid = Form("cCaloBragg");                TString hid = Form("cCaloBragg");      
273          TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));          TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
   
274          if ( tc ){          if ( tc ){
275  //       tc->Clear();  //       tc->Clear();
276          } else {          } else {
277           tc = new TCanvas(hid,hid);           tc = new TCanvas(hid,hid);
278           tc->Divide(1,2);  //       tc->Divide(1,2);
279          };          };
280          //          //
281          TString thid = Form("hCaloBragg");        //      TString thid = Form("hCaloBragg");      
282          TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));  //         TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
283          if ( th ) th->Delete();  //      if ( th ) th->Delete();
284  //       th->Clear();  //       th->Clear();
285  //       th->Reset();  //       th->Reset();
286  //      } else {  //      } else {
287          th = new TH2F(thid,thid,300,-0.5,300.,1000,0.,150.);  //      th = new TH2F(thid,thid,300,-0.5,300.,1000,0.,150.);
288          th->SetMarkerStyle(20);  //      th->SetMarkerStyle(20);
289  //      };  //      };
290          //          //
291          TString thid2 = Form("hCaloBragg2");              TString thid2 = Form("hCaloBragg2");    
# Line 285  void CaloBragg::Draw(){ Line 305  void CaloBragg::Draw(){
305    
306          tc->cd(1);          tc->cd(1);
307          //          //
308          for(Int_t i=0;i<=estremi[1][0];i++)th->Fill(Depth[i],dEpianimean[i]);  //      for(Int_t i=0;i<=estremi[1][0];i++)th->Fill(Depth[i],dEpianimean[i]);
309          for(Int_t i=0;i<=estremi[1][0];i++)th2->Fill(Depth[i],calorimetro[i][1]*MIP);          for(Int_t i=0;i<=estremi[1][0];i++)th2->Fill(Depth[i],calorimetro[i][1]*MIP);
310          th->Draw();  //      th->Draw();
311          th2->Draw("same");          th2->Draw("same");
312    
313          tc->cd(2);          tc->cd(2);
# Line 393  void CaloBragg::conversione(){ Line 413  void CaloBragg::conversione(){
413    Float_t ytgy = 0;    Float_t ytgy = 0;
414    Float_t a = 0.;    Float_t a = 0.;
415    
   
416    /*silicio*/    /*silicio*/
417    ytgx = ySi * L2->GetCaloLevel2()->tanx[0];    ytgx = ySi * L2->GetCaloLevel2()->tanx[0];
418    ytgy = ySi * L2->GetCaloLevel2()->tany[0];    ytgy = ySi * L2->GetCaloLevel2()->tany[0];
# Line 403  void CaloBragg::conversione(){ Line 422  void CaloBragg::conversione(){
422        
423    spessore[0] = SiCross/10. * rhoSi; //spessore silicio in g/cm2    spessore[0] = SiCross/10. * rhoSi; //spessore silicio in g/cm2
424        
     
425    /*tungsteno*/    /*tungsteno*/
426    ytgx = yW * L2->GetCaloLevel2()->tanx[0];    ytgx = yW * L2->GetCaloLevel2()->tanx[0];
427    ytgy = yW * L2->GetCaloLevel2()->tany[0];    ytgy = yW * L2->GetCaloLevel2()->tany[0];
# Line 417  void CaloBragg::conversione(){ Line 435  void CaloBragg::conversione(){
435    //       (g/cm2W)/(g/cm2Si)    //       (g/cm2W)/(g/cm2Si)
436    spessore[1] =  a;    spessore[1] =  a;
437    
     
438    //riscala mip allo spessore attraversato    //riscala mip allo spessore attraversato
439    spessore[2] = MIP*(SiCross/ySi);      spessore[2] = MIP*(SiCross/ySi);  
440    
# Line 457  void CaloBragg::BetheBloch(Float_t *x, F Line 474  void CaloBragg::BetheBloch(Float_t *x, F
474      (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(ISi,3.);
475        
476    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);      
   //   spessorecm  x ??/massSi x Zsi  
   // Energia = (*x)* rhoSi * 0.307/28.09 * 14. *SQ(*z)/SQ(*Bet)*(0.5*log(lg) - SQ(*Bet) - C/14.);  
     
477    
478    /*se ho gia' lo spessore in g/cm2 non se mejo???*/     Energia = (*x) * 0.307/28.09 * 14. *SQ(*z)/SQ(*Bet)*(0.5*log(lg) - SQ(*Bet) - C/14.);
   Energia = (*x) * 0.307/28.09 * 14. *SQ(*z)/SQ(*Bet)*(0.5*log(lg) - SQ(*Bet) - C/14.);  
479    
480    *out =Energia;//out    *out =Energia;//out
481    
# Line 482  void CaloBragg::ELOSS(Float_t *dx, Int_t Line 495  void CaloBragg::ELOSS(Float_t *dx, Int_t
495    
496    
497    Float_t Q=0.;    Float_t Q=0.;
498    Float_t v=0.;    Float_t v=0.;                                                        
                                                         //  Float_t Mass=0.;  
499    Float_t gam=0.;    Float_t gam=0.;
500    Float_t Bet=0.;    Float_t Bet=0.;
501    Float_t dEP=0.;    Float_t dEP=0.;
# Line 501  void CaloBragg::ELOSS(Float_t *dx, Int_t Line 513  void CaloBragg::ELOSS(Float_t *dx, Int_t
513    
514    //perdita energia per un protone    //perdita energia per un protone
515    Float_t protone =1.;    Float_t protone =1.;
516    Float_t Mass=MassP;    Float_t Mass=(elem[*Z-1]*MassP);
517    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);
518        
519    *out= (SQ(Q)*(dEP));//*dx;    *out= (SQ(Q)*(dEP));//*dx;
# Line 532  void CaloBragg::Enetrack(Int_t* Z, Float Line 544  void CaloBragg::Enetrack(Int_t* Z, Float
544        
545    Float_t Massa = (elem[(*Z)-1] * MassP);    Float_t Massa = (elem[(*Z)-1] * MassP);
546    
547      for( Int_t ipla=((int)(*primo)); ipla<= ((int)(*ultimo)); ipla++){
   //loop piani      (dal primo in cui entra)  
   //  for( Int_t ipla=((int)(*primo)); ipla< (2*NPLA); ipla++){  
  for( Int_t ipla=((int)(*primo)); ipla<= ((int)(*ultimo)); ipla++){  
548      dE=0.;      dE=0.;
   
549      //spessore silicio corretto x inclinazione, z, energia, out:rilascio      //spessore silicio corretto x inclinazione, z, energia, out:rilascio
550      ELOSS(&spessore[0], Z, &Ezero, &dE);//spessore in g/cm2!!      ELOSS(&spessore[0], Z, &Ezero, &dE);//spessore in g/cm2!!
551      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
# Line 548  void CaloBragg::Enetrack(Int_t* Z, Float Line 556  void CaloBragg::Enetrack(Int_t* Z, Float
556        out[ipla] = dE; //MeV        out[ipla] = dE; //MeV
557        Ezero = Ezero - dE;//energia residua        Ezero = Ezero - dE;//energia residua
558      };      };
   
559      //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
560      if(ipla%2 == 0){      if(ipla%2 == 0){
561        /*tungsteno*/        /*tungsteno*/
562        dE=0.;        dE=0.;
563        Float_t sp= spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato  in g/cm2        Float_t sp= spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato  in g/cm2
564        ELOSS(&sp, Z, &Ezero, &dE);        ELOSS(&sp, Z, &Ezero, &dE);
       //      cout<<"perdita per piano di W ="<<dE<<endl;  
565        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
566          return;          return;
567        }else{        }else{
568          Ezero = Ezero -dE;//energia residua          Ezero = Ezero -dE;//energia residua
         //      cout<<"w calc "<<dE<<endl;  
569        };        };
570      };          };    
   };//fine loop piani  
571    
572  //   for(Int_t i=0;i< 44; i++)cout<<"deposito energetico (teorico) per il piano i"<<i<<", "<<out[i]<<endl;    };//fine loop piani
     
 };//end Enetrack  
573    
574    
575    };//end Enetrack
576    
577    
578    
# Line 590  void CaloBragg::chiquadro(Float_t dE[], Line 593  void CaloBragg::chiquadro(Float_t dE[],
593    Float_t w,wi;    Float_t w,wi;
594    
595    for(Int_t ipla=0; ipla<2*NPLA; ipla++){    for(Int_t ipla=0; ipla<2*NPLA; ipla++){
     
596      //tutti i piani attraversati dalla traiettoria      //tutti i piani attraversati dalla traiettoria
597      if(calorimetro[ipla][0] != -1.){      if(calorimetro[ipla][0] != -1.){ //
598        w=0.; //normalizzazione;        w=0.; //normalizzazione;
599        wi=1.;//peso        wi=1.;//peso
600    
# Line 604  void CaloBragg::chiquadro(Float_t dE[], Line 606  void CaloBragg::chiquadro(Float_t dE[],
606        if((ipla>estremi[1][0])) wi=0.;        if((ipla>estremi[1][0])) wi=0.;
607    
608        //normalizzazione        //normalizzazione
609        if (calorimetro[ipla][1] != 0.)  w=1./(calorimetro[ipla][1]* MIP);                if (calorimetro[ipla][1] != 0.)  w=1./(calorimetro[ipla][1]* MIP);    //  
610                
611        //tolgo piani con rilasci inferiori al 30% del precedente        //tolgo piani con rilasci inferiori al 30% del precedente
612        if(calorimetro[ipla][1] <= (0.7*PianoPrecedente)){        if(calorimetro[ipla][1] < (0.7*PianoPrecedente)){ // cosi' i piani senza rilascio non vengono considerati nel calcolo del chi2
613          wi=0.;          wi=0.;
614          //se sono piani intermedi (non si è fermta) li considero non buoni          //se sono piani intermedi (non si e' fermta) li considero non buoni
615          if( (ipla <= estremi[1][0]) && (calorimetro[ipla][1] !=0.)){          if( (ipla <= estremi[1][0]) && (calorimetro[ipla][1] !=0.)){//
616            badplane+=1.;            badplane+=1.;
617            badplanetot+=1.;            badplanetot+=1.;
618          };          };
# Line 622  void CaloBragg::chiquadro(Float_t dE[], Line 624  void CaloBragg::chiquadro(Float_t dE[],
624        Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));        Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));
625    
626        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));
627          if(debug){
628            printf("dedx  calcolata  %f e reale  %f  \n",dE[ipla],(calorimetro[ipla][1] * MIP));
629          }
630        //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
631        if(wi != 0.){        if(wi != 0.){//
632          PianoPrecedente= calorimetro[ipla][1];//tengo piano precedente          PianoPrecedente= calorimetro[ipla][1];//tengo piano precedente
633          badplane = 0.;//azzero contatore piani scartati consecutivi          badplane = 0.;//azzero contatore piani scartati consecutivi
634        };        };
635      };      };
636              
637        //da Emi
638      if(badplane > 2){            if(badplane > 2){      
639        out[1] =79.;        out[1] =79.;
640        break;        break;
641      };        };
642                
643    };//fine loop piani    };//fine loop piani
644    //chi2,frammentato,pskip      //chi2,frammentato,pskip  
645    out[0]=sum;    out[0]=sum;
646    out[2]=badplanetot;    out[2]=badplanetot;
   
 //  if(out[1] ==79.)cout<<"frammentato !!!!!"<<endl;  
647        
648  };//end chiquadro  };//end chiquadro
649    
650    
651    
 //loopze(Integrale,&zero,&zmin, &zmax, spess1 , estr3, calo1, bestchi2);  
652  void CaloBragg::loopze( Float_t step, Float_t E0,Float_t Zstart, Float_t Zlimite){  void CaloBragg::loopze( Float_t step, Float_t E0,Float_t Zstart, Float_t Zlimite){
653  //  //
654    //loop su z ed energie per trovare miglior z (ed energia)    //loop su z ed energie per trovare miglior z (ed energia)
# Line 659  void CaloBragg::loopze( Float_t step, Fl Line 661  void CaloBragg::loopze( Float_t step, Fl
661    //    //
662    
663    
   
664    Float_t dEplan[2*NPLA];//energia rilasciata calcolata    Float_t dEplan[2*NPLA];//energia rilasciata calcolata
665    memset(dEplan,0,2*NPLA*sizeof(Float_t));    memset(dEplan,0,2*NPLA*sizeof(Float_t));
666    
# Line 667  void CaloBragg::loopze( Float_t step, Fl Line 668  void CaloBragg::loopze( Float_t step, Fl
668    
669    Float_t Massa = 0.;    Float_t Massa = 0.;
670    
671    Float_t Stepint =(step)/100.;//passo per il calcolo di energia //era 1000    Float_t Stepint =(step)/1000.;//passo per il calcolo di energia
672    
673    Float_t energia =0.;//energia del loop    Float_t energia =0.;//energia del loop
674    
# Line 691  void CaloBragg::loopze( Float_t step, Fl Line 692  void CaloBragg::loopze( Float_t step, Fl
692      Massa = elem[inucl-1]*MassP;      Massa = elem[inucl-1]*MassP;
693            
694      //loop energia      //loop energia
695      for(Int_t iene= 0; iene<1000; iene++){      for(Int_t iene= 0; iene<1000; iene++){// da non cambiare in base a Stepint altrimenti cambia la matrice bestchi2!!!cosi' non raggiungo mai integrale!!!!! mettere <=??
696            
697        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..
698    
# Line 700  void CaloBragg::loopze( Float_t step, Fl Line 701  void CaloBragg::loopze( Float_t step, Fl
701        //calcolo chi2        //calcolo chi2
702        chiquadro(dEplan,chi2);        chiquadro(dEplan,chi2);
703                
704        if( (chi2[1] != 79) ){//salto quelli che frammentano        if( (chi2[1] != 79.) ){//salto quelli che frammentano
705          matrixchi2[inucl][iene][0]=chi2[0];//valore chi2 per questo z a questa energia          matrixchi2[inucl][iene][0]=chi2[0];//valore chi2 per questo z a questa energia
706          matrixchi2[inucl][iene][1]=energia;//energia per questo chi2          matrixchi2[inucl][iene][1]=energia;//energia per questo chi2
707          matrixchi2[inucl][iene][2]=chi2[2];//piani saltati nel chi2          matrixchi2[inucl][iene][2]=chi2[2];//piani saltati nel chi2
# Line 715  void CaloBragg::loopze( Float_t step, Fl Line 716  void CaloBragg::loopze( Float_t step, Fl
716    
717   };//fine loop z   };//fine loop z
718    
719    //  for (Int_t nu=0; nu<colmax; nu++){  
720    //    for (Int_t en=0; en<rowmax; en++){    //Emi
721    for (Int_t nu=(int)(Zstart); nu<max; nu++){    for (Int_t nu=(int)(Zstart); nu<max; nu++){
722      for (Int_t en=0; en<1000; en++){      for (Int_t en=0; en<1000; en++){
723        if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] !=0.)){        if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] >0.)){
724          bestchi2[0]= matrixchi2[nu][en][0];// chi2          bestchi2[0]= matrixchi2[nu][en][0];// chi2
725          bestchi2[1]= (Float_t)nu; // z          bestchi2[1]= (Float_t)nu; // z
726          bestchi2[2]= matrixchi2[nu][en][1];//energia;          bestchi2[2]= matrixchi2[nu][en][1];//energia;
727          bestchi2[3]=    matrixchi2[nu][en][2];// totale piani saltati          bestchi2[3]= matrixchi2[nu][en][2];// totale piani saltati
728        }        }
729      }      }
730    }    }
   //==========================//  
   //  Int_t testz = (Int_t)bestchi2[1];  
   //  Enetrack(&testz, &bestchi2[2], &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani  
   //  for(Int_t i=0;i<=estremi[1][0];i++){  
    // cout<<"dEplan "<<dEplan[i]*10<<endl;  
    // cout<<"calorimetro "<<calorimetro[i][1]<<endl;  
   //  }  
   //==========================//  
 };//endloopze  
   
   
   
   
   
 void CaloBragg::mediatroncata(){  
   //calcolo Z con media troncata e utilizzo questo Z per trovare l'energia migliore  
   //in: ordplane[44]   => array con energia dei piani  
   //    spess[3]    => conversioni spessore di silicio, w,  mip  
   //    estr[2][2]  => primo[0][0] e ultimo[1][0] piano attraversati ed energie[][1]  
   //    calo[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano  
   //    integrale      => energia totale nel calorimetro considerando il W  
   //  
   // out[4] chi2,z,Etot,Pskip  
   
   Float_t ordplane[44];//mi serve per la media troncata  
   memset(ordplane,0,44*sizeof(Float_t));  
731    
   for(Int_t ipla=0; ipla< 2*NPLA; ipla++)  ordplane[ipla]=calorimetro[ipla][1]; //energia del piano  
   
     
   //ordino tutte le energie dei piani in ordine crescente  
   //Int_t ii=0;  
732    
733      Long64_t work[200];  };//endloopze
     Int_t ind = 0;  
     //    Int_t l = 0;  
     Int_t RN = 0;  
     Float_t sum4 = 0.;  
     Float_t qm = 0.;  
     //  
     //Float_t qmt = ethr*0.8; // *0.9  
     //  
     //Int_t uplim = TMath::Max(3,N);  
     //  
     while ( RN < 4 && ind < 44 ){  
       qm = TMath::KOrdStat(44,ordplane,ind,work);  
       if (qm >= 0.7 ){    
         if ( RN < 4 ){  
           sum4 += qm;  
           RN++;  
         };  
         //      l++;  
         //      if ( debug ) printf(" value no %i qm %f sum4 %f \n",l,qm,sum4);  
       };  
       ind++;  
     };  
     //  
     sum4 /= (Float_t)RN;  
     Float_t Zmean = (sqrt((sum4*MIP)/(4.*spessore[2])));  
     if(Zmean ==0.) Zmean=1.;  
     if ( Zmean < 1. ) Zmean = 1.;  
     //  
     //    Zmean = 1.;  
734    
 //     cout<<"sum4="<<sum4<<endl;  
 //     cout<<"Zmean="<<Zmean<<endl;  
735    
     
   //calcolo energia migliore per Z trovato con media troncata  
     
   Float_t zmin=Zmean;  
     
   bestchi2[0]=10000.;  
   bestchi2[1]=0.;  
   bestchi2[2]=0.;  
   bestchi2[3]=0.;  
   Float_t zero=0.;  
736    
  // cout<<"inizio media troncata"<<endl;  
  // cout<<"input: step "<<Integrale<<", zero "<<zero<<", zmin "<<zmin<<", zmax"<<zmin<<endl;  
737    
738    
739      //    step  energia   zstart  zstop    // void CaloBragg::mediatroncata(){
740    loopze(Integrale,zero,zmin,zmin);  //   //calcolo Z con media troncata e utilizzo questo Z per trovare l'energia migliore
741      //   //in: ordplane[44]   => array con energia dei piani
742   // cout<<"fine media troncata"<<endl;  //   //    spess[3]    => conversioni spessore di silicio, w,  mip
743   // cout<<"output: chi2 "<<bestchi2[0]<<", z"<<bestchi2[1]<<", Etot "<<bestchi2[2]<<", Pskip "<<bestchi2[3]<<endl;  //   //    estr[2][2]  => primo[0][0] e ultimo[1][0] piano attraversati ed energie[][1]
744    //   //    calo[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano
745    //   //    integrale      => energia totale nel calorimetro considerando il W
746    //   //
747    //   // out[4] chi2,z,Etot,Pskip
748    
749    //   Float_t ordplane[44];//mi serve per la media troncata
750    //   memset(ordplane,0,44*sizeof(Float_t));
751    
752    //   for(Int_t ipla=0; ipla< 2*NPLA; ipla++)  ordplane[ipla]=calorimetro[ipla][1]; //energia del piano
753    
754      
755    //   //ordino tutte le energie dei piani in ordine crescente
756    
757    //     Long64_t work[200];
758    //     Int_t ind = 0;
759    //     //Int_t l = 0;
760    //     Int_t RN = 0;
761    //     Float_t sum4 = 0.;
762    //     Float_t qm = 0.;
763    //     //
764    //     //Float_t qmt = ethr*0.8; // *0.9
765    //     //
766    //     //Int_t uplim = TMath::Max(3,N);
767    //     //
768    //     while ( RN < 4 && ind < 44 ){
769    //       qm = TMath::KOrdStat(44,ordplane,ind,work);
770    //       if (qm >= 0.7 ){      
771    //      if ( RN < 4 ){
772    //        sum4 += qm;
773    //        RN++;
774    //      };
775    // //   l++;
776    // //   if ( debug ) printf(" value no %i qm %f sum4 %f \n",l,qm,sum4);
777    //       };
778    //       ind++;
779    //     };
780    //     //
781    //     sum4 /= (Float_t)RN;
782    //     Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2])));//ma non e'/1??
783    //     if(Zmean ==0.) Zmean=1.;
784    //     if ( Zmean < 1. ) Zmean = 1.;
785    
786        
787    // //     Zmean =round(Zmean);
788    // //     if(Zmean <1.) Zmean=1.;
789    
790    // //     if(Zmean >0.)Zmean =round(Zmean);
791    
792    qtchi2=bestchi2[0];  //     //======== per i nuclei=======
793    qtz=bestchi2[1];  //     if (Zmean >=2.){
794    qtetot=bestchi2[2];  //     ind = 0;
795    qtpskip=bestchi2[3];  //     RN = 0;
796  };//end mediatroncata  //     sum4 = 0.;
797    //  cout<<"z media troncata ok"<<endl;  //     qm = 0.;
798    //     while ( RN < 4 && ind < 44 ){
799    //       qm = TMath::KOrdStat(44,ordplane,ind,work);
800    //       if (qm >= (Zmean*Zmean)-Zmean*Zmean*0.2 ){    
801    //      if ( RN < 4 ){
802    //        sum4 += qm;
803    //        RN++;
804    //      };
805    //       };
806    //       ind++;
807    //     };
808    //     //
809    //     sum4 /= (Float_t)RN;
810    //     Zmean = (sqrt((sum4*MIP)/(4.*spessore[2])));//ma non e' /1??
811    //     }
812    
813      
814    //   //calcolo energia migliore per Z trovato con media troncata
815    //   //  Float_t zmin=Zmean;
816    //     Float_t zmin=round(Zmean);
817      
818    //   bestchi2[0]=10000.;
819    //   bestchi2[1]=0.;
820    //   bestchi2[2]=0.;
821    //   bestchi2[3]=0.;
822    //   Float_t zero=0.;
823    
824    //   //    step   energia zstart zstop  
825    //   loopze(Integrale,zero,zmin,zmin);
826      
827    
828    //   qtchi2=bestchi2[0];
829    //   qtz=bestchi2[1];
830    //   qtetot=bestchi2[2];
831    //   qtpskip=bestchi2[3];
832    // };//end mediatroncata
833    
834    
835    
# Line 838  void CaloBragg::Zdaloop(){ Line 845  void CaloBragg::Zdaloop(){
845    
846    
847    /*z se particella fosse al minimo*/  //energia1piano/mip corretta    /*z se particella fosse al minimo*/  //energia1piano/mip corretta
848    Float_t zmax = (sqrt(estremi[0][1]/spessore[2]));    Float_t zmax = round(sqrt(estremi[0][1]/spessore[2]));
849    if(zmax<31)zmax=zmax+1;    if(zmax<31)zmax=zmax+1;
850        
851    /*calcolo Z ed E con loop sui vari elementi ed energie*/    /*calcolo Z ed E con loop sui vari elementi ed energie*/
# Line 851  void CaloBragg::Zdaloop(){ Line 858  void CaloBragg::Zdaloop(){
858    bestchi2[2]=0.;    bestchi2[2]=0.;
859    bestchi2[3]=0.;    bestchi2[3]=0.;
860    Float_t zero=0.;    Float_t zero=0.;
861    //primo loop      //------------primo loop   ----------------------
862    //     energia   ezero, zstart  zstop Si attrav  1 piano  energie piani  out    //     energia   ezero, zstart  zstop
   
  // cout<<"inizio primo loop"<<endl;  
  // cout<<"input: step "<<Integrale<<", zero "<<zero<<", zmin "<<zmin<<", zmax"<<zmax<<endl;  
   
863    loopze(Integrale,zero,zmin,zmax);    loopze(Integrale,zero,zmin,zmax);
864    
 //  cout<<"fine primo loop"<<endl;  
 //  cout<<"output: chi2 "<<bestchi2[0]<<", z"<<bestchi2[1]<<", Etot "<<bestchi2[2]<<", Pskip "<<bestchi2[3]<<endl;  
865        
866    //secondo loop    //------------secondo loop  ----------------------
867    for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i];    for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i];
868    bestchi2[0] = 10000.;    bestchi2[0] = 10000.;
869    bestchi2[1] = 0.;    bestchi2[1] = 0.;
870    bestchi2[2] = 0.;    bestchi2[2] = 0.;
871    bestchi2[3] = 0.;//riazzero    bestchi2[3] = 0.;//riazzero
872    
873    Float_t step = bestchitemp[2]/100.;//Integrale/10.;//era 1000    Float_t step = bestchitemp[2];//
874    //zero=bestchitemp[2]-step/2.;//energia da 1 giro - 1step    zero=0;  // qualsiasi altro valore peggiora le cose
875    //Float_t step = bestchitemp[2]/1000.;    zmin=zmax=bestchitemp[1];
   //zero=bestchitemp[2]-step;  
   zero=0;  
   zmin=bestchitemp[1]-2;// zda primo giro  
   if(zmin<1)zmin=1;  
   zmax=bestchitemp[1]+1;//  
   
 //  cout<<"inizio secondo loop"<<endl;  
 //  cout<<"input: step "<<step<<", zero "<<zero<<", zmin "<<zmin<<", zmax"<<zmax<<endl;  
   
876    loopze(step,zero,zmin,zmax); //    loopze(step,zero,zmin,zmax); //
877    
 //  cout<<"fine secondo loop"<<endl;  
 //  cout<<"output: chi2 "<<bestchi2[0]<<", z"<<bestchi2[1]<<", Etot "<<bestchi2[2]<<", Pskip "<<bestchi2[3]<<endl;  
   
   
 //   cout<<"z loop ok"<<endl;  
878        
879    //chi2,z,Etot,Pskip    //chi2,z,Etot,Pskip
880    lpchi2=bestchi2[0];    lpchi2=bestchi2[0];
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.23