/[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.6 by mocchiut, Fri Nov 28 16:01:38 2008 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 32  void CaloBragg::Clear(){ Line 32  void CaloBragg::Clear(){
32    //    //
33    tr = 0;    tr = 0;
34    sntr = 0;    sntr = 0;
35    qtchi2 = 0.;  //   qtchi2 = 0.;
36    qtz = 0.;  //   qtz = 0.;
37    qtetot = 0.;  //   qtetot = 0.;
38    qtpskip = 0.;  //   qtpskip = 0.;
39    lpchi2 = 0.;    lpchi2 = 0.;
40    lpz = 0.;    lpz = 0.;
41    lpetot = 0.;    lpetot = 0.;
# Line 44  void CaloBragg::Clear(){ Line 44  void CaloBragg::Clear(){
44    memset(spessore,0,3*sizeof(Float_t));    memset(spessore,0,3*sizeof(Float_t));
45    memset(estremi,0,2*2*sizeof(Float_t));    memset(estremi,0,2*2*sizeof(Float_t));
46    Integrale=0.;    Integrale=0.;
47    
48      for(Int_t l=0;l<44;l++){
49        calorimetro[l][0]=-1.;
50      }
51    
   //  
52  };  };
53    
54  void CaloBragg::Print(){  void CaloBragg::Print(){
# Line 55  void CaloBragg::Print(){ Line 58  void CaloBragg::Print(){
58    //    //
59    printf("========================================================================\n");    printf("========================================================================\n");
60    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);
61    printf(" chi 2 from truncated mean: %f \n", qtchi2);    printf(" first plane: %f \n", estremi[0][0]);
62    printf(" Z from truncated mean %f: \n", qtz);    printf(" last plane: %f \n", estremi[1][0]);
63    printf(" energy from truncated mean %f: \n", qtetot);  //   printf(" chi 2 from truncated mean: %f \n", qtchi2);
64    printf(" plane not used for truncated mean %f: \n", qtpskip);  //   printf(" Z from truncated mean %f: \n", qtz);
65    //   printf(" energy from truncated mean %f: \n", qtetot);
66    //   printf(" plane not used for truncated mean %f: \n", qtpskip);
67    printf(" chi 2 from loop %f:  \n", lpchi2);    printf(" chi 2 from loop %f:  \n", lpchi2);
68    printf(" Z from loop %f: \n", lpz);    printf(" Z from loop %f: \n", lpz);
69    printf(" energy from loop %f: \n", lpetot);    printf(" energy from loop %f: \n", lpetot);
# Line 109  void CaloBragg::Process(Int_t ntr){ Line 114  void CaloBragg::Process(Int_t ntr){
114    
115    //    //
116    //    //
   // Always calculate stdedx1  
117    //    //
118    Int_t view = 0;    Int_t view = 0;
119    Int_t plane = 0;    Int_t plane = 0;
# Line 134  void CaloBragg::Process(Int_t ntr){ Line 138  void CaloBragg::Process(Int_t ntr){
138        ptrack = L2->GetTrack(ntr);        ptrack = L2->GetTrack(ntr);
139        if ( ptrack ) track = ptrack->GetCaloTrack();        if ( ptrack ) track = ptrack->GetCaloTrack();
140      } else {      } else {
141        track = L2->GetCaloStoredTrack(ntr);        track = L2->GetCaloStoredTrack(ntr);  //al momento e' vera solo questa riga
142      };      };
143      //      //
144      if ( !track && ntr >= 0 ){      if ( !track && ntr >= 0 ){
# Line 155  void CaloBragg::Process(Int_t ntr){ Line 159  void CaloBragg::Process(Int_t ntr){
159    //    //
160    for(Int_t p=0; p<22; p++){    for(Int_t p=0; p<22; p++){
161      for(Int_t v=0; v<2; v++){      for(Int_t v=0; v<2; v++){
       cout<< L2->GetCaloLevel2()->cibar[p][v];      
162        /*per usare traccia non del calo camboare cibar*/        /*per usare traccia non del calo camboare cibar*/
163        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
164        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
165      };      };
166    };    };
167    
168      /*per ogni evento calcolo la conversione mip e w attraversato in equivalente Si*/
169      conversione(); // out: 1) g/cm2 Si , 2) spessoreW equivalente in Si, 3)Mip corretta per inclinazione
170    
171    /*per ogni evento clacolo la conversione mip e w attraversato in equivalente Si*/    /*settaggio della soglia per il loop sulla determinazione del piano di partenza */
172      Float_t ordplane[44];//mi serve per la media troncata
173      memset(ordplane,0,44*sizeof(Float_t));
174    
175    conversione(); // out: 1) g/cm2 Si , 2) spessoreW equivalente in Si, 3)Mip corretta per inclinazione    for(Int_t ipla=0; ipla< 2*NPLA; ipla++)  ordplane[ipla]=calorimetro[ipla][1]; //energia del piano
176    
177    cout<<"spessore"<<spessore[0]<<" gcm2si; "<<spessore[1]<<"W in Si; "<<spessore[2]<<"mip "<<endl;    
178      //ordino tutte le energie dei piani in ordine crescente
179    
180        Long64_t work[200];
181        Int_t ind = 0;
182        //Int_t l = 0;
183        Int_t RN = 0;
184        Float_t sum4 = 0.;
185        Float_t qm = 0.;
186        while ( RN < 4 && ind < 44 ){
187          qm = TMath::KOrdStat((Long64_t)44,ordplane,(Long64_t)ind,work);
188          if (qm >= 0.7 ){  
189            if ( RN < 4 ){
190              sum4 += qm;
191              RN++;
192            };
193          };
194          ind++;
195        };
196        //
197        //sum4 /= (Float_t)RN;
198        Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2])));
199        if(Zmean ==0.) Zmean=1.;
200        if ( Zmean < 1. ) Zmean = 1.;
201    
   /*trova primo e ultimo piano attraversati*/  
