/[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.2 by pamelats, Wed Jun 18 15:24:12 2008 UTC revision 1.13 by mocchiut, Thu Oct 6 09:27:39 2011 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      newchi2 = false;
30      usenewBB = false;
31      fzeta = -1.;
32    //    //
33  };  };
34    
35  void CaloBragg::Clear(){  void CaloBragg::Clear(){
36    //    //
37      ndf = 0;
38    tr = 0;    tr = 0;
39    sntr = 0;    sntr = 0;
40    qtchi2 = 0.;    //   qtchi2 = 0.;
41    qtz = 0.;    //   qtz = 0.;
42    qtetot = 0.;    //   qtetot = 0.;
43    qtpskip = 0.;    //   qtpskip = 0.;
44    lpchi2 = 0.;    lpchi2 = 0.;
45    lpz = 0.;    lpz = 0.;
46      lpisotope= 0.;
47    lpetot = 0.;    lpetot = 0.;
48    lppskip = 0.;    lppskip = 0.;
49    
50    memset(calorimetro,0,44*2*sizeof(Float_t));    memset(calorimetro,0,44*2*sizeof(Float_t));
51    memset(spessore,0,3*sizeof(Float_t));    memset(spessore,0,4*sizeof(Float_t));
52    memset(estremi,0,2*2*sizeof(Float_t));    memset(estremi,0,2*2*sizeof(Float_t));
53    Integrale=0.;    Integrale=0.;
54    
55      for(Int_t l=0;l<44;l++){
56        calorimetro[l][0]=-1.;
57      }
58    
   //  
59  };  };
60    
61  void CaloBragg::Print(){  void CaloBragg::Print(){
# Line 55  void CaloBragg::Print(){ Line 65  void CaloBragg::Print(){
65    //    //
66    printf("========================================================================\n");    printf("========================================================================\n");
67    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);
68    printf(" chi 2 from truncated mean: %f \n", qtchi2);    printf(" first plane: %f \n", estremi[0][0]);
69    printf(" Z from truncated mean %f: \n", qtz);    printf(" last plane: %f \n", estremi[1][0]);
70    printf(" energy from truncated mean %f: \n", qtetot);    //   printf(" chi 2 from truncated mean: %f \n", qtchi2);
71    printf(" plane not used for truncated mean %f: \n", qtpskip);    //   printf(" Z from truncated mean %f: \n", qtz);
72      //   printf(" energy from truncated mean %f: \n", qtetot);
73      //   printf(" plane not used for truncated mean %f: \n", qtpskip);
74    printf(" chi 2 from loop %f:  \n", lpchi2);    printf(" chi 2 from loop %f:  \n", lpchi2);
75    printf(" Z from loop %f: \n", lpz);    printf(" Z from loop %f: \n", lpz);
76      printf(" isotope from loop %f: \n", lpisotope);
77    printf(" energy from loop %f: \n", lpetot);    printf(" energy from loop %f: \n", lpetot);
78    printf(" plane not used for loop %f: \n", lppskip);    printf(" plane not used for loop %f: \n", lppskip);
79      printf(" ndf: %i \n",ndf);
80    printf("========================================================================\n");    printf("========================================================================\n");
81    //    //
82  };  };
# Line 77  void CaloBragg::Process(){ Line 91  void CaloBragg::Process(){
91    Process(-1);    Process(-1);
92  };  };
93    
94    
95    void CaloBragg::CleanPlanes(Float_t epiano[22][2]){
96      //  return;
97      Int_t hitplanes = 0;
98      for (Int_t i = 0; i<22; i++){
99        for (Int_t j = 1; j>=0; j--){
100          if ( epiano[i][j] > 0.7 ) hitplanes++;
101        };
102      };
103      Float_t lowlim = 0.85;
104      Float_t dedxone = 0.;
105      Float_t step1 = 0.8*L2->GetCaloLevel2()->qtot/(Float_t)hitplanes;
106      while ( dedxone < step1 ){
107        for (Int_t i = 0; i<22; i++){
108          for (Int_t j = 1; j>=0; j--){
109            if ( epiano[i][j] >= step1 && dedxone < 0.7 ) dedxone = epiano[i][j];
110          };
111        };    
112      }
113      if ( dedxone < 0.7 ){
114        for (Int_t i = 0; i<22; i++){
115          for (Int_t j = 1; j>=0; j--){
116            if ( epiano[i][j] > 0. && dedxone < 0.7 ) dedxone = epiano[i][j];
117          };
118        };    
119      }
120      //
121      //  printf(" dedxone = %f step1 %f  \n",dedxone,step1);
122      Bool_t revulsera = false;
123      Bool_t nullius = false;
124      Int_t nulliferus = 0;
125      for (Int_t i = 0; i<22; i++){
126        for (Int_t j = 1; j>=0; j--){
127          if ( epiano[i][j] < dedxone*lowlim ){
128            //        printf(" %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius);
129            epiano[i][j] = 0.;
130          } else {
131            //x        printf(" %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius);
132            nulliferus = 0;
133            revulsera = true;
134          };
135          if ( epiano[i][j] < 0.7 && revulsera ) nulliferus++;
136          if ( nulliferus > 10 ) nullius = true;
137          if ( nullius ) epiano[i][j] = 0.;
138        };
139      };
140    
141    }
142    
143    
144  void CaloBragg::Process(Int_t ntr){  void CaloBragg::Process(Int_t ntr){
145    //      //  
146    if ( !L2 ){    if ( !L2 ){
# Line 119  void CaloBragg::Process(Int_t ntr){ Line 183  void CaloBragg::Process(Int_t ntr){
183    for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){    for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
184      //      //
185      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
186        //
187        if ( !usepl18x && view==0 && plane==18 ) mip = 0.;
188        //
189      epiano[plane][view]+=mip;      epiano[plane][view]+=mip;
190      //      //
191      //      //
192    };    };
193    //    //
194      this->CleanPlanes(*&epiano);
195    //    //
196    PamTrack *ptrack = 0;    PamTrack *ptrack = 0;
197    CaloTrkVar *track = 0;    CaloTrkVar *track = 0;
# Line 155  void CaloBragg::Process(Int_t ntr){ Line 223  void CaloBragg::Process(Int_t ntr){
223    for(Int_t p=0; p<22; p++){    for(Int_t p=0; p<22; p++){
224      for(Int_t v=0; v<2; v++){      for(Int_t v=0; v<2; v++){
225        /*per usare traccia non del calo camboare cibar*/        /*per usare traccia non del calo camboare cibar*/
226        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      
227        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
228      };      };
229    };    };
230    
231    /*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*/
232    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;  
233    
234      /*settaggio della soglia per il loop sulla determinazione del piano di partenza */
235      Float_t ordplane[44];//mi serve per la media troncata
236      memset(ordplane,0,44*sizeof(Float_t));
237    
238      for(Int_t ipla=0; ipla< 2*NPLA; ipla++)  ordplane[ipla]=calorimetro[ipla][1]; //energia del piano
239    
240      
241      //ordino tutte le energie dei piani in ordine crescente
242    
243      Long64_t work[200];
244      Int_t ind = 0;
245      //Int_t l = 0;
246      Int_t RN = 0;
247      Float_t sum4 = 0.;
248      Float_t qm = 0.;
249      while ( RN < 4 && ind < 44 ){
250        qm = TMath::KOrdStat((Long64_t)44,ordplane,(Long64_t)ind,work);
251        if (qm >= 0.7 ){    
252          if ( RN < 4 ){
253            sum4 += qm;
254            RN++;
255          };
256        };
257        ind++;
258      };
259      //
260      //sum4 /= (Float_t)RN;
261      Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2])));
262      if(Zmean ==0.) Zmean=1.;
263      if ( Zmean < 1. ) Zmean = 1.;
264    
265    
266    /*trova primo e ultimo piano attraversati*/    /*trova primo e ultimo piano attraversati*/
267    Int_t p = 0;//contatore piani    Int_t p = 0;//contatore piani
268    //per il primo parte da 0 e va in giù    //per il primo parte da 0 e va in giu'
269    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;  
270      //    if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){      //    if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){
271      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
272        if( (calorimetro[p][0] >0) && (calorimetro[p][1]*MIP >Zmean*0.7)){ // 70% della MIP
273        estremi[0][0]=p;        estremi[0][0]=p;
274        estremi[0][1]=calorimetro[p][1] *MIP; //energia in MeV        estremi[0][1]=calorimetro[p][1] *MIP; //energia in MeV
275      };      };
276      p++;      p++;
277    };    };
278   //ultimo parte da 44 e sale  
279      //ultimo parte da 44 e sale
280    p=43;    p=43;
281    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.)){  
282      if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){      if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){
283        estremi[1][0]=p;//era p        estremi[1][0]=p;//
284        estremi[1][1]=calorimetro[p][1] *MIP;//energia in MeV        estremi[1][1]=calorimetro[p][1] *MIP;//energia in MeV
285      };      };
286      p = p-1;      p = p-1;
287    };    };
288    //    //
289    
290      Float_t lastok = 0.;
291      //  Bool_t goback = false;
292      for ( int o = 0; o < estremi[1][0]; o++ ){
293        //
294        if ( calorimetro[o][1] > 0.7 ) lastok = calorimetro[o][1];
295        if ( calorimetro[o][1] < 0.7 && lastok > 0. ) calorimetro[o][1] = lastok;    
296        //    if ( calorimetro[o][1] < 0.7 ) goback = true;
297        //
298      };
299      lastok = 0.;
300      //  if ( goback ){
301      for ( int o = estremi[1][0]; o >= 0;  o-- ){
302        //
303        //    printf(" goback1: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);
304        if ( o < estremi[1][0] && calorimetro[o][1] > calorimetro[o+1][1]*1.2 && lastok > 0. ) calorimetro[o][1] = lastok;  
305        if ( calorimetro[o][1] > 0.7 ) lastok = calorimetro[o][1];
306        if ( calorimetro[o][1] < 0.7 && lastok > 0. ) calorimetro[o][1] = lastok;    
307        //    printf(" goback2: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);
308        //
309      };
310      //  };
311    
312      if ( startZero ) {
313        estremi[0][0] = 0.;
314        //    estremi[0][1] = 0.;
315      }
316        
317    /*integrale: energia totale rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/    /*integrale: energia totale rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/
318    for(Int_t pl=0; pl<(2*NPLA); pl++){    for(Int_t pl=0; pl<(2*NPLA); pl++){
319        //    printf(" integrale: calorimetro %f  \n",calorimetro[pl][1]);
320      //calcolo intergale in unita di spessori di silicio      //calcolo intergale in unita di spessori di silicio
321      Integrale += calorimetro[pl][1] * MIP;//piano di silicio      Integrale += calorimetro[pl][1] * MIP;//piano di silicio
322      // 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          
323      if(pl%2!=0){                                                              //equival W in Si      if(pl%2!=0){                                                              //equival W in Si
324        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;  
325      };      };
326    };    };
327        //Integrale=24000;//Integrale*1000;
328      Integrale *= 1000.;
329        
330    /*z ed energia con media troncata*/    /*z ed energia con media troncata*/
331  //   cout<<"Media troncata"<<endl;    //  mediatroncata();  // out:  1)chi2, 2)z, 3)Etot, 4)Pskip
   mediatroncata();  // out:  1)chi2, 2)z, 3)Etot, 4)Pskip  
