/[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.1 by pamelats, Fri Jun 13 08:11:04 2008 UTC revision 1.2 by pamelats, Wed Jun 18 15:24:12 2008 UTC
# Line 109  void CaloBragg::Process(Int_t ntr){ Line 109  void CaloBragg::Process(Int_t ntr){
109    
110    //    //
111    //    //
   // Always calculate stdedx1  
112    //    //
113    Int_t view = 0;    Int_t view = 0;
114    Int_t plane = 0;    Int_t plane = 0;
# Line 134  void CaloBragg::Process(Int_t ntr){ Line 133  void CaloBragg::Process(Int_t ntr){
133        ptrack = L2->GetTrack(ntr);        ptrack = L2->GetTrack(ntr);
134        if ( ptrack ) track = ptrack->GetCaloTrack();        if ( ptrack ) track = ptrack->GetCaloTrack();
135      } else {      } else {
136        track = L2->GetCaloStoredTrack(ntr);        track = L2->GetCaloStoredTrack(ntr);  //al momento e' vera solo questa riga
137      };      };
138      //      //
139      if ( !track && ntr >= 0 ){      if ( !track && ntr >= 0 ){
# Line 155  void CaloBragg::Process(Int_t ntr){ Line 154  void CaloBragg::Process(Int_t ntr){
154    //    //
155    for(Int_t p=0; p<22; p++){    for(Int_t p=0; p<22; p++){
156      for(Int_t v=0; v<2; v++){      for(Int_t v=0; v<2; v++){
       cout<< L2->GetCaloLevel2()->cibar[p][v];      
157        /*per usare traccia non del calo camboare cibar*/        /*per usare traccia non del calo camboare cibar*/
158        calorimetro[(2*p)+1-v][0] = L2->GetCaloLevel2()->cibar[p][v];//strip attraversata        calorimetro[(2*p)+1-v][0] = L2->GetCaloLevel2()->cibar[p][v];//strip attraversata
159        calorimetro[(2*p)+1-v][1] = (epiano[p][v]); //energia del piano //(epiano[p][v])/0.89        calorimetro[(2*p)+1-v][1] = (epiano[p][v]); //energia del piano //(epiano[p][v])/0.89
160      };      };
161    };    };
162    
163      /*per ogni evento calcolo la conversione mip e w attraversato in equivalente Si*/
   /*per ogni evento clacolo la conversione mip e w attraversato in equivalente Si*/  
   
164    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
165      //cout<<"spessore= "<<spessore[0]<<" gcm2si; "<<spessore[1]<<"W in Si; "<<spessore[2]<<"mip "<<endl;
166    
   cout<<"spessore"<<spessore[0]<<" gcm2si; "<<spessore[1]<<"W in Si; "<<spessore[2]<<"mip "<<endl;  
167    
   /*trova primo e ultimo piano attraversati*/  
168    
169      /*trova primo e ultimo piano attraversati*/
170    Int_t p = 0;//contatore piani    Int_t p = 0;//contatore piani
171    //primo parte da 0 e va in giù    //per il primo parte da 0 e va in giù
172    while( estremi[0][1] == 0  &&  p<(2*NPLA) ){    while( estremi[0][1] == 0  &&  p<(2*NPLA) ){
173      cout<<"piano "<<p <<"; strip "<<calorimetro[p][0]<<"; en "<<calorimetro[p][1]<<endl;      //cout<<"piano "<<p <<"; strip "<<calorimetro[p][0]<<"; en "<<calorimetro[p][1]<<endl;
174      //    if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){      //    if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){
175      if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){ //0.7 soglia minima      if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){ //0.7 soglia minima
176        estremi[0][0]=p;        estremi[0][0]=p;
# Line 185  void CaloBragg::Process(Int_t ntr){ Line 181  void CaloBragg::Process(Int_t ntr){
181   //ultimo parte da 44 e sale   //ultimo parte da 44 e sale
182    p=43;    p=43;
183    while( (estremi[1][1] == 0)  &&  (p>(int)estremi[0][0]) ){    while( (estremi[1][1] == 0)  &&  (p>(int)estremi[0][0]) ){
184    cout<<"piano "<<p <<"; strip "<<calorimetro[p][0]<<"; en "<<calorimetro[p][1]<<endl;      //cout<<"piano "<<p <<"; strip "<<calorimetro[p][0]<<"; en "<<calorimetro[p][1]<<endl;
185    //    if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){    //    if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){
186    if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){      if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){
187        estremi[1][0]=p;        estremi[1][0]=p;//era p
188        estremi[1][1]=calorimetro[p][1] *MIP;//energia in MeV        estremi[1][1]=calorimetro[p][1] *MIP;//energia in MeV
189      };      };
190      p = p-1;      p = p-1;
191    };    };
192    //    //
   cout<<"estremi: in "<<estremi[0][0]<<"piano, "<<estremi[0][1]<<"energia; out "<<estremi[1][0]<<"piano, "<<estremi[1][1]<<"energia."<<endl;  
   
193    
194        
195    /*integrale: energia rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/    /*integrale: energia totale rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/
196      for(Int_t pl=0; pl<(2*NPLA); pl++){
   for(Int_t pl=0; pl<(2*NPLA); pl++){  
   
197      //calcolo intergale in unita di spessori di silicio      //calcolo intergale in unita di spessori di silicio
198      Integrale += calorimetro[pl][1] * MIP;//piano di silicio      Integrale += calorimetro[pl][1] * MIP;//piano di silicio
     // cout<<Integrale<<" piani "<<pl  
     cout<<"nel calo ho strip "<<calorimetro[pl][1]*MIP<<endl;  
199      // 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          
200      if(pl%2!=0){                                                             //equival W in Si      if(pl%2!=0){                                                              //equival W in Si
201        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]);
202        // cout<<Integrale<<" W "<<pl<<endl;        //      cout<<" W "<<0.5*((calorimetro[pl-1][1] * MIP)+(calorimetro[pl][1] * MIP))*(spessore[1])<<endl;
       cout<<" W "<<0.5*((calorimetro[pl-1][1] * MIP)+(calorimetro[pl][1] * MIP))*(spessore[1])<<endl;  
203      };      };
204    };    };
   
   //  cout<<"integrale ok"<<endl;  
205        
206        
   /*z se particella fosse al minimo*/  //energia1piano/mip corretta  
   //Float_t zmax = (sqrt(estremi[0][1]/spessore[2]));  
     
   
     
207    /*z ed energia con media troncata*/    /*z ed energia con media troncata*/
208    //  Float_t  bestchi2mean[4] = {0.,0.,0.,0.};//chi2,z,Etot,Pskip  //   cout<<"Media troncata"<<endl;
209    mediatroncata();      mediatroncata();  // out:  1)chi2, 2)z, 3)Etot, 4)Pskip
210        
211    /*z ed energia con loop*/    /*z ed energia con loop*/
212    //  Float_t  bestchi2loop[4] = {0.,0.,0.,0.};//chi2,z,Etot,Pskip//chi2,z,Etot,Pskip  //    cout<<"Zdaloop"<<endl;
213    Zdaloop();    Zdaloop(); // out:  1)chi2, 2)z, 3)Etot, 4)Pskip
214      
215    
216    /*energia rilasciata da z migliore*/    /*energia rilasciata da z migliore*/
217   //  Float_t dEpianimean[2*NPLA];   //  Float_t dEpianimean[2*NPLA];
# Line 242  void CaloBragg::Process(Int_t ntr){ Line 224  void CaloBragg::Process(Int_t ntr){
224    
225        
226    if ( debug ) this->Print();    if ( debug ) this->Print();
227    if ( debug ) printf(" esci \n");    if ( debug ) printf(" fine evento \n");
228    //    //
229  };  };
230    
# Line 256  void CaloBragg::Draw(){ Line 238  void CaloBragg::Draw(){
238   Float_t Depth[44];   Float_t Depth[44];
239   Int_t tz=(Int_t)qtz;   Int_t tz=(Int_t)qtz;
240   Int_t tz1=(Int_t)lpz;   Int_t tz1=(Int_t)lpz;
241   Enetrack(&tz, &qtetot, &estremi[0][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
242   Enetrack(&tz1, &lpetot, &estremi[0][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
243    
244   Float_t sp= spessore[0]*spessore[1];   Float_t sp= spessore[0]*spessore[1];
245   for(Int_t i=0;i<44;i++)Depth[i]=i*sp;   for(Int_t i=0;i<44;i++)Depth[i]=i*sp;
246    //    //
247    gStyle->SetLabelSize(0.04);    gStyle->SetLabelSize(0.04);
248    gStyle->SetNdivisions(510,"XY");    gStyle->SetNdivisions(510,"XY");
   //gStyle->SetLogy();  
249    //    //
250    
251          TString hid = Form("cCaloBragg");                TString hid = Form("cCaloBragg");      
# Line 304  void CaloBragg::Draw(){ Line 285  void CaloBragg::Draw(){
285    
286          tc->cd(1);          tc->cd(1);
287          //          //
288          for(Int_t i=0;i<44;i++)th->Fill(Depth[i],dEpianimean[i]);          for(Int_t i=0;i<=estremi[1][0];i++)th->Fill(Depth[i],dEpianimean[i]);
289          for(Int_t i=0;i<44;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);
290          th->Draw();          th->Draw();
291          th2->Draw("same");          th2->Draw("same");
292    
293          tc->cd(2);          tc->cd(2);
294          //          //
295          for(Int_t i=0;i<44;i++)th3->Fill(Depth[i],dEpianiloop[i]);          for(Int_t i=0;i<=estremi[1][0];i++)th3->Fill(Depth[i],dEpianiloop[i]);
296          th3->Draw();          th3->Draw();
297          th2->Draw("same");          th2->Draw("same");
298    
# Line 325  void CaloBragg::Draw(){ Line 306  void CaloBragg::Draw(){
306  };  };
307    
308    
 void CaloBragg::SWAP( Float_t *A, Float_t *B ){  
    Float_t  Tmp = *A;  
    *A = *B;  
    *B = Tmp;  
 };  
   
   
309    
310  void CaloBragg::LoadParam(){  void CaloBragg::LoadParam(){
311    
# Line 370  void CaloBragg::LoadParam(){ Line 344  void CaloBragg::LoadParam(){
344    elem[31] = 72.61; //Ge 32    elem[31] = 72.61; //Ge 32
345    
346    
   
347  //parametri calorimetro  //parametri calorimetro
348    NPLA = 22;    NPLA = 22;
349    NCHA = 96;    NCHA = 96;
# Line 472  void CaloBragg::BetheBloch(Float_t *x, F Line 445  void CaloBragg::BetheBloch(Float_t *x, F
445    Float_t lg =0.;    Float_t lg =0.;
446    Float_t Energia=0.;    Float_t Energia=0.;
447    Float_t C=0.;    Float_t C=0.;
   
448    
449    eta = (*gam)*(*Bet);    eta = (*gam)*(*Bet);
450    
# Line 481  void CaloBragg::BetheBloch(Float_t *x, F Line 453  void CaloBragg::BetheBloch(Float_t *x, F
453        
454    lg = 2.* Me * SQ(eta) * Wmax / SQ(ISi);    lg = 2.* Me * SQ(eta) * Wmax / SQ(ISi);
455                                                                                                 //  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;
456    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(ISi,2.) +
457      (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.);
458        
459    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);      
460    //   spessorecm  x ??/massSi x Zsi    //   spessorecm  x ??/massSi x Zsi
# Line 534  void CaloBragg::ELOSS(Float_t *dx, Int_t Line 506  void CaloBragg::ELOSS(Float_t *dx, Int_t
506        
507    *out= (SQ(Q)*(dEP));//*dx;    *out= (SQ(Q)*(dEP));//*dx;
508    
 };//end ELOSS  
   
509    
510    };//end ELOSS
511    
512    
513    
514    
515  void CaloBragg::Enetrack(Int_t* Z, Float_t* E0, Float_t* primo, Float_t out[]){  void CaloBragg::Enetrack(Int_t* Z, Float_t* E0, Float_t* primo,Float_t* ultimo, Float_t out[]){
516    
517    //calcola energia rilasciata sulla traccia (usa ELOSS)    //calcola energia rilasciata sulla traccia (usa ELOSS)
518    //  in : Z             =>carica    //  in : Z             =>carica
# Line 560  void CaloBragg::Enetrack(Int_t* Z, Float Line 531  void CaloBragg::Enetrack(Int_t* Z, Float
531    memset(out, 0, 2*NPLA*sizeof(Float_t));    memset(out, 0, 2*NPLA*sizeof(Float_t));
532        
533    Float_t Massa = (elem[(*Z)-1] * MassP);    Float_t Massa = (elem[(*Z)-1] * MassP);
     
534    
535    
536    //loop piani      (dal primo in cui entra)    //loop piani      (dal primo in cui entra)
537    for( Int_t ipla=((int)(*primo)); ipla< (2*NPLA); ipla++){    //  for( Int_t ipla=((int)(*primo)); ipla< (2*NPLA); ipla++){
538     for( Int_t ipla=((int)(*primo)); ipla<= ((int)(*ultimo)); ipla++){
539      dE=0.;      dE=0.;
540    
541      //spessore silicio corretto x inclinazione, z, energia, out:rilascio      //spessore silicio corretto x inclinazione, z, energia, out:rilascio
542      ELOSS(&spessore[0], Z, &Ezero, &dE);//spessore in g/cm2!!      ELOSS(&spessore[0], Z, &Ezero, &dE);//spessore in g/cm2!!
       
543      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
544        out[ipla] = Ezero - Massa; //MeV        out[ipla] = Ezero - Massa; //MeV
545        return;        return;
# Line 585  void CaloBragg::Enetrack(Int_t* Z, Float Line 555  void CaloBragg::Enetrack(Int_t* Z, Float
555        dE=0.;        dE=0.;
556        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
557        ELOSS(&sp, Z, &Ezero, &dE);        ELOSS(&sp, Z, &Ezero, &dE);
558                //      cout<<"perdita per piano di W ="<<dE<<endl;
559        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
560          return;          return;
561        }else{        }else{
# Line 594  void CaloBragg::Enetrack(Int_t* Z, Float Line 564  void CaloBragg::Enetrack(Int_t* Z, Float
564        };        };
565      };          };    
566    };//fine loop piani    };//fine loop piani
567    
568    //   for(Int_t i=0;i< 44; i++)cout<<"deposito energetico (teorico) per il piano i"<<i<<", "<<out[i]<<endl;
569        
570  };//end Enetrack  };//end Enetrack
571    
# Line 613  void CaloBragg::chiquadro(Float_t dE[], Line 585  void CaloBragg::chiquadro(Float_t dE[],
585    
586    Float_t sum = 0.;    Float_t sum = 0.;
587    Float_t PianoPrecedente=0.;    Float_t PianoPrecedente=0.;
   //  Float_t differenza =0.;  