202    
203      /*trova primo e ultimo piano attraversati*/
204    Int_t p = 0;//contatore piani    Int_t p = 0;//contatore piani
205    //primo parte da 0 e va in giĆ¹    //per il primo parte da 0 e va in giu'
206    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;  
207      //    if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){      //    if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){
208      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
209        if( (calorimetro[p][0] >0) && (calorimetro[p][1]*MIP >Zmean*0.7)){ // 70% della MIP
210        estremi[0][0]=p;        estremi[0][0]=p;
211        estremi[0][1]=calorimetro[p][1] *MIP; //energia in MeV        estremi[0][1]=calorimetro[p][1] *MIP; //energia in MeV
212      };      };
# Line 184  void CaloBragg::Process(Int_t ntr){ Line 214  void CaloBragg::Process(Int_t ntr){
214    };    };
215   //ultimo parte da 44 e sale   //ultimo parte da 44 e sale
216    p=43;    p=43;
217    while( (estremi[1][1] == 0)  &&  (p>(int)estremi[0][0]) ){    while( (estremi[1][1] <= 0.)  &&  (p>(int)estremi[0][0]) ){
218    cout<<"piano "<<p <<"; strip "<<calorimetro[p][0]<<"; en "<<calorimetro[p][1]<<endl;      if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){
219    //    if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){        estremi[1][0]=p;//
   if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){  
       estremi[1][0]=p;  
220        estremi[1][1]=calorimetro[p][1] *MIP;//energia in MeV        estremi[1][1]=calorimetro[p][1] *MIP;//energia in MeV
221      };      };
222      p = p-1;      p = p-1;
223    };    };
224    //    //
   cout<<"estremi: in "<<estremi[0][0]<<"piano, "<<estremi[0][1]<<"energia; out "<<estremi[1][0]<<"piano, "<<estremi[1][1]<<"energia."<<endl;  
   
   
225        
226    /*integrale: energia rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/    /*integrale: energia totale rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/
227      for(Int_t pl=0; pl<(2*NPLA); pl++){
   for(Int_t pl=0; pl<(2*NPLA); pl++){  
   
228      //calcolo intergale in unita di spessori di silicio      //calcolo intergale in unita di spessori di silicio
229      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;  
230      // 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          
231      if(pl%2!=0){                                                             //equival W in Si      if(pl%2!=0){                                                              //equival W in Si
232        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<<Integrale<<" W "<<pl<<endl;  
       cout<<" W "<<0.5*((calorimetro[pl-1][1] * MIP)+(calorimetro[pl][1] * MIP))*(spessore[1])<<endl;  
233      };      };
234    };    };
235      Integrale=24000;//Integrale*1000;
   //  cout<<"integrale ok"<<endl;  
     
     
   /*z se particella fosse al minimo*/  //energia1piano/mip corretta  
   //Float_t zmax = (sqrt(estremi[0][1]/spessore[2]));  
     
   