332        
333    /*z ed energia con loop*/    /*z ed energia con loop*/
334  //    cout<<"Zdaloop"<<endl;    if ( debug ) printf(" call Zdaloop with integrale %f \n",Integrale);
335    Zdaloop(); // out:  1)chi2, 2)z, 3)Etot, 4)Pskip    Zdaloop(); // out:  1)chi2, 2)z, 3)Etot, 4)Pskip
336    
   
   /*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  
   
337        
338    if ( debug ) this->Print();    if ( debug ) this->Print();
339    if ( debug ) printf(" fine evento \n");    if ( debug ) printf(" fine evento \n");
# Line 229  void CaloBragg::Process(Int_t ntr){ Line 341  void CaloBragg::Process(Int_t ntr){
341  };  };
342    
343    
344    Float_t CaloBragg::Integral(){
345      Process();
346    
347      Float_t dEpianiloop[44];
348      Int_t tz1=(Int_t)lpz;
349      Int_t ti1=(Int_t)lpisotope;
350    
351      Enetrack(&tz1, &ti1 , &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop
352    
353    
354      Float_t integ = 0.;
355      for(Int_t i=0;i<=estremi[1][0];i++){
356        //    integ += dEplan[i];
357        //printf(" step %i integ %f deplan %f \n",i,integ,dEplan[i]);
358        integ += dEpianiloop[i];
359        //    printf(" step %i integ %f deplan %f \n",i,integ,dEpianiloop[i]);
360      }
361      return integ;
362    }
363    
364    Float_t CaloBragg::LastIntegral(){
365      Process();
366    
367      Float_t integ = 0.;
368      for(Int_t i=0;i<=estremi[1][0];i++){
369        integ += dEplan[i];
370        //printf(" step %i integ %f deplan %f \n",i,integ,dEplan[i]);
371      }
372      return integ;
373    }
374    
375    
376  void CaloBragg::Draw(){  void CaloBragg::Draw(){
377    
378    Process();    Process();
379    
380   Float_t dEpianimean[44];    this->Draw(0.,0.,0.);
381   Float_t dEpianiloop[44];  
382   Float_t Depth[44];  }
383   Int_t tz=(Int_t)qtz;  
384   Int_t tz1=(Int_t)lpz;  void CaloBragg::Draw(Int_t Z, Int_t isotope, Float_t enetot){
385   Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata  
386   Enetrack(&tz1, &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop    //  Float_t dEpianimean[44];
387      Float_t dEpianiloop[44];
388      Float_t Depth[44];
389      //  Int_t tz=(Int_t)qtz;
390      Int_t tz1= Z;
391      Int_t ti1= isotope;
392      Float_t enet = enetot;
393      //  Float_t enet = lpetot;
394      
395      if ( Z > 0. && enetot > 0. ){
396        estremi[0][0] = 0;
397        estremi[1][0] = 43;
398    
399    
400   Float_t sp= spessore[0]*spessore[1];      Float_t    ytgx = 0.;
401   for(Int_t i=0;i<44;i++)Depth[i]=i*sp;      Float_t    ytgy = 0.;
402        
403        //lunghezza effettiva di silicio attraversata (mm)
404        Float_t    SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy));
405        
406        spessore[0] = (SiCross/10.) * rhoSi; //spessore silicio in g/cm2
407        
408        /*tungsteno*/
409        
410        //rapporto tra rilasci energetici nei due materiali
411        Float_t WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad
412        //gcm2W = WCross/10. * rhoW;
413        
414        //       (g/cm2W)/(g/cm2Si)
415        spessore[3] =  (WCross/10.) * rhoW;
416        Float_t a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664);  //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664);
417        spessore[1] =  a;
418        //riscala mip allo spessore attraversato
419        spessore[2] = MIP*(SiCross/ySi);  
420    
421      } else {
422        tz1=(Int_t)lpz;
423        ti1=(Int_t)lpisotope;
424        enet = lpetot;
425        //  Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata
426      
427      }
428      Enetrack(&tz1, &ti1, &enet, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop
429    
430      Float_t sp= spessore[0]*spessore[1];
431      for(Int_t i=0;i<44;i++)Depth[i]=i*sp;
432    //    //
433    gStyle->SetLabelSize(0.04);    gStyle->SetLabelSize(0.04);
434    gStyle->SetNdivisions(510,"XY");    gStyle->SetNdivisions(510,"XY");
435    //    //
436      TString hid = Form("cCaloBragg");    
437      TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
438      if ( tc ){
439        //   tc->Clear();
440      } else {
441        tc = new TCanvas(hid,hid);
442        //   tc->Divide(1,2);
443      };
444      //
445      //    TString thid = Form("hCaloBragg");      
446      //         TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
447      //    if ( th ) th->Delete();
448      //     th->Clear();
449      //     th->Reset();
450      //    } else {
451      //    th = new TH2F(thid,thid,300,-0.5,300.,1000,0.,150.);
452      //    th->SetMarkerStyle(20);
453      //    };
454      //
455      tc->cd();
456      TString thid2 = Form("hCaloBragg2");  
457      TH2F *th2  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid2));
458      if ( th2 ) th2->Delete();
459      th2 = new TH2F(thid2,thid2,300,-0.5,300.,1000,0.,150.);
460      th2->SetMarkerStyle(20);
461      th2->SetMarkerColor(kRed);
462      //
463      TString thid3 = Form("hCaloBragg3");  
464      TH2F *th3  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid3));
465      if ( th3 ) th3->Delete();
466      th3 = new TH2F(thid3,thid3,300,-0.5,300.,1000,0.,150.);
467      th3->SetMarkerStyle(20);
468      th3->SetMarkerColor(kBlue);
469    
470    
471      //  tc->cd(1);
472      //
473      //    for(Int_t i=0;i<=estremi[1][0];i++)th->Fill(Depth[i],dEpianimean[i]);
474      for(Int_t i=0;i<=estremi[1][0];i++)th2->Fill(Depth[i],calorimetro[i][1]*MIP);
475      //    th->Draw();
476      th2->Draw("same");
477    
478          TString hid = Form("cCaloBragg");          //  tc->cd(2);
479          TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));    tc->cd();
480      //
481          if ( tc ){    for(Int_t i=0;i<=estremi[1][0];i++){
482  //       tc->Clear();      th3->Fill(Depth[i],dEpianiloop[i]);
483          } else {      //    printf(" i %i Depth %f depianiloop %f \n",i,Depth[i],dEpianiloop[i]);
484           tc = new TCanvas(hid,hid);    }
485           tc->Divide(1,2);    th3->Draw();
486          };    th2->Draw("same");
         //  
         TString thid = Form("hCaloBragg");        
         TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));  
         if ( th ) th->Delete();  
 //       th->Clear();  
 //       th->Reset();  
 //      } else {  
         th = new TH2F(thid,thid,300,-0.5,300.,1000,0.,150.);  
         th->SetMarkerStyle(20);  
 //      };  
         //  
         TString thid2 = Form("hCaloBragg2");      
         TH2F *th2  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid2));  
         if ( th2 ) th2->Delete();  
         th2 = new TH2F(thid2,thid2,300,-0.5,300.,1000,0.,150.);  
         th2->SetMarkerStyle(20);  
         th2->SetMarkerColor(kRed);  
         //  
         TString thid3 = Form("hCaloBragg3");      
         TH2F *th3  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid3));  
         if ( th3 ) th3->Delete();  
         th3 = new TH2F(thid3,thid3,300,-0.5,300.,1000,0.,150.);  
         th3->SetMarkerStyle(20);  
         th3->SetMarkerColor(kBlue);  
   
   
         tc->cd(1);  
         //  
         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++)th2->Fill(Depth[i],calorimetro[i][1]*MIP);  
         th->Draw();  
         th2->Draw("same");  
   
         tc->cd(2);  
         //  
         for(Int_t i=0;i<=estremi[1][0];i++)th3->Fill(Depth[i],dEpianiloop[i]);  
         th3->Draw();  
         th2->Draw("same");  