588    Float_t badplane=0.;    Float_t badplane=0.;
589    Float_t badplanetot=0.;    Float_t badplanetot=0.;
590    Float_t w,wi;    Float_t w,wi;
# Line 629  void CaloBragg::chiquadro(Float_t dE[], Line 600  void CaloBragg::chiquadro(Float_t dE[],
600        if (ipla<estremi[0][0])  wi=0.;        if (ipla<estremi[0][0])  wi=0.;
601                
602        //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
603        //    if((ipla>estremi[1][0]) && (calorimetro[ipla][1] >0.) ) wi=0.;        //if((ipla>estremi[1][0]) && (calorimetro[ipla][1] >0.) ) wi=0.;
604          if((ipla>estremi[1][0])) wi=0.;
605    
606        //normalizzazione        //normalizzazione
607        if (calorimetro[ipla][1] != 0.)  w=1./(calorimetro[ipla][1]* MIP);                if (calorimetro[ipla][1] != 0.)  w=1./(calorimetro[ipla][1]* MIP);        
608                
609        //tolgo piani con rilasci inferiori al 30% del precedente        //tolgo piani con rilasci inferiori al 30% del precedente
610        if(calorimetro[ipla][1] <= (0.7*PianoPrecedente)){        if(calorimetro[ipla][1] <= (0.7*PianoPrecedente)){
           
611          wi=0.;          wi=0.;
612          //se sono piani intermedi (non si è fermta) li considero non buoni          //se sono piani intermedi (non si è fermta) li considero non buoni
613          if( (ipla < estremi[1][0]) && (calorimetro[ipla][1] !=0.)){          if( (ipla <= estremi[1][0]) && (calorimetro[ipla][1] !=0.)){
614            badplane+=1.;            badplane+=1.;
615            badplanetot+=1.;            badplanetot+=1.;
616          };          };
617        };        };
         
   
618        
619        //meno peso ai piani con rilasci maggiori di 1000 MIP        //meno peso ai piani con rilasci maggiori di 1000 MIP
620          if(calorimetro[ipla][1] > 1000) wi=0.5;          if(calorimetro[ipla][1] > 1000) wi=0.5;
621                            
         //do peso maggiore alle ultime 6 misure  
         //      if((ipla >= estr2[1][0]-12)) wi=0.;  
   