236        
237    /*z ed energia con media troncata*/    /*z ed energia con media troncata*/
238    //  Float_t  bestchi2mean[4] = {0.,0.,0.,0.};//chi2,z,Etot,Pskip    //  mediatroncata();  // out:  1)chi2, 2)z, 3)Etot, 4)Pskip
   mediatroncata();    
239        
240    /*z ed energia con loop*/    /*z ed energia con loop*/
241    //  Float_t  bestchi2loop[4] = {0.,0.,0.,0.};//chi2,z,Etot,Pskip//chi2,z,Etot,Pskip      Zdaloop(); // out:  1)chi2, 2)z, 3)Etot, 4)Pskip
   Zdaloop();  
     
   
   /*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  
242    
243        
244    if ( debug ) this->Print();    if ( debug ) this->Print();
245    if ( debug ) printf(" esci \n");    if ( debug ) printf(" fine evento \n");
246    //    //
247  };  };
248    
# Line 251  void CaloBragg::Draw(){ Line 251  void CaloBragg::Draw(){
251    
252    Process();    Process();
253    
254   Float_t dEpianimean[44];  //  Float_t dEpianimean[44];
255   Float_t dEpianiloop[44];   Float_t dEpianiloop[44];
256   Float_t Depth[44];   Float_t Depth[44];
257   Int_t tz=(Int_t)qtz;  //  Int_t tz=(Int_t)qtz;
258   Int_t tz1=(Int_t)lpz;   Int_t tz1=(Int_t)lpz;
259   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
260   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
261    
262   Float_t sp= spessore[0]*spessore[1];   Float_t sp= spessore[0]*spessore[1];
263   for(Int_t i=0;i<44;i++)Depth[i]=i*sp;   for(Int_t i=0;i<44;i++)Depth[i]=i*sp;
264    //    //
265    gStyle->SetLabelSize(0.04);    gStyle->SetLabelSize(0.04);
266    gStyle->SetNdivisions(510,"XY");    gStyle->SetNdivisions(510,"XY");
   //gStyle->SetLogy();  
267    //    //
   
268          TString hid = Form("cCaloBragg");                TString hid = Form("cCaloBragg");      
269          TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));          TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
   
270          if ( tc ){          if ( tc ){
271  //       tc->Clear();  //       tc->Clear();
272          } else {          } else {
273           tc = new TCanvas(hid,hid);           tc = new TCanvas(hid,hid);
274           tc->Divide(1,2);  //       tc->Divide(1,2);
275          };          };
276          //          //
277          TString thid = Form("hCaloBragg");        //      TString thid = Form("hCaloBragg");      
278          TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));  //         TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
279          if ( th ) th->Delete();  //      if ( th ) th->Delete();
280  //       th->Clear();  //       th->Clear();
281  //       th->Reset();  //       th->Reset();
282  //      } else {  //      } else {
283          th = new TH2F(thid,thid,300,-0.5,300.,1000,0.,150.);  //      th = new TH2F(thid,thid,300,-0.5,300.,1000,0.,150.);
284          th->SetMarkerStyle(20);  //      th->SetMarkerStyle(20);
285  //      };  //      };
286          //          //
287          TString thid2 = Form("hCaloBragg2");              TString thid2 = Form("hCaloBragg2");    
# Line 304  void CaloBragg::Draw(){ Line 301  void CaloBragg::Draw(){
301    
302          tc->cd(1);          tc->cd(1);
303          //          //
304          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]);
305          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);
306          th->Draw();  //      th->Draw();
307          th2->Draw("same");          th2->Draw("same");
308    
309          tc->cd(2);          tc->cd(2);
310          //          //
311          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]);
312          th3->Draw();          th3->Draw();
313          th2->Draw("same");          th2->Draw("same");
314    
# Line 325  void CaloBragg::Draw(){ Line 322  void CaloBragg::Draw(){
322  };  };
323    
324    
 void CaloBragg::SWAP( Float_t *A, Float_t *B ){  
    Float_t  Tmp = *A;  
    *A = *B;  
    *B = Tmp;  
 };  
   
   
325    
326  void CaloBragg::LoadParam(){  void CaloBragg::LoadParam(){
327    
# Line 370  void CaloBragg::LoadParam(){ Line 360  void CaloBragg::LoadParam(){
360    elem[31] = 72.61; //Ge 32    elem[31] = 72.61; //Ge 32
361    
362    
   
363  //parametri calorimetro  //parametri calorimetro
364    NPLA = 22;    NPLA = 22;
365    NCHA = 96;    NCHA = 96;
# Line 420  void CaloBragg::conversione(){ Line 409  void CaloBragg::conversione(){
409    Float_t ytgy = 0;    Float_t ytgy = 0;
410    Float_t a = 0.;    Float_t a = 0.;
411    
   
412    /*silicio*/    /*silicio*/
413    ytgx = ySi * L2->GetCaloLevel2()->tanx[0];    ytgx = ySi * L2->GetCaloLevel2()->tanx[0];
414    ytgy = ySi * L2->GetCaloLevel2()->tany[0];    ytgy = ySi * L2->GetCaloLevel2()->tany[0];
# Line 430  void CaloBragg::conversione(){ Line 418  void CaloBragg::conversione(){
418        
419    spessore[0] = SiCross/10. * rhoSi; //spessore silicio in g/cm2    spessore[0] = SiCross/10. * rhoSi; //spessore silicio in g/cm2
420        
     