487    
488          tc->Modified();    tc->Modified();
489          tc->Update();    tc->Update();
490    
491    //    //
492    gStyle->SetLabelSize(0);    gStyle->SetLabelSize(0);
# Line 309  void CaloBragg::Draw(){ Line 498  void CaloBragg::Draw(){
498    
499  void CaloBragg::LoadParam(){  void CaloBragg::LoadParam(){
500    
501    //    // elem[Z-1][isotop] 0 is the most common one
502    elem[0] = 1.00794; //H  1    //
503    elem[1] = 4.0026;  //He 2  
504    elem[2] = 6.941;   //Li 3    elem[0][0] = 1.00782;   //H  1
505    elem[3] = 9.012182;//Be 4    elem[0][1] = 2.01410;   // 2H (Isotope)
506    elem[4] = 10.811;  //B  5    elem[0][2] = -1.;
507    elem[5] = 12.0107; //C  6    elem[0][3] = -1.;
508    elem[6] = 14.00674;//N  7    elem[0][4] = -1.;
509    elem[7] = 15.9994; //O  8    elem[0][5] = -1.;
510    elem[8] = 18.9984; //F  9    elem[0][6] = -1.;
511    elem[9] = 20.1797; //Ne 10  
512    elem[10] = 22.98977;//Na 11    elem[1][0] = 4.002602;  //He 2
513    elem[11] = 24.3050; //Mg 12    elem[1][1] = 3.016029;  // 3He (Isotope)
514    elem[12] = 26.9815; //Al 13    elem[1][2] = -1.;
515    elem[13] = 28.0855; //Si 14    elem[1][3] = -1.;
516    elem[14] = 30.974;  //P  15    elem[1][4] = -1.;
517    elem[15] = 32.066;  //S  16    elem[1][5] = -1.;
518    elem[16] = 35.4527; //Cl 17    elem[1][6] = -1.;
519    elem[17] = 39.948;  //Ar 18  
520    elem[18] = 39.0983; //K  19    elem[2][0] = 7.016004;  //Li 3
521    elem[19] = 40.078;  //Ca 20    elem[2][1] = 6.015123;  //6Li (Isotope)
522    elem[20] = 44.95591;//Sc 21    elem[2][2] = -1.;
523    elem[21] = 47.867;  //Ti 22    elem[2][3] = -1.;
524    elem[22] = 50.9415; //V  23    elem[2][4] = -1.;
525    elem[23] = 51.9961; //Cr 24    elem[2][5] = -1.;
526    elem[24] = 54.938049;//Mn 25    elem[2][6] = -1.;
527    elem[25] = 55.845;  //Fe 26  
528    elem[26] = 58.9332; //Co 27    elem[3][0] = 9.012182;  //Be 4
529    elem[27] = 58.6934; //Ni 28    elem[3][1] = 10.01353;  //10Be (Isotope) (most stable)
530    elem[28] = 63.546; //Cu 29    elem[3][2] = -1.;
531    elem[29] = 65.39; //Zn 30    elem[3][3] = -1.;
532    elem[30] = 69.723; //Ga 31    elem[3][4] = -1.;
533    elem[31] = 72.61; //Ge 32    elem[3][5] = -1.;
534      elem[3][6] = -1.;
535    
536      elem[4][0] = 11.0093;   //B 5
537      elem[4][1] = 10.01294;  //10B  (Isotope)
538      elem[4][2] = -1.;
539      elem[4][3] = -1.;
540      elem[4][4] = -1.;
541      elem[4][5] = -1.;
542      elem[4][5] = -1.;
543    
544      elem[5][0] = 12.0107;   //C  6
545      elem[5][1] = 13.00335;  //13C (Isotope)
546      elem[5][2] = -1.;
547      elem[5][3] = -1.;
548      elem[5][4] = -1.;
549      elem[5][5] = -1.;
550      elem[5][5] = -1.;
551    
552      elem[6][0] = 14.00674;  //N  7
553      elem[6][1] = 15.00011;  //15N  (Isotope)
554      elem[6][2] = -1.;
555      elem[6][3] = -1.;
556      elem[6][4] = -1.;
557      elem[6][5] = -1.;
558      elem[6][5] = -1.;
559    
560      elem[7][0] = 15.99491;  //O  8
561      elem[7][1] = 17.99916;  //18O  (Isotope)
562      elem[7][2] = 16.99916;  //17O  (Isotope)
563      elem[7][3] = -1.;
564      elem[7][4] = -1.;
565      elem[7][5] = -1.;
566      elem[7][5] = -1.;
567    
568      elem[8][0] = 18.99840;  //F  9
569      elem[8][1] = -1.;
570      elem[8][2] = -1.;
571      elem[8][3] = -1.;
572      elem[8][4] = -1.;
573      elem[8][5] = -1.;
574      elem[8][5] = -1.;
575    
576      elem[9][0] = 19.99244;  //Ne 10
577      elem[9][1] = 21.99138;  //22Ne (Isotope)
578      elem[9][2] = 20.99384;  //21Ne 10
579      elem[9][3] = -1.;
580      elem[9][4] = -1.;
581      elem[9][5] = -1.;
582      elem[9][6] = -1.;
583    
584      elem[10][0] = 22.98977; //Na 11
585      elem[10][1] = 21.99444; //22Na (Isotope) (most stable)
586      elem[10][2] = -1.;
587      elem[10][3] = -1.;
588      elem[10][4] = -1.;
589      elem[10][5] = -1.;
590      elem[10][6] = -1.;
591    
592      elem[11][0] = 23.98504; //Mg 12
593      elem[11][1] = 25.98259; //26Mg (Isotope)
594      elem[11][2] = 24.98504; //25Mg (Isotope)
595      elem[11][3] = -1.;
596      elem[11][4] = -1.;
597      elem[11][5] = -1.;
598      elem[11][6] = -1.;
599    
600      elem[12][0] = 26.98154; //Al 13
601      elem[12][1] = 25.98489; //26Al (Isotope) (most stable)
602      elem[12][2] = -1.;
603      elem[12][3] = -1.;
604      elem[12][4] = -1.;
605      elem[12][5] = -1.;
606      elem[12][6] = -1.;
607    
608      elem[13][0] = 27.97692; //Si 14
609      elem[13][1] = 28.97649; //29Si (Isotope)
610      elem[13][2] = 29.97377; //30Si (Isotope)
611      elem[13][3] = -1.;
612      elem[13][4] = -1.;
613      elem[13][5] = -1.;
614      elem[13][6] = -1.;
615    
616      elem[14][0] = 30.97376;  //P  15
617      elem[14][1] = -1.;
618      elem[14][2] = -1.;
619      elem[14][3] = -1.;
620      elem[14][4] = -1.;
621      elem[14][5] = -1.;
622      elem[14][6] = -1.;
623    
624      elem[15][0] = 31.97207;  //S  16
625      elem[15][1] = 33.96787;  //34S  (Isotope)
626      elem[15][2] = 32.97146;  //33S  (Isotope)
627      elem[15][3] = 35.96708;  //36S  (Isotope)
628      elem[15][4] = -1.;
629      elem[15][5] = -1.;
630      elem[15][6] = -1.;
631    
632      elem[16][0] = 34.96885; //Cl 17
633      elem[16][1] = 36.96831; //37Cl 17
634      elem[16][2] = 35.96890; //36Cl (Isotope)
635      elem[16][3] = -1.;
636      elem[16][4] = -1.;
637      elem[16][5] = -1.;
638      elem[16][6] = -1.;
639    
640      elem[17][0] = 39.962383;  //Ar 18
641      elem[17][1] = 35.967545;  //36Ar (Isotope)
642      elem[17][2] = 37.962732;  //38Ar (Isotope)
643      elem[17][3] = 38.964313;  //39Ar (Isotope)
644      elem[17][4] = -1.;
645      elem[17][5] = -1.;
646      elem[17][6] = -1.;
647    
648      elem[18][0] = 38.963707; //K  19
649      elem[18][1] = 40.961825; //41K  (Isotope)
650      elem[18][2] = 39.963998; //40K  (Isotope)
651      elem[18][3] = -1.;
652      elem[18][4] = -1.;
653      elem[18][5] = -1.;
654      elem[18][6] = -1.;
655    
656      elem[19][0] = 39.962590;  //Ca 20
657      elem[19][1] = 43.955482;  //44Ca (Isotope)
658      elem[19][2] = 41.958618;  //42Ca (Isotope)
659      elem[19][3] = 42.958767;  //43Ca (Isotope)
660      elem[19][4] = 45.953693;  //46Ca (Isotope)
661      elem[19][5] = 40.962278;  //41Ca (Isotope)
662      elem[19][6] = -1.;
663    
664      elem[20][0] = 44.955912;//Sc 21
665      elem[20][1] = -1.;
666      elem[20][2] = -1.;
667      elem[20][3] = -1.;
668      elem[20][4] = -1.;
669      elem[20][5] = -1.;
670      elem[20][6] = -1.;
671    
672      elem[21][0] = 47.947946;  //Ti 22
673      elem[21][1] = 45.952632;  //46Ti (Isotope)
674      elem[21][2] = 46.951763;  //47Ti (Isotope)
675      elem[21][3] = 48.947870;  //49Ti (Isotope)
676      elem[21][4] = 49.944791;  //50Ti (Isotope)
677      elem[21][5] = 43.959690;  //44Ti (Isotope) (half life 60y)
678      elem[21][6] = -1.;
679    
680      elem[22][0] = 50.943960; //V  23
681      elem[22][1] = 49.947158; //50V  (Isotope)
682      elem[22][2] = -1.;
683      elem[22][3] = -1.;
684      elem[22][4] = -1.;
685      elem[22][5] = -1.;
686      elem[22][6] = -1.;
687    
688      elem[23][0] = 51.940507; //Cr 24
689      elem[23][1] = 52.940649; //53Cr (Isotope)
690      elem[23][2] = 49.946044; //50Cr (Isotope)
691      elem[23][3] = 53.938880; //54Cr (Isotope)
692      elem[23][4] = -1.;
693      elem[23][5] = -1.;
694      elem[23][6] = -1.;
695    
696      elem[24][0] = 54.938049;//Mn 25
697      elem[24][1] = 52.941290;//53Mn (Isotope)
698      elem[24][2] = -1.;
699      elem[24][3] = -1.;
700      elem[24][4] = -1.;
701      elem[24][5] = -1.;
702      elem[24][6] = -1.;
703    
704      elem[25][0] = 55.934937;  //Fe 26
705      elem[25][1] = 53.939610;  //54Fe (Isotope)
706      elem[25][2] = 56.935394;  //57Fe (Isotope)
707      elem[25][3] = 57.933276;  //58Fe (Isotope)
708      elem[25][4] = 59.934072;  //58Fe (Isotope)
709    
710      elem[26][0] = 58.933195; //Co 27
711      elem[26][1] = 59.933817; //60Co (Isotope)
712      elem[26][2] = -1.;
713      elem[26][3] = -1.;
714      elem[26][4] = -1.;
715      elem[26][5] = -1.;
716      elem[26][6] = -1.;
717    
718    
719      elem[27][0] = 57.935343; //Ni 28
720      elem[27][1] = 61.928345; //62Ni (Isotope)
721      elem[27][2] = 59.930786; //60Ni (Isotope)
722      elem[27][3] = 60.931056; //61Ni (Isotope)
723      elem[27][4] = 63.927966; //64Ni (Isotope)
724      elem[27][5] = 58.934346; //59Ni (Isotope)
725      elem[27][6] = -1.;
726    
727      elem[28][0] = 62.929597; //Cu 29
728      elem[28][1] = 64.927789; //65Cu (Isotope)
729      elem[28][2] = -1.;
730      elem[28][3] = -1.;
731      elem[28][4] = -1.;
732      elem[28][5] = -1.;
733      elem[28][6] = -1.;
734    
735      elem[29][0] = 63.929142; //Zn 30
736      elem[29][1] = 65.926033; //66Zn (Isotope)
737      elem[29][2] = 67.924844; //68Zn (Isotope)
738      elem[29][3] = 66.927127; //67Zn (Isotope)
739      elem[29][4] = 69.925319; //70Zn (Isotope)
740      elem[29][5] = -1.;
741      elem[29][6] = -1.;
742    
743      elem[30][0] = 68.925573; //Ga 31
744      elem[30][1] = 70.924701; //71Ga (Isotope)
745      elem[30][2] = -1.;
746      elem[30][3] = -1.;
747      elem[30][4] = -1.;
748      elem[30][5] = -1.;
749      elem[30][6] = -1.;
750    
751      elem[31][0] = 73.921177; //Ge 32
752      elem[31][1] = 71.922075; //72Ge (Isotope)
753      elem[31][2] = 69.924247; //70Ge (Isotope)
754      elem[31][3] = 75.921403; //76Ge (Isotope)
755      elem[31][4] = 73.923459; //73Ge (Isotope)
756      elem[31][5] = -1.;
757      elem[31][6] = -1.;
758    
759    
760  //parametri calorimetro    //parametri calorimetro
761    NPLA = 22;    NPLA = 22;
762    NCHA = 96;    NCHA = 96;
763    nView = 2;    nView = 2;
# Line 365  void CaloBragg::LoadParam(){ Line 778  void CaloBragg::LoadParam(){
778    pigr = 3.1415;    pigr = 3.1415;
779    Na = 6.02e-23;    Na = 6.02e-23;
780    ZA = 0.49; /*Z/A per Si*/    ZA = 0.49; /*Z/A per Si*/
781    ISi =182e-06; /*MeV*/    //  ISi =182e-06; /*MeV*/
782      ISi = 171e-06; /*MeV*/
783      IW  = 735e-06; /*MeV*/
784      //  ISi =0.0001059994; /*GeV!!*/ no era giusto!!
785    Me = 0.511; /* MeV*/    Me = 0.511; /* MeV*/
786    MassP = 931.27;/*MeV*/    MassP = 931.27;/*MeV*/
787    r2 = 7.95e-26; /*ro*ro in cm */    r2 = 7.95e-26; /*ro*ro in cm */
# Line 393  void CaloBragg::conversione(){ Line 809  void CaloBragg::conversione(){
809    Float_t ytgy = 0;    Float_t ytgy = 0;
810    Float_t a = 0.;    Float_t a = 0.;
811    
   
812    /*silicio*/    /*silicio*/
813    ytgx = ySi * L2->GetCaloLevel2()->tanx[0];    ytgx = ySi * L2->GetCaloLevel2()->tanx[0];
814    ytgy = ySi * L2->GetCaloLevel2()->tany[0];    ytgy = ySi * L2->GetCaloLevel2()->tany[0];
# Line 401  void CaloBragg::conversione(){ Line 816  void CaloBragg::conversione(){
816    //lunghezza effettiva di silicio attraversata (mm)    //lunghezza effettiva di silicio attraversata (mm)
817    SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy));    SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy));
818        
819    spessore[0] = SiCross/10. * rhoSi; //spessore silicio in g/cm2    spessore[0] = (SiCross/10.) * rhoSi; //spessore silicio in g/cm2
     