622        Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));        Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));
623    
624        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));
# Line 669  void CaloBragg::chiquadro(Float_t dE[], Line 635  void CaloBragg::chiquadro(Float_t dE[],
635    };//fine loop piani    };//fine loop piani
636    //chi2,frammentato,pskip      //chi2,frammentato,pskip  
637    out[0]=sum;    out[0]=sum;
   
638    out[2]=badplanetot;    out[2]=badplanetot;
639    
640      if(out[1] ==79.)cout<<"frammentato !!!!!"<<endl;
641        
642  };//end chiquadro  };//end chiquadro
643    
# Line 684  void CaloBragg::loopze( Float_t step, Fl Line 651  void CaloBragg::loopze( Float_t step, Fl
651    //     E0        => energia iniziale (intergale)            //     E0        => energia iniziale (intergale)        
652    //     Zstart    => minimo z da cui patire    //     Zstart    => minimo z da cui patire
653    //     Zlimite   => z a cui fermarsi (z al minimo di ionizz sul 1o piano)    //     Zlimite   => z a cui fermarsi (z al minimo di ionizz sul 1o piano)
   //     spessore  => array conversione spessore Si, mip, W  
   //     estr1   => array con primo[0][0] e ultimo[1][0] piano attraversati ed energie([x][1])  
   //     calo2[44][2]=> [0]strip attraversata [1]energia misurata su ognuno dei 44 piani  