421    /*tungsteno*/    /*tungsteno*/
422    ytgx = yW * L2->GetCaloLevel2()->tanx[0];    ytgx = yW * L2->GetCaloLevel2()->tanx[0];
423    ytgy = yW * L2->GetCaloLevel2()->tany[0];    ytgy = yW * L2->GetCaloLevel2()->tany[0];
# Line 444  void CaloBragg::conversione(){ Line 431  void CaloBragg::conversione(){
431    //       (g/cm2W)/(g/cm2Si)    //       (g/cm2W)/(g/cm2Si)
432    spessore[1] =  a;    spessore[1] =  a;
433    
     
434    //riscala mip allo spessore attraversato    //riscala mip allo spessore attraversato
435    spessore[2] = MIP*(SiCross/ySi);      spessore[2] = MIP*(SiCross/ySi);  
436    
# Line 472  void CaloBragg::BetheBloch(Float_t *x, F Line 458  void CaloBragg::BetheBloch(Float_t *x, F
458    Float_t lg =0.;    Float_t lg =0.;
459    Float_t Energia=0.;    Float_t Energia=0.;
460    Float_t C=0.;    Float_t C=0.;
   
461    
462    eta = (*gam)*(*Bet);    eta = (*gam)*(*Bet);
463    
# Line 481  void CaloBragg::BetheBloch(Float_t *x, F Line 466  void CaloBragg::BetheBloch(Float_t *x, F
466        
467    lg = 2.* Me * SQ(eta) * Wmax / SQ(ISi);    lg = 2.* Me * SQ(eta) * Wmax / SQ(ISi);
468                                                                                                 //  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;
469    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.) +
470      (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.);
471        
472    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.);  
     
473    
474    /*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.);  
475    
476    *out =Energia;//out    *out =Energia;//out
477    
# Line 510  void CaloBragg::ELOSS(Float_t *dx, Int_t Line 491  void CaloBragg::ELOSS(Float_t *dx, Int_t
491    
492    
493    Float_t Q=0.;    Float_t Q=0.;
494    Float_t v=0.;    Float_t v=0.;                                                        
                                                         //  Float_t Mass=0.;  
495    Float_t gam=0.;    Float_t gam=0.;
496    Float_t Bet=0.;    Float_t Bet=0.;
497    Float_t dEP=0.;    Float_t dEP=0.;
# Line 529  void CaloBragg::ELOSS(Float_t *dx, Int_t Line 509  void CaloBragg::ELOSS(Float_t *dx, Int_t
509    
510    //perdita energia per un protone    //perdita energia per un protone
511    Float_t protone =1.;    Float_t protone =1.;
512    Float_t Mass=MassP;    Float_t Mass=(elem[*Z-1]*MassP);
513    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);
514        
515    *out= (SQ(Q)*(dEP));//*dx;    *out= (SQ(Q)*(dEP));//*dx;
516    
 };//end ELOSS  
   
517    
518    };//end ELOSS
519    
520    
521    
522    
523  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[]){
524    
525    //calcola energia rilasciata sulla traccia (usa ELOSS)    //calcola energia rilasciata sulla traccia (usa ELOSS)
526    //  in : Z             =>carica    //  in : Z             =>carica
# Line 560  void CaloBragg::Enetrack(Int_t* Z, Float Line 539  void CaloBragg::Enetrack(Int_t* Z, Float
539    memset(out, 0, 2*NPLA*sizeof(Float_t));    memset(out, 0, 2*NPLA*sizeof(Float_t));
540        
541    Float_t Massa = (elem[(*Z)-1] * MassP);    Float_t Massa = (elem[(*Z)-1] * MassP);
     
   
542    
543    //loop piani      (dal primo in cui entra)    for( Int_t ipla=((int)(*primo)); ipla<= ((int)(*ultimo)); ipla++){
   for( Int_t ipla=((int)(*primo)); ipla< (2*NPLA); ipla++){  
544      dE=0.;      dE=0.;
   
545      //spessore silicio corretto x inclinazione, z, energia, out:rilascio      //spessore silicio corretto x inclinazione, z, energia, out:rilascio
546      ELOSS(&spessore[0], Z, &Ezero, &dE);//spessore in g/cm2!!      ELOSS(&spessore[0], Z, &Ezero, &dE);//spessore in g/cm2!!
       
547      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
548        out[ipla] = Ezero - Massa; //MeV        out[ipla] = Ezero - Massa; //MeV
549        return;        return;
# Line 578  void CaloBragg::Enetrack(Int_t* Z, Float Line 552  void CaloBragg::Enetrack(Int_t* Z, Float
552        out[ipla] = dE; //MeV        out[ipla] = dE; //MeV
553        Ezero = Ezero - dE;//energia residua        Ezero = Ezero - dE;//energia residua
554      };      };
   
555      //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
556      if(ipla%2 == 0){      if(ipla%2 == 0){
557        /*tungsteno*/        /*tungsteno*/
558        dE=0.;        dE=0.;
559        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
560        ELOSS(&sp, Z, &Ezero, &dE);        ELOSS(&sp, Z, &Ezero, &dE);
         