820        
821    /*tungsteno*/    /*tungsteno*/
822    ytgx = yW * L2->GetCaloLevel2()->tanx[0];    ytgx = yW * L2->GetCaloLevel2()->tanx[0];
# Line 412  void CaloBragg::conversione(){ Line 826  void CaloBragg::conversione(){
826    WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad    WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad
827    //gcm2W = WCross/10. * rhoW;    //gcm2W = WCross/10. * rhoW;
828    
  a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664);  //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664);  
   
829    //       (g/cm2W)/(g/cm2Si)    //       (g/cm2W)/(g/cm2Si)
830      spessore[3] =  (WCross/10.) * rhoW;
831      a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664);  //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664);
832    spessore[1] =  a;    spessore[1] =  a;
   
     
833    //riscala mip allo spessore attraversato    //riscala mip allo spessore attraversato
834    spessore[2] = MIP*(SiCross/ySi);      spessore[2] = MIP*(SiCross/ySi);  
   
835  };//end conversione  };//end conversione
836    
837    
838    
839    
840    
841  void CaloBragg::BetheBloch(Float_t *x, Float_t *z, Float_t *Mass, Float_t *gam, Float_t *Bet, Float_t *out){  void CaloBragg::BetheBloch(Float_t *x, Float_t *z, Float_t *Mass, Float_t *gam, Float_t *Bet, Float_t *out, Float_t II){
842    
843    //rilascio energetico con bethe bloch con correzioni    //rilascio energetico con bethe bloch con correzioni
844    //in:    x: g/cm2    //in:    x: g/cm2
# Line 445  void CaloBragg::BetheBloch(Float_t *x, F Line 856  void CaloBragg::BetheBloch(Float_t *x, F
856    Float_t lg =0.;    Float_t lg =0.;
857    Float_t Energia=0.;    Float_t Energia=0.;
858    Float_t C=0.;    Float_t C=0.;
859      Float_t INo = ISi;
860    
861      if ( usenewBB ) INo = II;
862    
863    eta = (*gam)*(*Bet);    eta = (*gam)*(*Bet);
864    
865    //Bet=3/gam;  SQ(*gam) * SQ(*Bet)    //Bet=3/gam;  SQ(*gam) * SQ(*Bet)
866    Wmax = 2.* Me * SQ(eta) / (1. + 2.*(*gam)*Me/(*Mass) + SQ(Me)/SQ(*Mass));    Wmax = 2.* Me * SQ(eta) / (1. + 2.*(*gam)*Me/(*Mass) + SQ(Me)/SQ(*Mass));
867        
868    lg = 2.* Me * SQ(eta) * Wmax / SQ(ISi);    lg = 2.* Me * SQ(eta) * Wmax / SQ(INo);
869                                                                                                 //  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;
870    C=(0.42237*pow(eta,-2.) + 0.0304*pow(eta,-4.) - 0.00038*pow(eta,-6.))*pow(10.,-6.)* pow(ISi,2.) +    C=(0.42237*pow(eta,-2.) + 0.0304*pow(eta,-4.) - 0.00038*pow(eta,-6.))*pow(10.,-6.)* pow(INo,2.) +
871      (3.858*pow(eta,-2.) - 0.1668*pow(eta,-4.) + 0.00158*pow(eta,-6.))*pow(10.,-9.)*pow(ISi,3.);      (3.858*pow(eta,-2.) - 0.1668*pow(eta,-4.) + 0.00158*pow(eta,-6.))*pow(10.,-9.)*pow(INo,3.);
872        
873    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.);  
     
874    
   /*se ho gia' lo spessore in g/cm2 non se mejo???*/  
875    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.);
876    
877    *out =Energia;//out    *out =Energia;//out
# Line 471  void CaloBragg::BetheBloch(Float_t *x, F Line 881  void CaloBragg::BetheBloch(Float_t *x, F
881    
882    
883    
884  void CaloBragg::ELOSS(Float_t *dx, Int_t *Z, Float_t *Etot, Float_t *out){  void CaloBragg::ELOSS(Float_t *dx, Int_t *Z, Int_t *isotope, Float_t *Etot, Float_t *out, Float_t II){
885        
886    /*perdita di energia per ioni pesanti (come da routine geant)*/    /*perdita di energia per ioni pesanti (come da routine geant)*/
887    //  in : dx    => spessore g/cm2    //  in : dx    => spessore g/cm2
# Line 482  void CaloBragg::ELOSS(Float_t *dx, Int_t Line 892  void CaloBragg::ELOSS(Float_t *dx, Int_t
892    
893    
894    Float_t Q=0.;    Float_t Q=0.;
895    Float_t v=0.;    Float_t v=0.;                                                        
                                                         //  Float_t Mass=0.;  
896    Float_t gam=0.;    Float_t gam=0.;
897    Float_t Bet=0.;    Float_t Bet=0.;
898    Float_t dEP=0.;    Float_t dEP=0.;
899        
900    // gamma  //  Mass = A * MassP; /*in Mev/c2*/    // gamma  //  Mass = A * MassP; /*in Mev/c2*/
901    gam =  (*Etot)/(elem[*Z-1]*MassP); // E = gamma M c2    gam =  (*Etot)/(elem[*Z-1][*isotope]*MassP); // E = gamma M c2
902        
903    
904    Bet = sqrt((SQ(gam) -1.)/SQ(gam));    Bet = sqrt((SQ(gam) -1.)/SQ(gam));
905        
906    v= 121.4139*(Bet/pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet/pow((*Z),(2./3.))));    //  v= 121.4139*(Bet/pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet/pow((*Z),(2./3.))));
907      v= 121.4139*(Bet*pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet*pow((*Z),(2./3.)))); // EMI AAAAGGH!!
908    
909    //carica effettiva    //carica effettiva
910    Q= (*Z)*(1- (1.034 - 0.1777*exp(-0.08114*(*Z)))*exp(-v));    Q= (*Z)*(1- (1.034 - 0.1777*exp(-0.08114*(*Z)))*exp(-v));
911    
912    //perdita energia per un protone    //perdita energia per un protone
913    Float_t protone =1.;    Float_t protone =1.;
914    Float_t Mass=MassP;    //  Float_t Mass=(elem[*Z-1]*MassP); //EMI
915    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);
916    
917      BetheBloch(dx, &protone, &MassP, &gam, &Bet, &dEP, II);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP); //EMI
918        
919    *out= (SQ(Q)*(dEP));//*dx;    *out= (SQ(Q)*(dEP));//*dx;
920    
# Line 512  void CaloBragg::ELOSS(Float_t *dx, Int_t Line 924  void CaloBragg::ELOSS(Float_t *dx, Int_t
924    
925    
926    
927  void CaloBragg::Enetrack(Int_t* Z, Float_t* E0, Float_t* primo,Float_t* ultimo, Float_t out[]){  void CaloBragg::Enetrack(Int_t* Z, Int_t* isotope, Float_t* E0, Float_t* primo,Float_t* ultimo, Float_t out[]){
928    
929    //calcola energia rilasciata sulla traccia (usa ELOSS)    //calcola energia rilasciata sulla traccia (usa ELOSS)
930    //  in : Z             =>carica    //  in : Z             =>carica
# Line 530  void CaloBragg::Enetrack(Int_t* Z, Float Line 942  void CaloBragg::Enetrack(Int_t* Z, Float
942    //azzero energia rilasciata sui piani    //azzero energia rilasciata sui piani
943    memset(out, 0, 2*NPLA*sizeof(Float_t));    memset(out, 0, 2*NPLA*sizeof(Float_t));
944        
945    Float_t Massa = (elem[(*Z)-1] * MassP);    Float_t Massa = (elem[(*Z)-1][*isotope] * MassP);
   
946    
947    //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++){  
  for( Int_t ipla=((int)(*primo)); ipla<= ((int)(*ultimo)); ipla++){  
948      dE=0.;      dE=0.;
   
949      //spessore silicio corretto x inclinazione, z, energia, out:rilascio      //spessore silicio corretto x inclinazione, z, energia, out:rilascio
950      ELOSS(&spessore[0], Z, &Ezero, &dE);//spessore in g/cm2!!      ELOSS(&spessore[0], Z , isotope , &Ezero, &dE, ISi);//spessore in g/cm2!!
951      if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop  
952            if(dE!=dE) return; //controlla che non sia un NaN
953        
954            if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop
955        out[ipla] = Ezero - Massa; //MeV        out[ipla] = Ezero - Massa; //MeV
956        return;        return;
957                
958      }else{      }else{
959        out[ipla] = dE; //MeV        out[ipla] = dE; //MeV
960        Ezero = Ezero - dE;//energia residua        Ezero = Ezero - dE;//energia residua
961          if ( debug ) printf(" zompa %i out %f dE %f ezero %f \n",ipla,out[ipla],dE,Ezero);
962      };      };
   
963      //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
964      if(ipla%2 == 0){      if(ipla%2 == 0){
965        /*tungsteno*/        /*tungsteno*/
966        dE=0.;        dE=0.;
967        Float_t sp= spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato  in g/cm2        Float_t sp = 0.;
968        ELOSS(&sp, Z, &Ezero, &dE);        Float_t II = ISi;
969        //      cout<<"perdita per piano di W ="<<dE<<endl;        if ( usenewBB ){
970            sp = spessore[3];        
971            II = IW;
972          } else {
973            sp = spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato  in g/cm2
974          }
975          //      printf(" sp %f II %f \n",sp,II);
976          ELOSS(&sp, Z, isotope , &Ezero, &dE,II);
977        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
978          return;          return;
979        }else{        }else{
980          Ezero = Ezero -dE;//energia residua          Ezero = Ezero -dE;//energia residua
         //      cout<<"w calc "<<dE<<endl;  
981        };        };
982      };          };    
   };//fine loop piani  