654    //    //
655    //out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2    //out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2
656    //    //
657    
658    
659    
660    Float_t dEplan[2*NPLA];//energia rilasciata calcolata    Float_t dEplan[2*NPLA];//energia rilasciata calcolata
661    memset(dEplan,0,2*NPLA*sizeof(Float_t));    memset(dEplan,0,2*NPLA*sizeof(Float_t));
662    
# Line 698  void CaloBragg::loopze( Float_t step, Fl Line 664  void CaloBragg::loopze( Float_t step, Fl
664    
665    Float_t Massa = 0.;    Float_t Massa = 0.;
666    
667    Float_t Stepint =(step)/1000.;//passo per il calcolo di energia    Float_t Stepint =(step)/100.;//passo per il calcolo di energia //era 1000
668    
669    Float_t energia =0.;//energia del loop    Float_t energia =0.;//energia del loop
670    
# Line 722  void CaloBragg::loopze( Float_t step, Fl Line 688  void CaloBragg::loopze( Float_t step, Fl
688      Massa = elem[inucl-1]*MassP;      Massa = elem[inucl-1]*MassP;
689            
690      //loop energia      //loop energia
691      for(Int_t iene= 0; iene<3000; iene++){      for(Int_t iene= 0; iene<1000; iene++){
692            
693        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..
694        //cout<<"folse"<<estremi[0][0];  
695        Enetrack(&Z, &energia, &estremi[0][0], dEplan);//calcola rilascio energetico sui piani        Enetrack(&Z, &energia, &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani
696                
697        //calcolo chi2        //calcolo chi2
698        chiquadro(dEplan,chi2);        chiquadro(dEplan,chi2);
# Line 734  void CaloBragg::loopze( Float_t step, Fl Line 700  void CaloBragg::loopze( Float_t step, Fl
700        if( (chi2[1] != 79) ){//salto quelli che frammentano        if( (chi2[1] != 79) ){//salto quelli che frammentano
701          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
702          matrixchi2[inucl][iene][1]=energia;//energia per questo chi2          matrixchi2[inucl][iene][1]=energia;//energia per questo chi2
703          matrixchi2[inucl][iene][2]=chi2[2];//piani saltati ne chi2          matrixchi2[inucl][iene][2]=chi2[2];//piani saltati nel chi2
704        } else {          } else {  
705          matrixchi2[inucl][iene][0]=1000;//valore chi2 per questo z a questa energia          matrixchi2[inucl][iene][0]=1000;//valore chi2 per questo z a questa energia
706          matrixchi2[inucl][iene][1]=1000;//energia per questo chi2          matrixchi2[inucl][iene][1]=1000;//energia per questo chi2
707          matrixchi2[inucl][iene][2]=1000;//piani saltati ne chi2          matrixchi2[inucl][iene][2]=1000;//piani saltati nel chi2
708        }        }
709      }//fine loop energia      }//fine loop energia
710    
# Line 749  void CaloBragg::loopze( Float_t step, Fl Line 715  void CaloBragg::loopze( Float_t step, Fl
715      for (Int_t en=0; en<rowmax; en++){      for (Int_t en=0; en<rowmax; en++){
716        if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] !=0.)){        if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] !=0.)){
717          bestchi2[0]= matrixchi2[nu][en][0];// chi2          bestchi2[0]= matrixchi2[nu][en][0];// chi2
718          bestchi2[1]= (Float_t)nu;          bestchi2[1]= (Float_t)nu; // z
719          bestchi2[2]= matrixchi2[nu][en][1];//energia;          bestchi2[2]= matrixchi2[nu][en][1];//energia;
720          bestchi2[3]=    matrixchi2[nu][en][2];// totale piani saltati          bestchi2[3]=    matrixchi2[nu][en][2];// totale piani saltati
721        }        }
722      }      }
723    }    }
724      //==========================//
725      Int_t testz = (Int_t)bestchi2[1];
726    
727      Enetrack(&testz, &bestchi2[2], &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani
728      for(Int_t i=0;i<=estremi[1][0];i++){
729        cout<<"dEplan "<<dEplan[i]*10<<endl;
730        cout<<"calorimetro "<<calorimetro[i][1]<<endl;
731      }
732      //==========================//
733  };//endloopze  };//endloopze
734    
735    
# Line 770  void CaloBragg::mediatroncata(){ Line 744  void CaloBragg::mediatroncata(){
744    //    calo[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano    //    calo[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano
745    //    integrale      => energia totale nel calorimetro considerando il W    //    integrale      => energia totale nel calorimetro considerando il W
746    //    //
747    // out[6] chi2,z,Etot,Pskip,flagmaxene,step    // out[4] chi2,z,Etot,Pskip
748    Float_t ordplane[44];//mi serve per la media troncat  
749      Float_t ordplane[44];//mi serve per la media troncata
750    memset(ordplane,0,44*sizeof(Float_t));    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    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    //ordino tutte le energie dei piani in ordine crescente
756    Int_t ii=0;    //Int_t ii=0;
757  //   while( ii < NPLA-1 ){    
758            Long64_t work[200];
759  //     if(ordplane[ii+1] >= ordplane[ii]){      Int_t ind = 0;
760  //       ii++;      Int_t l = 0;
761  //     }else{      Int_t RN = 0;
762  //       SWAP( &(ordplane[ii]), &(ordplane[ii+1]));      Float_t sum4 = 0.;
763  //       ii=0;      Float_t qm = 0.;
764  //     };      //
765  //   };      //Float_t qmt = ethr*0.8; // *0.9
766        //
767        //Int_t uplim = TMath::Max(3,N);
768    //scelgo 4 piani minimo 'sensati'      //
769    Float_t sum4=0.;      while ( l < 4 && ind < 44 ){
770    Int_t pi=0;        qm = TMath::KOrdStat(44,ordplane,ind,work);
771    cout<<"sum4="<<sum4<<endl;        if (qm >= 0.7 ){  
772    while(pi<2*NPLA){          if ( l < 4 ){
773      if ((ordplane[pi] >0.7)&& sum4==0.){// (ordplane[pi] >50)            sum4 += qm;
774        sum4=(ordplane[pi] + ordplane[pi+1] + ordplane[pi+2] + ordplane[pi+3]);            RN++;
775        pi=2*NPLA;          };
776            l++;
777            if ( debug ) printf(" value no %i qm %f sum4 %f \n",l,qm,sum4);
778          };
779          ind++;
780      };      };
781      pi++;      //
782    };      sum4 /= (Float_t)RN;
783    cout<<"sum4="<<sum4<<endl;      Float_t Zmean = (sqrt((sum4*MIP)/(4.*spessore[2])));
784    Float_t Zmean = (sqrt((sum4*MIP)/(4.*spessore[2])));      if(Zmean ==0.) Zmean=1.;
785    cout<<"Zmean="<<Zmean<<endl;  
786    if(Zmean ==0.) Zmean=1.;  //     cout<<"sum4="<<sum4<<endl;
787      //     cout<<"Zmean="<<Zmean<<endl;
   /*z se particella fosse al minimo*/  //energia1piano/mip corretta  
   //  Float_t zmax = (sqrt(estremi[0][1]/spessore[2]));  