561        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
562          return;          return;
563        }else{        }else{
564          Ezero = Ezero -dE;//energia residua          Ezero = Ezero -dE;//energia residua
         //      cout<<"w calc "<<dE<<endl;  
565        };        };
566      };          };    
567    
568    };//fine loop piani    };//fine loop piani
     
 };//end Enetrack  
569    
570    
571    };//end Enetrack
572    
573    
574    
# Line 613  void CaloBragg::chiquadro(Float_t dE[], Line 584  void CaloBragg::chiquadro(Float_t dE[],
584    
585    Float_t sum = 0.;    Float_t sum = 0.;
586    Float_t PianoPrecedente=0.;    Float_t PianoPrecedente=0.;
   //  Float_t differenza =0.;  
587    Float_t badplane=0.;    Float_t badplane=0.;
588    Float_t badplanetot=0.;    Float_t badplanetot=0.;
589    Float_t w,wi;    Float_t w,wi;
590    
591    for(Int_t ipla=0; ipla<2*NPLA; ipla++){    for(Int_t ipla=0; ipla<2*NPLA; ipla++){
     
592      //tutti i piani attraversati dalla traiettoria      //tutti i piani attraversati dalla traiettoria
593      if(calorimetro[ipla][0] != -1.){      if(calorimetro[ipla][0] != -1.){ //
594        w=0.; //normalizzazione;        w=0.; //normalizzazione;
595        wi=1.;//peso        wi=1.;//peso
596    
# Line 629  void CaloBragg::chiquadro(Float_t dE[], Line 598  void CaloBragg::chiquadro(Float_t dE[],
598        if (ipla<estremi[0][0])  wi=0.;        if (ipla<estremi[0][0])  wi=0.;
599                
600        //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
601        //    if((ipla>estremi[1][0]) && (calorimetro[ipla][1] >0.) ) wi=0.;        //if((ipla>estremi[1][0]) && (calorimetro[ipla][1] >0.) ) wi=0.;
602          if((ipla>estremi[1][0])) wi=0.;
603    
604        //normalizzazione        //normalizzazione
605        if (calorimetro[ipla][1] != 0.)  w=1./(calorimetro[ipla][1]* MIP);                if (calorimetro[ipla][1] != 0.)  w=1./(calorimetro[ipla][1]* MIP);    //  
606                
607        //tolgo piani con rilasci inferiori al 30% del precedente        //tolgo piani con rilasci inferiori al 30% del precedente
608        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
           
609          wi=0.;          wi=0.;
610          //se sono piani intermedi (non si ĆØ fermta) li considero non buoni          //se sono piani intermedi (non si e' fermta) li considero non buoni
611          if( (ipla < estremi[1][0]) && (calorimetro[ipla][1] !=0.)){          if( (ipla <= estremi[1][0]) && (calorimetro[ipla][1] !=0.)){//
612            badplane+=1.;            badplane+=1.;
613            badplanetot+=1.;            badplanetot+=1.;
614          };          };
615        };        };
         
   
616        
617        //meno peso ai piani con rilasci maggiori di 1000 MIP        //meno peso ai piani con rilasci maggiori di 1000 MIP
618          if(calorimetro[ipla][1] > 1000) wi=0.5;          if(calorimetro[ipla][1] > 1000) wi=0.5;
619                            
         //do peso maggiore alle ultime 6 misure  
         //      if((ipla >= estr2[1][0]-12)) wi=0.;  
   
620        Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));        Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));
621    
622        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));
623          if(debug){
624            printf("dedx  calcolata  %f e reale  %f  \n",dE[ipla],(calorimetro[ipla][1] * MIP));
625          }
626        //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
627        if(wi != 0.){        if(wi != 0.){//
628          PianoPrecedente= calorimetro[ipla][1];//tengo piano precedente          PianoPrecedente= calorimetro[ipla][1];//tengo piano precedente
629          badplane = 0.;//azzero contatore piani scartati consecutivi          badplane = 0.;//azzero contatore piani scartati consecutivi
630        };        };
631      };      };
632                        
633      if(badplane > 2) out[1] =79.;      //da Emi
634                  if(badplane > 2){      
635          out[1] =79.;
636          break;
637        };
638      
639    };//fine loop piani    };//fine loop piani
640    //chi2,frammentato,pskip      //chi2,frammentato,pskip  
641    out[0]=sum;    out[0]=sum;
   