983    
984  //   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  
985    
986    
987    };//end Enetrack
988    
989    
990    
# Line 588  void CaloBragg::chiquadro(Float_t dE[], Line 1003  void CaloBragg::chiquadro(Float_t dE[],
1003    Float_t badplane=0.;    Float_t badplane=0.;
1004    Float_t badplanetot=0.;    Float_t badplanetot=0.;
1005    Float_t w,wi;    Float_t w,wi;
1006      //
1007    for(Int_t ipla=0; ipla<2*NPLA; ipla++){    if ( newchi2 ){
1008          ndf = 0;
1009      //tutti i piani attraversati dalla traiettoria      sum = 0.;
1010      if(calorimetro[ipla][0] != -1.){      for( Int_t ipla=((int)(estremi[0][0])); ipla<= ((int)(estremi[1][0])); ipla++){
1011        w=0.; //normalizzazione;        sum += pow((dE[ipla] - (calorimetro[ipla][1] * spessore[2]))/(0.05*dE[ipla]),2.);
1012        wi=1.;//peso        //      printf(" quiqui: dE %f calor %f spessore[2] %f \n",dE[ipla],spessore[2]*calorimetro[ipla][1],spessore[2]);
1013          ndf++;
1014        }
1015        ndf -= 2;
1016        if ( ndf > 0 ) sum /= (float)ndf;
1017        out[0] = sum;
1018        out[1] = 0.;
1019        out[2] = (int)(estremi[1][0])-ndf;
1020        //    printf(" sum %f ndf %i \n ",sum,ndf);
1021      } else {
1022        for(Int_t ipla=0; ipla<2*NPLA; ipla++){
1023          //tutti i piani attraversati dalla traiettoria
1024          if(calorimetro[ipla][0] != -1.){ //
1025            w=0.; //normalizzazione;
1026            wi=1.;//peso
1027    
1028        //tolgo piani attraversati dalla traccia ma precedenti il piano individuato come ingresso          //tolgo piani attraversati dalla traccia ma precedenti il piano individuato come ingresso
1029        if (ipla<estremi[0][0])  wi=0.;          if (ipla<estremi[0][0])  wi=0.;
1030                
1031        //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
1032        //if((ipla>estremi[1][0]) && (calorimetro[ipla][1] >0.) ) wi=0.;          //if((ipla>estremi[1][0]) && (calorimetro[ipla][1] >0.) ) wi=0.;
1033        if((ipla>estremi[1][0])) wi=0.;          if((ipla>estremi[1][0])) wi=0.;
1034    
1035        //normalizzazione          //normalizzazione      
1036        if (calorimetro[ipla][1] != 0.)  w=1./(calorimetro[ipla][1]* MIP);                  if (calorimetro[ipla][1] != 0.)  w=1./(calorimetro[ipla][1]* MIP);    //        
1037                
1038        //tolgo piani con rilasci inferiori al 30% del precedente          //tolgo piani con rilasci inferiori al 30% del precedente
1039        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
1040          wi=0.;            wi=0.;
1041          //se sono piani intermedi (non si è fermta) li considero non buoni            //se sono piani intermedi (non si e' fermta) li considero non buoni
1042          if( (ipla <= estremi[1][0]) && (calorimetro[ipla][1] !=0.)){            if( (ipla <= estremi[1][0]) && (calorimetro[ipla][1] !=0.)){//
1043            badplane+=1.;              badplane+=1.;
1044            badplanetot+=1.;              badplanetot+=1.;
1045          };            };
1046        };          };
1047        
1048        //meno peso ai piani con rilasci maggiori di 1000 MIP          //meno peso ai piani con rilasci maggiori di 1000 MIP
1049          if(calorimetro[ipla][1] > 1000) wi=0.5;          //      if(calorimetro[ipla][1] > 1000) wi=0.5;
1050                        if(calorimetro[ipla][1] > 1200.) wi=0.5;
1051        Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));             if(debug) printf("chiquadro start \n ");
1052            Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));
1053        sum += SQ(arg); // w*wi*(dEpiani[p][v]-(eplane[p][v]*MIP))));//( dEpiani[p][v] - (eplane[p][v]*MIP));  
1054            sum += SQ(arg); // w*wi*(dEpiani[p][v]-(eplane[p][v]*MIP))));//( dEpiani[p][v] - (eplane[p][v]*MIP));
1055        //se trovo piano non buono (tolto quindi wi=0) non modifico il piano precedente          if(debug){
1056        if(wi != 0.){            printf("dedx  calcolata  %f e reale  %f  \n",dE[ipla],(calorimetro[ipla][1] * MIP));
1057          PianoPrecedente= calorimetro[ipla][1];//tengo piano precedente          }
1058          badplane = 0.;//azzero contatore piani scartati consecutivi          //se trovo piano non buono (tolto quindi wi=0) non modifico il piano precedente
1059            if(wi != 0.){//
1060              PianoPrecedente= calorimetro[ipla][1];//tengo piano precedente
1061              badplane = 0.;//azzero contatore piani scartati consecutivi
1062            };
1063        };        };
     };  
             
     if(badplane > 2) out[1] =79.;  