788    
789        
790    //calcolo energia migliore per Z trovato con media troncata    //calcolo energia migliore per Z trovato con media troncata
791    //Int_t step = 0;    
792    Float_t zmin=Zmean;    Float_t zmin=Zmean;
793    //  Float_t bestchi2[4] = {10000,0,0,0};//chi2,z,Etot,Pskip //era bestchi2[5]    
794    bestchi2[0]=10000.;    bestchi2[0]=10000.;
795    bestchi2[1]=0.;    bestchi2[1]=0.;
796    bestchi2[2]=0.;    bestchi2[2]=0.;
797    bestchi2[3]=0.;    bestchi2[3]=0.;
798    Float_t zero=0.;    Float_t zero=0.;
799      //    step  energia   zstart  zstop Si attrav  1 piano  energie piani  out  
800      cout<<"inizio media troncata"<<endl;
801      cout<<"input: step "<<Integrale<<", zero "<<zero<<", zmin "<<zmin<<", zmax"<<zmin<<endl;
802    
803    
804        //    step  energia   zstart  zstop  
805    loopze(Integrale,zero,zmin,zmin);    loopze(Integrale,zero,zmin,zmin);
806        
807      cout<<"fine media troncata"<<endl;
808      cout<<"output: chi2 "<<bestchi2[0]<<", z"<<bestchi2[1]<<", Etot "<<bestchi2[2]<<", Pskip "<<bestchi2[3]<<endl;
809    
810    
811    
812    qtchi2=bestchi2[0];    qtchi2=bestchi2[0];
813    qtz=bestchi2[1];    qtz=bestchi2[1];
814    qtetot=bestchi2[2];    qtetot=bestchi2[2];
# Line 839  void CaloBragg::Zdaloop(){ Line 826  void CaloBragg::Zdaloop(){
826    //    calo1[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano    //    calo1[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano
827    //    integrale=> energia totale nel calorimetro considerando il W    //    integrale=> energia totale nel calorimetro considerando il W
828    //    //
829    // out[6] chi2,z,Etot,Pskip    // out[4] chi2,z,Etot,Pskip
830    
831    
832    /*z se particella fosse al minimo*/  //energia1piano/mip corretta    /*z se particella fosse al minimo*/  //energia1piano/mip corretta
833    Float_t zmax = (sqrt(estremi[0][1]/spessore[2]));    Float_t zmax = (sqrt(estremi[0][1]/spessore[2]));
834      if(zmax<31)zmax=zmax+1;
835        
836    /*calcolo Z ed E con loop sui vari elementi ed energie*/    /*calcolo Z ed E con loop sui vari elementi ed energie*/
837    //Int_t step = 0;  
838    Float_t zmin=1.;    Float_t zmin=1.;
839    Float_t bestchitemp[4] = {0,0,0,0};    Float_t bestchitemp[4] = {0,0,0,0};
840    //  Float_t bestchi2[4] = {10000,0,0,0};//chi2,z,Etot,Pskip  
841    bestchi2[0]=10000.;    bestchi2[0]=10000.;
842    bestchi2[1]=0.;    bestchi2[1]=0.;
843    bestchi2[2]=0.;    bestchi2[2]=0.;
# Line 857  void CaloBragg::Zdaloop(){ Line 845  void CaloBragg::Zdaloop(){
845    Float_t zero=0.;    Float_t zero=0.;
846    //primo loop      //primo loop  
847    //     energia   ezero, zstart  zstop Si attrav  1 piano  energie piani  out    //     energia   ezero, zstart  zstop Si attrav  1 piano  energie piani  out
848    
849      cout<<"inizio primo loop"<<endl;
850      cout<<"input: step "<<Integrale<<", zero "<<zero<<", zmin "<<zmin<<", zmax"<<zmax<<endl;
851    
852    loopze(Integrale,zero,zmin,zmax);    loopze(Integrale,zero,zmin,zmax);
853    
854      cout<<"fine primo loop"<<endl;
855      cout<<"output: chi2 "<<bestchi2[0]<<", z"<<bestchi2[1]<<", Etot "<<bestchi2[2]<<", Pskip "<<bestchi2[3]<<endl;
856        
857    //secondo loop    //secondo loop
858    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 866  void CaloBragg::Zdaloop(){ Line 861  void CaloBragg::Zdaloop(){
861    bestchi2[2] = 0.;    bestchi2[2] = 0.;
862    bestchi2[3] = 0.;//riazzero    bestchi2[3] = 0.;//riazzero
863    
864    Float_t step = Integrale/1000.;    Float_t step = bestchitemp[2]/100.;//Integrale/10.;//era 1000
865    zero=bestchitemp[2]-step;//energia da 1 giro - 1step    //zero=bestchitemp[2]-step/2.;//energia da 1 giro - 1step
866      //Float_t step = bestchitemp[2]/1000.;
867      //zero=bestchitemp[2]-step;
868      zero=0;
869    zmin=bestchitemp[1]-2;// zda primo giro    zmin=bestchitemp[1]-2;// zda primo giro
870    if(zmin<1)zmin=1;    if(zmin<1)zmin=1;
871    zmax=bestchitemp[1]+1;    zmax=bestchitemp[1]+1;//
872    
873      cout<<"inizio secondo loop"<<endl;
874      cout<<"input: step "<<step<<", zero "<<zero<<", zmin "<<zmin<<", zmax"<<zmax<<endl;
875    
876      loopze(step,zero,zmin,zmax); //
877    
878      cout<<"fine secondo loop"<<endl;
879      cout<<"output: chi2 "<<bestchi2[0]<<", z"<<bestchi2[1]<<", Etot "<<bestchi2[2]<<", Pskip "<<bestchi2[3]<<endl;
880    
   loopze(step,zero,zmin,zmax);  
881    
882    //  cout<<"z loop ok"<<endl;  //   cout<<"z loop ok"<<endl;
883        
884    //chi2,z,Etot,Pskip    //chi2,z,Etot,Pskip
885    lpchi2=bestchi2[0];    lpchi2=bestchi2[0];
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.23