642    out[2]=badplanetot;    out[2]=badplanetot;
643        
644  };//end chiquadro  };//end chiquadro
645    
646    
647    
 //loopze(Integrale,&zero,&zmin, &zmax, spess1 , estr3, calo1, bestchi2);  
648  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){
649  //  //
650    //loop su z ed energie per trovare miglior z (ed energia)    //loop su z ed energie per trovare miglior z (ed energia)
# Line 684  void CaloBragg::loopze( Float_t step, Fl Line 652  void CaloBragg::loopze( Float_t step, Fl
652    //     E0        => energia iniziale (intergale)            //     E0        => energia iniziale (intergale)        
653    //     Zstart    => minimo z da cui patire    //     Zstart    => minimo z da cui patire
654    //     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  
655    //    //
656    //out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2    //out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2
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)/1000.;//passo per il calcolo di energia
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++){// da non cambiare in base a Stepint altrimenti cambia la matrice bestchi2!!!cosi' non raggiungo mai integrale!!!!! mettere <=??
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);
699                
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            break;
709        }        }
710      }//fine loop energia      }//fine loop energia
711    
712    
713   };//fine loop z   };//fine loop z
714    
715    for (Int_t nu=0; nu<colmax; nu++){  
716      for (Int_t en=0; en<rowmax; en++){    //Emi
717        if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] !=0.)){    for (Int_t nu=(int)(Zstart); nu<max; nu++){
718        for (Int_t en=0; en<1000; en++){
719          if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] >0.)){
720          bestchi2[0]= matrixchi2[nu][en][0];// chi2          bestchi2[0]= matrixchi2[nu][en][0];// chi2
721          bestchi2[1]= (Float_t)nu;          bestchi2[1]= (Float_t)nu; // z
722          bestchi2[2]= matrixchi2[nu][en][1];//energia;          bestchi2[2]= matrixchi2[nu][en][1];//energia;
723          bestchi2[3]=    matrixchi2[nu][en][2];// totale piani saltati          bestchi2[3]= matrixchi2[nu][en][2];// totale piani saltati
724        }        }
725      }      }
726    }    }
727    
 };//endloopze  
728    
729    };//endloopze
730    
731    
732    
733    
 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[6] chi2,z,Etot,Pskip,flagmaxene,step  
   Float_t ordplane[44];//mi serve per la media troncat  
   memset(ordplane,0,44*sizeof(Float_t));  