1064                        
1065    };//fine loop piani        //da Emi
1066    //chi2,frammentato,pskip          if(badplane > 2){      
1067    out[0]=sum;          //      printf(" AAAAAAAAAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG\n");
1068    out[2]=badplanetot;          out[1] =79.;
1069            break;
1070    if(out[1] ==79.)cout<<"frammentato !!!!!"<<endl;        };
1071        
1072        };//fine loop piani
1073        //chi2,frammentato,pskip  
1074        out[0]=sum;
1075        out[2]=badplanetot;
1076      }  
1077  };//end chiquadro  };//end chiquadro
1078    
1079    
1080    
1081  //loopze(Integrale,&zero,&zmin, &zmax, spess1 , estr3, calo1, bestchi2);  void CaloBragg::loopze( Float_t step, Float_t E0,Float_t Zstart, Float_t Zlimite, Int_t nostep = 1000){
1082  void CaloBragg::loopze( Float_t step, Float_t E0,Float_t Zstart, Float_t Zlimite){    //
 //  
1083    //loop su z ed energie per trovare miglior z (ed energia)    //loop su z ed energie per trovare miglior z (ed energia)
1084    //in:  nloop     => energia massima da provare (nloop x E0)    //in:  nloop     => energia massima da provare (nloop x E0)
1085    //     E0        => energia iniziale (intergale)            //     E0        => energia iniziale (intergale)        
# Line 655  void CaloBragg::loopze( Float_t step, Fl Line 1089  void CaloBragg::loopze( Float_t step, Fl
1089    //out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2    //out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2
1090    //    //
1091    
1092     //printf("entrato");
   
   Float_t dEplan[2*NPLA];//energia rilasciata calcolata  
1093    memset(dEplan,0,2*NPLA*sizeof(Float_t));    memset(dEplan,0,2*NPLA*sizeof(Float_t));
1094    
1095    Int_t Z = 0;// z iniziale    Int_t Z = 0;// z iniziale
1096      
1097      Int_t isotope=0;
1098    
1099    Float_t Massa = 0.;    Float_t Massa = 0.;
1100    
1101    Float_t Stepint =(step)/100.;//passo per il calcolo di energia //era 1000    Float_t Stepint =(step)/(Float_t)nostep;//passo per il calcolo di energia
1102    
1103    Float_t energia =0.;//energia del loop    Float_t energia =0.;//energia del loop
1104    
1105    Float_t chi2[3] = {0,0,0};//out dal calcolo chi2: chi2, piani consecutivi saltati, piani totali saltati    Float_t chi2[3] = {0,0,0};//out dal calcolo chi2: chi2, piani consecutivi saltati, piani totali saltati
1106      
1107      Int_t zmin = (int)Zstart;  
1108    Int_t max=32;//max z di cui so la massa :P    Int_t max=32;//max z di cui so la massa :P
1109    if((Zlimite)<=31) max=(int)(Zlimite) + 1;    if((Zlimite)<=31) max=(int)(Zlimite) + 1;
1110    
1111    if(debug)  printf("loopze inizio max  %d \n",max);
1112      if ( fzeta > 0. ){
1113        zmin = fzeta;
1114        max = fzeta+1;
1115      }
1116        
1117    Int_t colmax=32;    Int_t colmax=32;
1118    Int_t rowmax=3000;    Int_t rowmax=3000;
1119      Int_t isomax=7;
1120    
1121    Float_t matrixchi2[colmax][rowmax][3];    Float_t matrixchi2[colmax][isomax][rowmax][3];
1122    memset(matrixchi2, 0, colmax*rowmax*3*sizeof(Float_t));    memset(matrixchi2, 0, colmax*isomax*rowmax*3*sizeof(Float_t));
1123    
1124      Int_t imin = 1-nostep/2;
1125      Int_t imax = nostep/2;
1126    
1127    //loop elementi      //loop elementi  
1128    for(Int_t inucl=(int)(Zstart); inucl<max; inucl++){    for(Int_t inucl=zmin; inucl<max; inucl++){
1129            
1130      Z= inucl;      Z= inucl;
1131      
1132      Massa = elem[inucl-1]*MassP;    //loop isotopi
1133          while ( elem[inucl-1][isotope] > 0. ){  
1134        Massa = elem[inucl-1][isotope]*MassP;
1135    
1136      //loop energia      //loop energia
1137      for(Int_t iene= 0; iene<1000; iene++){      Int_t iene2 = -1;
       
       energia=  Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica..  
1138    
1139        Enetrack(&Z, &energia, &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani      //    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 <=??
1140          
1141        //calcolo chi2      for(Int_t iene= imin; iene<imax; iene++){// da non cambiare in base a Stepint altrimenti cambia la matrice bestchi2!!!cosi' non raggiungo mai integrale!!!!! mettere <=??
1142        chiquadro(dEplan,chi2);           iene2++;
1143                     energia=  Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica..
1144            
1145            
1146                     if( fene > 0. ) energia=fene;   //forza l'energia
1147                     if (debug) printf("loopze energia %f, z %d, isotopo %d ,iene %d\n",energia,Z,isotope,iene);
1148    //              printf(" energia %f , forzata %f \n",energia,fene);
1149                     Enetrack(&Z, &isotope, &energia, &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani
1150                
1151        if( (chi2[1] != 79) ){//salto quelli che frammentano           chiquadro(dEplan,chi2); //calcolo chi2
1152          matrixchi2[inucl][iene][0]=chi2[0];//valore chi2 per questo z a questa energia                   if (debug) printf("loopze chi %f \n",chi2[0]);
1153          matrixchi2[inucl][iene][1]=energia;//energia per questo chi2       if(debug && TMath::Finite(chi2[0])==1 && (TMath::IsNaN(chi2[0])!=1) )    printf("loopze    fin mat %f \n",chi2[0]);
1154          matrixchi2[inucl][iene][2]=chi2[2];//piani saltati nel chi2        //      printf(" last deplan from: Z = %i iene %i energia %f chi2 %f \n",inucl,iene,energia,chi2[0]);
1155        } else {             if( (chi2[1] != 79.) ){//salto quelli che frammentano
1156          matrixchi2[inucl][iene][0]=1000;//valore chi2 per questo z a questa energia                  matrixchi2[inucl][isotope][iene2][0]=chi2[0];//valore chi2 per questo z a questa energia
1157          matrixchi2[inucl][iene][1]=1000;//energia per questo chi2                  matrixchi2[inucl][isotope][iene2][1]=energia;//energia per questo chi2
1158          matrixchi2[inucl][iene][2]=1000;//piani saltati nel chi2                          matrixchi2[inucl][isotope][iene2][2]=chi2[2];//piani saltati nel chi2
1159                if( fene > 0. ) break;
1160             } else {      
1161                            matrixchi2[inucl][isotope][iene2][0]=1000.;//valore chi2 per questo z a questa energia
1162                            matrixchi2[inucl][isotope][iene2][1]=1000.;//energia per questo chi2
1163                            matrixchi2[inucl][isotope][iene2][2]=1000.;//piani saltati nel chi2
1164                    break;
1165        }        }
     }//fine loop energia  
1166    
1167        }//fine loop energia
1168    
1169      isotope++; //incremento il contatore isotopi
1170      }//fine loop isotopi
1171      isotope=0; //riazzero il contatore isotopi
1172    
1173   };//fine loop z    }//fine loop z
1174    
1175    for (Int_t nu=0; nu<colmax; nu++){    isotope=0;//non dovrebbe servire
1176      for (Int_t en=0; en<rowmax; en++){  
1177        if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] !=0.)){    //Emi
1178          bestchi2[0]= matrixchi2[nu][en][0];// chi2    for (Int_t nu=zmin; nu<max; nu++){
1179          bestchi2[1]= (Float_t)nu; // z      while(elem[nu-1][isotope]> 0.){
1180          bestchi2[2]= matrixchi2[nu][en][1];//energia;          for (Int_t en=0; en<nostep; en++){
1181          bestchi2[3]=    matrixchi2[nu][en][2];// totale piani saltati          if((matrixchi2[nu][isotope][en][0]<bestchi2[0]) && (matrixchi2[nu][isotope][en][0] >0.)){
1182              bestchi2[0]= matrixchi2[nu][isotope][en][0];// chi2
1183              bestchi2[1]= (Float_t)nu; // z
1184              bestchi2[2]= matrixchi2[nu][isotope][en][1];//energia;
1185              bestchi2[3]= matrixchi2[nu][isotope][en][2];// totale piani saltati
1186          bestchi2[4]= (Float_t)isotope; //isotopo
1187            }
1188        }        }
1189        isotope++;
1190      }      }
1191       isotope=0;
1192    }    }
   //==========================//  
   Int_t testz = (Int_t)bestchi2[1];  