734    
735    for(Int_t ipla=0; ipla< 2*NPLA; ipla++)  ordplane[ipla]=calorimetro[ipla][1]; //energia del piano  // void CaloBragg::mediatroncata(){
736    //   //calcolo Z con media troncata e utilizzo questo Z per trovare l'energia migliore
737    //   //in: ordplane[44]   => array con energia dei piani
738    //   //    spess[3]    => conversioni spessore di silicio, w,  mip
739    //   //    estr[2][2]  => primo[0][0] e ultimo[1][0] piano attraversati ed energie[][1]
740    //   //    calo[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano
741    //   //    integrale      => energia totale nel calorimetro considerando il W
742    //   //
743    //   // out[4] chi2,z,Etot,Pskip
744    
745    //   Float_t ordplane[44];//mi serve per la media troncata
746    //   memset(ordplane,0,44*sizeof(Float_t));
747    
748    //   for(Int_t ipla=0; ipla< 2*NPLA; ipla++)  ordplane[ipla]=calorimetro[ipla][1]; //energia del piano
749    
750      
751    //   //ordino tutte le energie dei piani in ordine crescente
752    
753    //     Long64_t work[200];
754    //     Int_t ind = 0;
755    //     //Int_t l = 0;
756    //     Int_t RN = 0;
757    //     Float_t sum4 = 0.;
758    //     Float_t qm = 0.;
759    //     //
760    //     //Float_t qmt = ethr*0.8; // *0.9
761    //     //
762    //     //Int_t uplim = TMath::Max(3,N);
763    //     //
764    //     while ( RN < 4 && ind < 44 ){
765    //       qm = TMath::KOrdStat(44,ordplane,ind,work);
766    //       if (qm >= 0.7 ){      
767    //      if ( RN < 4 ){
768    //        sum4 += qm;
769    //        RN++;
770    //      };
771    // //   l++;
772    // //   if ( debug ) printf(" value no %i qm %f sum4 %f \n",l,qm,sum4);
773    //       };
774    //       ind++;
775    //     };
776    //     //
777    //     sum4 /= (Float_t)RN;
778    //     Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2])));//ma non e'/1??
779    //     if(Zmean ==0.) Zmean=1.;
780    //     if ( Zmean < 1. ) Zmean = 1.;
781    
     
   //ordino tutte le energie dei piani in ordine crescente  
   Int_t ii=0;  
 //   while( ii < NPLA-1 ){    
782            
783  //     if(ordplane[ii+1] >= ordplane[ii]){  // //     Zmean =round(Zmean);
784  //       ii++;  // //     if(Zmean <1.) Zmean=1.;
 //     }else{  
 //       SWAP( &(ordplane[ii]), &(ordplane[ii+1]));  
 //       ii=0;  
 //     };  
 //   };  