1193    
   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;  
   }  
   //==========================//  
1194  };//endloopze  };//endloopze
1195    
1196    
1197    
1198    
1199    
1200  void CaloBragg::mediatroncata(){  // void CaloBragg::mediatroncata(){
1201    //calcolo Z con media troncata e utilizzo questo Z per trovare l'energia migliore  //   //calcolo Z con media troncata e utilizzo questo Z per trovare l'energia migliore
1202    //in: ordplane[44]   => array con energia dei piani  //   //in: ordplane[44]   => array con energia dei piani
1203    //    spess[3]    => conversioni spessore di silicio, w,  mip  //   //    spess[3]    => conversioni spessore di silicio, w,  mip
1204    //    estr[2][2]  => primo[0][0] e ultimo[1][0] piano attraversati ed energie[][1]  //   //    estr[2][2]  => primo[0][0] e ultimo[1][0] piano attraversati ed energie[][1]
1205    //    calo[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano  //   //    calo[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano
1206    //    integrale      => energia totale nel calorimetro considerando il W  //   //    integrale      => energia totale nel calorimetro considerando il W
1207    //  //   //
1208    // out[4] chi2,z,Etot,Pskip  //   // out[4] chi2,z,Etot,Pskip
1209    
1210    Float_t ordplane[44];//mi serve per la media troncata  //   Float_t ordplane[44];//mi serve per la media troncata
1211    memset(ordplane,0,44*sizeof(Float_t));  //   memset(ordplane,0,44*sizeof(Float_t));
1212    
1213    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
1214    
1215        
1216    //ordino tutte le energie dei piani in ordine crescente  //   //ordino tutte le energie dei piani in ordine crescente
1217    //Int_t ii=0;  
1218    //     Long64_t work[200];
1219      Long64_t work[200];  //     Int_t ind = 0;
1220      Int_t ind = 0;  //     //Int_t l = 0;
1221      Int_t l = 0;  //     Int_t RN = 0;
1222      Int_t RN = 0;  //     Float_t sum4 = 0.;
1223      Float_t sum4 = 0.;  //     Float_t qm = 0.;
1224      Float_t qm = 0.;  //     //
1225      //  //     //Float_t qmt = ethr*0.8; // *0.9
1226      //Float_t qmt = ethr*0.8; // *0.9  //     //
1227      //  //     //Int_t uplim = TMath::Max(3,N);
1228      //Int_t uplim = TMath::Max(3,N);  //     //
1229      //  //     while ( RN < 4 && ind < 44 ){
1230      while ( l < 4 && ind < 44 ){  //       qm = TMath::KOrdStat(44,ordplane,ind,work);
1231        qm = TMath::KOrdStat(44,ordplane,ind,work);  //       if (qm >= 0.7 ){      
1232        if (qm >= 0.7 ){    //      if ( RN < 4 ){
1233          if ( l < 4 ){  //        sum4 += qm;
1234            sum4 += qm;  //        RN++;
1235            RN++;  //      };
1236          };  // //   l++;
1237          l++;  // //   if ( debug ) printf(" value no %i qm %f sum4 %f \n",l,qm,sum4);
1238          if ( debug ) printf(" value no %i qm %f sum4 %f \n",l,qm,sum4);  //       };
1239        };  //       ind++;
1240        ind++;  //     };
1241      };  //     //
1242      //  //     sum4 /= (Float_t)RN;
1243      sum4 /= (Float_t)RN;  //     Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2])));//ma non e'/1??
1244      Float_t Zmean = (sqrt((sum4*MIP)/(4.*spessore[2])));  //     if(Zmean ==0.) Zmean=1.;
1245      if(Zmean ==0.) Zmean=1.;  //     if ( Zmean < 1. ) Zmean = 1.;
   
 //     cout<<"sum4="<<sum4<<endl;  
 //     cout<<"Zmean="<<Zmean<<endl;  
   
     
   //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.;  
   
   cout<<"inizio media troncata"<<endl;  
   cout<<"input: step "<<Integrale<<", zero "<<zero<<", zmin "<<zmin<<", zmax"<<zmin<<endl;  
   
   
     //    step  energia   zstart  zstop    
   loopze(Integrale,zero,zmin,zmin);  
     
   cout<<"fine media troncata"<<endl;  
   cout<<"output: chi2 "<<bestchi2[0]<<", z"<<bestchi2[1]<<", Etot "<<bestchi2[2]<<", Pskip "<<bestchi2[3]<<endl;  
1246    
1247        
1248    // //     Zmean =round(Zmean);
1249    // //     if(Zmean <1.) Zmean=1.;
1250    
1251    // //     if(Zmean >0.)Zmean =round(Zmean);
1252    
1253    qtchi2=bestchi2[0];  //     //======== per i nuclei=======
1254    qtz=bestchi2[1];  //     if (Zmean >=2.){
1255    qtetot=bestchi2[2];  //     ind = 0;
1256    qtpskip=bestchi2[3];  //     RN = 0;
1257  };//end mediatroncata  //     sum4 = 0.;
1258    //  cout<<"z media troncata ok"<<endl;  //     qm = 0.;
1259    //     while ( RN < 4 && ind < 44 ){
1260    //       qm = TMath::KOrdStat(44,ordplane,ind,work);
1261    //       if (qm >= (Zmean*Zmean)-Zmean*Zmean*0.2 ){    
1262    //      if ( RN < 4 ){
1263    //        sum4 += qm;
1264    //        RN++;
1265    //      };
1266    //       };
1267    //       ind++;
1268    //     };
1269    //     //
1270    //     sum4 /= (Float_t)RN;
1271    //     Zmean = (sqrt((sum4*MIP)/(4.*spessore[2])));//ma non e' /1??
1272    //     }
1273    
1274      
1275    //   //calcolo energia migliore per Z trovato con media troncata
1276    //   //  Float_t zmin=Zmean;
1277    //     Float_t zmin=round(Zmean);
1278      
1279    //   bestchi2[0]=10000.;
1280    //   bestchi2[1]=0.;
1281    //   bestchi2[2]=0.;
1282    //   bestchi2[3]=0.;
1283    //   Float_t zero=0.;
1284    
1285    //   //    step   energia zstart zstop  
1286    //   loopze(Integrale,zero,zmin,zmin);
1287      
1288    
1289    //   qtchi2=bestchi2[0];
1290    //   qtz=bestchi2[1];
1291    //   qtetot=bestchi2[2];
1292    //   qtpskip=bestchi2[3];
1293    // };//end mediatroncata
1294    
1295    
1296    
# Line 830  void CaloBragg::Zdaloop(){ Line 1306  void CaloBragg::Zdaloop(){
1306    
1307    
1308    /*z se particella fosse al minimo*/  //energia1piano/mip corretta    /*z se particella fosse al minimo*/  //energia1piano/mip corretta
1309    Float_t zmax = (sqrt(estremi[0][1]/spessore[2]));    //  Float_t zmax = round(sqrt(estremi[0][1]/spessore[2]));
1310    if(zmax<31)zmax=zmax+1;    //  if(zmax<31)zmax=zmax+1;
1311        
1312    /*calcolo Z ed E con loop sui vari elementi ed energie*/    /*calcolo Z ed E con loop sui vari elementi ed energie*/
1313    
1314    Float_t zmin=1.;    Float_t zmin=1.;
1315    Float_t bestchitemp[4] = {0,0,0,0};    Float_t zmax=32.;
1316      Float_t bestchitemp[5] = {0,0,0,0,0};
1317    
1318    bestchi2[0]=10000.;    bestchi2[0]=numeric_limits<Float_t>::max();
1319    bestchi2[1]=0.;    bestchi2[1]=0.;
1320    bestchi2[2]=0.;    bestchi2[2]=0.;
1321    bestchi2[3]=0.;    bestchi2[3]=0.;
1322      bestchi2[4]=0.;
1323    Float_t zero=0.;    Float_t zero=0.;
1324    //primo loop      //------------primo loop   ----------------------
1325    //     energia   ezero, zstart  zstop Si attrav  1 piano  energie piani  out    //     energia   ezero, zstart  zstop
1326      //  loopze(Integrale,zero,zmin,zmax);
1327    cout<<"inizio primo loop"<<endl;  
1328    cout<<"input: step "<<Integrale<<", zero "<<zero<<", zmin "<<zmin<<", zmax"<<zmax<<endl;    //->  loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,50);
1329      loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,200);
1330    loopze(Integrale,zero,zmin,zmax);  
1331      //  loopze(Integrale*2.,Integrale/100.,zmin,zmax);
1332    cout<<"fine primo loop"<<endl;    if ( debug) printf("Zdaloop start Integrale %f , outene %f \n",Integrale,bestchi2[2]);
1333    cout<<"output: chi2 "<<bestchi2[0]<<", z"<<bestchi2[1]<<", Etot "<<bestchi2[2]<<", Pskip "<<bestchi2[3]<<endl;    
1334        //------------secondo loop  ----------------------
1335    //secondo loop    for(Int_t i=0;i<5;i++) bestchitemp[i]=bestchi2[i];
1336    for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i];    bestchi2[0]=numeric_limits<Float_t>::max();
   bestchi2[0] = 10000.;  
1337    bestchi2[1] = 0.;    bestchi2[1] = 0.;
1338    bestchi2[2] = 0.;    bestchi2[2] = 0.;
1339    bestchi2[3] = 0.;//riazzero    bestchi2[3] = 0.;
1340      bestchi2[4] = 0.;//riazzero
1341    
1342    Float_t step = bestchitemp[2]/100.;//Integrale/10.;//era 1000    Float_t step = bestchitemp[2];//
1343    //zero=bestchitemp[2]-step/2.;//energia da 1 giro - 1step    zero=0.;  // qualsiasi altro valore peggiora le cose
1344    //Float_t step = bestchitemp[2]/1000.;    //  zmin=zmax=bestchitemp[1];
1345    //zero=bestchitemp[2]-step;    zmin=bestchitemp[1]-1;
1346    zero=0;    zmax=bestchitemp[1]+1;
1347    zmin=bestchitemp[1]-2;// zda primo giro    //loopze(step,zero,zmin,zmax); //
   if(zmin<1)zmin=1;  
   zmax=bestchitemp[1]+1;//  
1348    
1349    cout<<"inizio secondo loop"<<endl;    //->  loopze(step,step/2.,zmin,zmax,200); //
1350    cout<<"input: step "<<step<<", zero "<<zero<<", zmin "<<zmin<<", zmax"<<zmax<<endl;    loopze(step,step/2.,zmin,zmax,500); //
1351    
1352    loopze(step,zero,zmin,zmax); //    //step = bestchitemp[2];//
1353    
1354    cout<<"fine secondo loop"<<endl;    //loopze(step/2,3*step/4.,zmin,zmax,500); //
   cout<<"output: chi2 "<<bestchi2[0]<<", z"<<bestchi2[1]<<", Etot "<<bestchi2[2]<<", Pskip "<<bestchi2[3]<<endl;  
1355    
1356      if ( debug ) printf("Zdaloop Integrale2 %f , outene %f step %f \n",Integrale,bestchi2[2],step);
1357    
 //   cout<<"z loop ok"<<endl;  
     
1358    //chi2,z,Etot,Pskip    //chi2,z,Etot,Pskip
1359    lpchi2=bestchi2[0];    lpchi2=bestchi2[0];
1360    lpz=bestchi2[1];    lpz=bestchi2[1];
1361    lpetot=bestchi2[2];    lpetot=bestchi2[2];
1362    lppskip=bestchi2[3];    lppskip=bestchi2[3];
1363      lpisotope=bestchi2[4];
1364  };//endZdaloop  };//endZdaloop
1365    
1366    
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.23