785    
786    // //     if(Zmean >0.)Zmean =round(Zmean);
787    
788    //scelgo 4 piani minimo 'sensati'  //     //======== per i nuclei=======
789    Float_t sum4=0.;  //     if (Zmean >=2.){
790    Int_t pi=0;  //     ind = 0;
791    cout<<"sum4="<<sum4<<endl;  //     RN = 0;
792    while(pi<2*NPLA){  //     sum4 = 0.;
793      if ((ordplane[pi] >0.7)&& sum4==0.){// (ordplane[pi] >50)  //     qm = 0.;
794        sum4=(ordplane[pi] + ordplane[pi+1] + ordplane[pi+2] + ordplane[pi+3]);  //     while ( RN < 4 && ind < 44 ){
795        pi=2*NPLA;  //       qm = TMath::KOrdStat(44,ordplane,ind,work);
796      };  //       if (qm >= (Zmean*Zmean)-Zmean*Zmean*0.2 ){    
797      pi++;  //      if ( RN < 4 ){
798    };  //        sum4 += qm;
799    cout<<"sum4="<<sum4<<endl;  //        RN++;
800    Float_t Zmean = (sqrt((sum4*MIP)/(4.*spessore[2])));  //      };
801    cout<<"Zmean="<<Zmean<<endl;  //       };
802    if(Zmean ==0.) Zmean=1.;  //       ind++;
803      //     };
804    /*z se particella fosse al minimo*/  //energia1piano/mip corretta  //     //
805    //  Float_t zmax = (sqrt(estremi[0][1]/spessore[2]));  //     sum4 /= (Float_t)RN;
806    //     Zmean = (sqrt((sum4*MIP)/(4.*spessore[2])));//ma non e' /1??
807      //     }
808    //calcolo energia migliore per Z trovato con media troncata  
809    //Int_t step = 0;    
810    Float_t zmin=Zmean;  //   //calcolo energia migliore per Z trovato con media troncata
811    //  Float_t bestchi2[4] = {10000,0,0,0};//chi2,z,Etot,Pskip //era bestchi2[5]  //   //  Float_t zmin=Zmean;
812    bestchi2[0]=10000.;  //     Float_t zmin=round(Zmean);
813    bestchi2[1]=0.;    
814    bestchi2[2]=0.;  //   bestchi2[0]=10000.;
815    bestchi2[3]=0.;  //   bestchi2[1]=0.;
816    Float_t zero=0.;  //   bestchi2[2]=0.;
817      //    step  energia   zstart  zstop Si attrav  1 piano  energie piani  out  //   bestchi2[3]=0.;
818    loopze(Integrale,zero,zmin,zmin);  //   Float_t zero=0.;
819      
820    qtchi2=bestchi2[0];  //   //    step   energia zstart zstop  
821    qtz=bestchi2[1];  //   loopze(Integrale,zero,zmin,zmin);
822    qtetot=bestchi2[2];    
823    qtpskip=bestchi2[3];  
824  };//end mediatroncata  //   qtchi2=bestchi2[0];
825    //  cout<<"z media troncata ok"<<endl;  //   qtz=bestchi2[1];
826    //   qtetot=bestchi2[2];
827    //   qtpskip=bestchi2[3];
828    // };//end mediatroncata
829    
830    
831    
# Line 839  void CaloBragg::Zdaloop(){ Line 837  void CaloBragg::Zdaloop(){
837    //    calo1[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano    //    calo1[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano
838    //    integrale=> energia totale nel calorimetro considerando il W    //    integrale=> energia totale nel calorimetro considerando il W
839    //    //
840    // out[6] chi2,z,Etot,Pskip    // out[4] chi2,z,Etot,Pskip
841    
842    
843    /*z se particella fosse al minimo*/  //energia1piano/mip corretta    /*z se particella fosse al minimo*/  //energia1piano/mip corretta
844    Float_t zmax = (sqrt(estremi[0][1]/spessore[2]));    Float_t zmax = round(sqrt(estremi[0][1]/spessore[2]));
845      if(zmax<31)zmax=zmax+1;
846        
847    /*calcolo Z ed E con loop sui vari elementi ed energie*/    /*calcolo Z ed E con loop sui vari elementi ed energie*/
848    //Int_t step = 0;  
849    Float_t zmin=1.;    Float_t zmin=1.;
850    Float_t bestchitemp[4] = {0,0,0,0};    Float_t bestchitemp[4] = {0,0,0,0};
851    //  Float_t bestchi2[4] = {10000,0,0,0};//chi2,z,Etot,Pskip  
852    bestchi2[0]=10000.;    bestchi2[0]=10000.;
853    bestchi2[1]=0.;    bestchi2[1]=0.;
854    bestchi2[2]=0.;    bestchi2[2]=0.;
855    bestchi2[3]=0.;    bestchi2[3]=0.;
856    Float_t zero=0.;    Float_t zero=0.;
857    //primo loop      //------------primo loop   ----------------------
858    //     energia   ezero, zstart  zstop Si attrav  1 piano  energie piani  out    //     energia   ezero, zstart  zstop
859    loopze(Integrale,zero,zmin,zmax);    loopze(Integrale,zero,zmin,zmax);
860    
861        
862    //secondo loop    //------------secondo loop  ----------------------
863    for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i];    for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i];
864    bestchi2[0] = 10000.;    bestchi2[0] = 10000.;
865    bestchi2[1] = 0.;    bestchi2[1] = 0.;
866    bestchi2[2] = 0.;    bestchi2[2] = 0.;
867    bestchi2[3] = 0.;//riazzero    bestchi2[3] = 0.;//riazzero
868    
869    Float_t step = Integrale/1000.;    Float_t step = bestchitemp[2];//
870    zero=bestchitemp[2]-step;//energia da 1 giro - 1step    zero=0;  // qualsiasi altro valore peggiora le cose
871    zmin=bestchitemp[1]-2;// zda primo giro    zmin=zmax=bestchitemp[1];
872    if(zmin<1)zmin=1;    loopze(step,zero,zmin,zmax); //
   zmax=bestchitemp[1]+1;  
   
   loopze(step,zero,zmin,zmax);  
873    
   //  cout<<"z loop ok"<<endl;  
874        
875    //chi2,z,Etot,Pskip    //chi2,z,Etot,Pskip
876    lpchi2=bestchi2[0];    lpchi2=bestchi2[0];

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.6

  ViewVC Help
Powered by ViewVC 1.1.23