/[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.12 by mocchiut, Fri Jun 17 11:56:42 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    lpetot = 0.;    lpetot = 0.;
47    lppskip = 0.;    lppskip = 0.;
48    memset(calorimetro,0,44*2*sizeof(Float_t));    memset(calorimetro,0,44*2*sizeof(Float_t));
49    memset(spessore,0,3*sizeof(Float_t));    memset(spessore,0,4*sizeof(Float_t));
50    memset(estremi,0,2*2*sizeof(Float_t));    memset(estremi,0,2*2*sizeof(Float_t));
51    Integrale=0.;    Integrale=0.;
52    
53      for(Int_t l=0;l<44;l++){
54        calorimetro[l][0]=-1.;
55      }
56    
   //  
57  };  };
58    
59  void CaloBragg::Print(){  void CaloBragg::Print(){
# Line 55  void CaloBragg::Print(){ Line 63  void CaloBragg::Print(){
63    //    //
64    printf("========================================================================\n");    printf("========================================================================\n");
65    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);
66    printf(" chi 2 from truncated mean: %f \n", qtchi2);    printf(" first plane: %f \n", estremi[0][0]);
67    printf(" Z from truncated mean %f: \n", qtz);    printf(" last plane: %f \n", estremi[1][0]);
68    printf(" energy from truncated mean %f: \n", qtetot);    //   printf(" chi 2 from truncated mean: %f \n", qtchi2);
69    printf(" plane not used for truncated mean %f: \n", qtpskip);    //   printf(" Z from truncated mean %f: \n", qtz);
70      //   printf(" energy from truncated mean %f: \n", qtetot);
71      //   printf(" plane not used for truncated mean %f: \n", qtpskip);
72    printf(" chi 2 from loop %f:  \n", lpchi2);    printf(" chi 2 from loop %f:  \n", lpchi2);
73    printf(" Z from loop %f: \n", lpz);    printf(" Z from loop %f: \n", lpz);
74    printf(" energy from loop %f: \n", lpetot);    printf(" energy from loop %f: \n", lpetot);
75    printf(" plane not used for loop %f: \n", lppskip);    printf(" plane not used for loop %f: \n", lppskip);
76      printf(" ndf: %i \n",ndf);
77    printf("========================================================================\n");    printf("========================================================================\n");
78    //    //
79  };  };
# Line 77  void CaloBragg::Process(){ Line 88  void CaloBragg::Process(){
88    Process(-1);    Process(-1);
89  };  };
90    
91    
92    void CaloBragg::CleanPlanes(Float_t epiano[22][2]){
93      //  return;
94      Int_t hitplanes = 0;
95      for (Int_t i = 0; i<22; i++){
96        for (Int_t j = 1; j>=0; j--){
97          if ( epiano[i][j] > 0.7 ) hitplanes++;
98        };
99      };
100      Float_t lowlim = 0.85;
101      Float_t dedxone = 0.;
102      Float_t step1 = 0.8*L2->GetCaloLevel2()->qtot/(Float_t)hitplanes;
103      while ( dedxone < step1 ){
104        for (Int_t i = 0; i<22; i++){
105          for (Int_t j = 1; j>=0; j--){
106            if ( epiano[i][j] >= step1 && dedxone < 0.7 ) dedxone = epiano[i][j];
107          };
108        };    
109      }
110      if ( dedxone < 0.7 ){
111        for (Int_t i = 0; i<22; i++){
112          for (Int_t j = 1; j>=0; j--){
113            if ( epiano[i][j] > 0. && dedxone < 0.7 ) dedxone = epiano[i][j];
114          };
115        };    
116      }
117      //
118      //  printf(" dedxone = %f step1 %f  \n",dedxone,step1);
119      Bool_t revulsera = false;
120      Bool_t nullius = false;
121      Int_t nulliferus = 0;
122      for (Int_t i = 0; i<22; i++){
123        for (Int_t j = 1; j>=0; j--){
124          if ( epiano[i][j] < dedxone*lowlim ){
125            //        printf(" %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius);
126            epiano[i][j] = 0.;
127          } else {
128            //x        printf(" %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius);
129            nulliferus = 0;
130            revulsera = true;
131          };
132          if ( epiano[i][j] < 0.7 && revulsera ) nulliferus++;
133          if ( nulliferus > 10 ) nullius = true;
134          if ( nullius ) epiano[i][j] = 0.;
135        };
136      };
137    
138    }
139    
140    
141  void CaloBragg::Process(Int_t ntr){  void CaloBragg::Process(Int_t ntr){
142    //      //  
143    if ( !L2 ){    if ( !L2 ){
# Line 109  void CaloBragg::Process(Int_t ntr){ Line 170  void CaloBragg::Process(Int_t ntr){
170    
171    //    //
172    //    //
   // Always calculate stdedx1  
173    //    //
174    Int_t view = 0;    Int_t view = 0;
175    Int_t plane = 0;    Int_t plane = 0;
# Line 120  void CaloBragg::Process(Int_t ntr){ Line 180  void CaloBragg::Process(Int_t ntr){
180    for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){    for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
181      //      //
182      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
183        //
184        if ( !usepl18x && view==0 && plane==18 ) mip = 0.;
185        //
186      epiano[plane][view]+=mip;      epiano[plane][view]+=mip;
187      //      //
188      //      //
189    };    };
190    //    //
191      this->CleanPlanes(*&epiano);
192    //    //
193    PamTrack *ptrack = 0;    PamTrack *ptrack = 0;
194    CaloTrkVar *track = 0;    CaloTrkVar *track = 0;
# Line 134  void CaloBragg::Process(Int_t ntr){ Line 198  void CaloBragg::Process(Int_t ntr){
198        ptrack = L2->GetTrack(ntr);        ptrack = L2->GetTrack(ntr);
199        if ( ptrack ) track = ptrack->GetCaloTrack();        if ( ptrack ) track = ptrack->GetCaloTrack();
200      } else {      } else {
201        track = L2->GetCaloStoredTrack(ntr);        track = L2->GetCaloStoredTrack(ntr);  //al momento e' vera solo questa riga
202      };      };
203      //      //
204      if ( !track && ntr >= 0 ){      if ( !track && ntr >= 0 ){
# Line 155  void CaloBragg::Process(Int_t ntr){ Line 219  void CaloBragg::Process(Int_t ntr){
219    //    //
220    for(Int_t p=0; p<22; p++){    for(Int_t p=0; p<22; p++){
221      for(Int_t v=0; v<2; v++){      for(Int_t v=0; v<2; v++){
       cout<< L2->GetCaloLevel2()->cibar[p][v];      
222        /*per usare traccia non del calo camboare cibar*/        /*per usare traccia non del calo camboare cibar*/
223        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      
224        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
225      };      };
226    };    };
227    
228      /*per ogni evento calcolo la conversione mip e w attraversato in equivalente Si*/
229      conversione(); // out: 1) g/cm2 Si , 2) spessoreW equivalente in Si, 3)Mip corretta per inclinazione
230    
231    /*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 */
232      Float_t ordplane[44];//mi serve per la media troncata
233      memset(ordplane,0,44*sizeof(Float_t));
234    
235    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
236    
237    cout<<"spessore"<<spessore[0]<<" gcm2si; "<<spessore[1]<<"W in Si; "<<spessore[2]<<"mip "<<endl;    
238      //ordino tutte le energie dei piani in ordine crescente
239    
240      Long64_t work[200];
241      Int_t ind = 0;
242      //Int_t l = 0;
243      Int_t RN = 0;
244      Float_t sum4 = 0.;
245      Float_t qm = 0.;
246      while ( RN < 4 && ind < 44 ){
247        qm = TMath::KOrdStat((Long64_t)44,ordplane,(Long64_t)ind,work);
248        if (qm >= 0.7 ){    
249          if ( RN < 4 ){
250            sum4 += qm;
251            RN++;
252          };
253        };
254        ind++;
255      };
256      //
257      //sum4 /= (Float_t)RN;
258      Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2])));
259      if(Zmean ==0.) Zmean=1.;
260      if ( Zmean < 1. ) Zmean = 1.;
261    
   /*trova primo e ultimo piano attraversati*/  
262    
263      /*trova primo e ultimo piano attraversati*/
264    Int_t p = 0;//contatore piani    Int_t p = 0;//contatore piani
265    //primo parte da 0 e va in giù    //per il primo parte da 0 e va in giu'
266    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;  
267      //    if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){      //    if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){
268      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
269        if( (calorimetro[p][0] >0) && (calorimetro[p][1]*MIP >Zmean*0.7)){ // 70% della MIP
270        estremi[0][0]=p;        estremi[0][0]=p;
271        estremi[0][1]=calorimetro[p][1] *MIP; //energia in MeV        estremi[0][1]=calorimetro[p][1] *MIP; //energia in MeV
272      };      };
273      p++;      p++;
274    };    };
275   //ultimo parte da 44 e sale  
276      //ultimo parte da 44 e sale
277    p=43;    p=43;
278    while( (estremi[1][1] == 0)  &&  (p>(int)estremi[0][0]) ){    while( (estremi[1][1] <= 0.)  &&  (p>(int)estremi[0][0]) ){
279    cout<<"piano "<<p <<"; strip "<<calorimetro[p][0]<<"; en "<<calorimetro[p][1]<<endl;      if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){
280    //    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;  
281        estremi[1][1]=calorimetro[p][1] *MIP;//energia in MeV        estremi[1][1]=calorimetro[p][1] *MIP;//energia in MeV
282      };      };
283      p = p-1;      p = p-1;
284    };    };
285    //    //
   cout<<"estremi: in "<<estremi[0][0]<<"piano, "<<estremi[0][1]<<"energia; out "<<estremi[1][0]<<"piano, "<<estremi[1][1]<<"energia."<<endl;  
286    
287      Float_t lastok = 0.;
288      //  Bool_t goback = false;
289      for ( int o = 0; o < estremi[1][0]; o++ ){
290        //
291        if ( calorimetro[o][1] > 0.7 ) lastok = calorimetro[o][1];
292        if ( calorimetro[o][1] < 0.7 && lastok > 0. ) calorimetro[o][1] = lastok;    
293        //    if ( calorimetro[o][1] < 0.7 ) goback = true;
294        //
295      };
296      lastok = 0.;
297      //  if ( goback ){
298      for ( int o = estremi[1][0]; o >= 0;  o-- ){
299        //
300        //    printf(" goback1: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);
301        if ( o < estremi[1][0] && calorimetro[o][1] > calorimetro[o+1][1]*1.2 && lastok > 0. ) calorimetro[o][1] = lastok;  
302        if ( calorimetro[o][1] > 0.7 ) lastok = calorimetro[o][1];
303        if ( calorimetro[o][1] < 0.7 && lastok > 0. ) calorimetro[o][1] = lastok;    
304        //    printf(" goback2: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);
305        //
306      };
307      //  };
308    
309      if ( startZero ) {
310        estremi[0][0] = 0.;
311        //    estremi[0][1] = 0.;
312      }
313        
314    /*integrale: energia rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/    /*integrale: energia totale rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/
315      for(Int_t pl=0; pl<(2*NPLA); pl++){
316    for(Int_t pl=0; pl<(2*NPLA); pl++){      //    printf(" integrale: calorimetro %f  \n",calorimetro[pl][1]);
   
317      //calcolo intergale in unita di spessori di silicio      //calcolo intergale in unita di spessori di silicio
318      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;  
319      // 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          
320      if(pl%2!=0){                                                             //equival W in Si      if(pl%2!=0){                                                              //equival W in Si
321        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;  
322      };      };
323    };    };
324      //Integrale=24000;//Integrale*1000;
325    //  cout<<"integrale ok"<<endl;    Integrale *= 1000.;
     
     
   /*z se particella fosse al minimo*/  //energia1piano/mip corretta  
   //Float_t zmax = (sqrt(estremi[0][1]/spessore[2]));  
     
   
326        
327    /*z ed energia con media troncata*/    /*z ed energia con media troncata*/
328    //  Float_t  bestchi2mean[4] = {0.,0.,0.,0.};//chi2,z,Etot,Pskip    //  mediatroncata();  // out:  1)chi2, 2)z, 3)Etot, 4)Pskip
   mediatroncata();    
329        
330    /*z ed energia con loop*/    /*z ed energia con loop*/
331    //  Float_t  bestchi2loop[4] = {0.,0.,0.,0.};//chi2,z,Etot,Pskip//chi2,z,Etot,Pskip    if ( debug ) printf(" call Zdaloop with integrale %f \n",Integrale);
332    Zdaloop();    Zdaloop(); // out:  1)chi2, 2)z, 3)Etot, 4)Pskip
     
   
   /*energia rilasciata da z migliore*/  
  //  Float_t dEpianimean[2*NPLA];  
   //  Int_t zet=(int)bestchi2mean[1];  
   //Enetrack(&zet, &bestchi2mean[2], &estremi[0][0], dEpianimean);//calcola rilascio energetico sui piani  
   
   // Float_t dEpianiloop[2*NPLA];  
   //zet=(int)bestchi2loop[1];  
   //Enetrack(&zet, &bestchi2loop[2], &estremi[0][0], dEpianiloop);//calcola rilascio energetico sui piani  
333    
334        
335    if ( debug ) this->Print();    if ( debug ) this->Print();
336    if ( debug ) printf(" esci \n");    if ( debug ) printf(" fine evento \n");
337    //    //
338  };  };
339    
340    
341    Float_t CaloBragg::Integral(){
342      Process();
343    
344      Float_t dEpianiloop[44];
345      Int_t tz1=(Int_t)lpz;
346      Enetrack(&tz1, &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop
347    
348    
349      Float_t integ = 0.;
350      for(Int_t i=0;i<=estremi[1][0];i++){
351        //    integ += dEplan[i];
352        //printf(" step %i integ %f deplan %f \n",i,integ,dEplan[i]);
353        integ += dEpianiloop[i];
354        //    printf(" step %i integ %f deplan %f \n",i,integ,dEpianiloop[i]);
355      }
356      return integ;
357    }
358    
359    Float_t CaloBragg::LastIntegral(){
360      Process();
361    
362      Float_t integ = 0.;
363      for(Int_t i=0;i<=estremi[1][0];i++){
364        integ += dEplan[i];
365        //printf(" step %i integ %f deplan %f \n",i,integ,dEplan[i]);
366      }
367      return integ;
368    }
369    
370    
371  void CaloBragg::Draw(){  void CaloBragg::Draw(){
372    
373    Process();    Process();
374    
375   Float_t dEpianimean[44];    this->Draw(0.,0.);
376   Float_t dEpianiloop[44];  
377   Float_t Depth[44];  }
378   Int_t tz=(Int_t)qtz;  
379   Int_t tz1=(Int_t)lpz;  void CaloBragg::Draw(Int_t Z, Float_t enetot){
380   Enetrack(&tz, &qtetot, &estremi[0][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata  
381   Enetrack(&tz1, &lpetot, &estremi[0][0], dEpianiloop);//calcola rilascio energetico sui piani da loop    //  Float_t dEpianimean[44];
382      Float_t dEpianiloop[44];
383      Float_t Depth[44];
384      //  Int_t tz=(Int_t)qtz;
385      Int_t tz1= Z;
386      Float_t enet = enetot;
387      //  Float_t enet = lpetot;
388      
389      if ( Z > 0. && enetot > 0. ){
390        estremi[0][0] = 0;
391        estremi[1][0] = 43;
392    
393    
394   Float_t sp= spessore[0]*spessore[1];      Float_t    ytgx = 0.;
395   for(Int_t i=0;i<44;i++)Depth[i]=i*sp;      Float_t    ytgy = 0.;
396        
397        //lunghezza effettiva di silicio attraversata (mm)
398        Float_t    SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy));
399        
400        spessore[0] = (SiCross/10.) * rhoSi; //spessore silicio in g/cm2
401        
402        /*tungsteno*/
403        
404        //rapporto tra rilasci energetici nei due materiali
405        Float_t WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad
406        //gcm2W = WCross/10. * rhoW;
407        
408        //       (g/cm2W)/(g/cm2Si)
409        spessore[3] =  (WCross/10.) * rhoW;
410        Float_t a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664);  //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664);
411        spessore[1] =  a;
412        //riscala mip allo spessore attraversato
413        spessore[2] = MIP*(SiCross/ySi);  
414    
415      } else {
416        tz1=(Int_t)lpz;
417        enet = lpetot;
418        //  Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata
419      
420      }
421      Enetrack(&tz1, &enet, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop
422    
423      Float_t sp= spessore[0]*spessore[1];
424      for(Int_t i=0;i<44;i++)Depth[i]=i*sp;
425    //    //
426    gStyle->SetLabelSize(0.04);    gStyle->SetLabelSize(0.04);
427    gStyle->SetNdivisions(510,"XY");    gStyle->SetNdivisions(510,"XY");
   //gStyle->SetLogy();  
428    //    //
429      TString hid = Form("cCaloBragg");    
430      TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
431      if ( tc ){
432        //   tc->Clear();
433      } else {
434        tc = new TCanvas(hid,hid);
435        //   tc->Divide(1,2);
436      };
437      //
438      //    TString thid = Form("hCaloBragg");      
439      //         TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
440      //    if ( th ) th->Delete();
441      //     th->Clear();
442      //     th->Reset();
443      //    } else {
444      //    th = new TH2F(thid,thid,300,-0.5,300.,1000,0.,150.);
445      //    th->SetMarkerStyle(20);
446      //    };
447      //
448      tc->cd();
449      TString thid2 = Form("hCaloBragg2");  
450      TH2F *th2  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid2));
451      if ( th2 ) th2->Delete();
452      th2 = new TH2F(thid2,thid2,300,-0.5,300.,1000,0.,150.);
453      th2->SetMarkerStyle(20);
454      th2->SetMarkerColor(kRed);
455      //
456      TString thid3 = Form("hCaloBragg3");  
457      TH2F *th3  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid3));
458      if ( th3 ) th3->Delete();
459      th3 = new TH2F(thid3,thid3,300,-0.5,300.,1000,0.,150.);
460      th3->SetMarkerStyle(20);
461      th3->SetMarkerColor(kBlue);
462    
463    
464      //  tc->cd(1);
465      //
466      //    for(Int_t i=0;i<=estremi[1][0];i++)th->Fill(Depth[i],dEpianimean[i]);
467      for(Int_t i=0;i<=estremi[1][0];i++)th2->Fill(Depth[i],calorimetro[i][1]*MIP);
468      //    th->Draw();
469      th2->Draw("same");
470    
471      //  tc->cd(2);
472      tc->cd();
473      //
474      for(Int_t i=0;i<=estremi[1][0];i++){
475        th3->Fill(Depth[i],dEpianiloop[i]);
476        //    printf(" i %i Depth %f depianiloop %f \n",i,Depth[i],dEpianiloop[i]);
477      }
478      th3->Draw();
479      th2->Draw("same");
480    
481          TString hid = Form("cCaloBragg");          tc->Modified();
482          TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));    tc->Update();
   
         if ( tc ){  
 //       tc->Clear();  
         } else {  
          tc = new TCanvas(hid,hid);  
          tc->Divide(1,2);  
         };  
         //  
         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<44;i++)th->Fill(Depth[i],dEpianimean[i]);  
         for(Int_t i=0;i<44;i++)th2->Fill(Depth[i],calorimetro[i][1]*MIP);  
         th->Draw();  
         th2->Draw("same");  
   
         tc->cd(2);  
         //  
         for(Int_t i=0;i<44;i++)th3->Fill(Depth[i],dEpianiloop[i]);  
         th3->Draw();  
         th2->Draw("same");  
   
         tc->Modified();  
         tc->Update();  
483    
484    //    //
485    gStyle->SetLabelSize(0);    gStyle->SetLabelSize(0);
# Line 325  void CaloBragg::Draw(){ Line 488  void CaloBragg::Draw(){
488  };  };
489    
490    
 void CaloBragg::SWAP( Float_t *A, Float_t *B ){  
    Float_t  Tmp = *A;  
    *A = *B;  
    *B = Tmp;  
 };  
   
   
491    
492  void CaloBragg::LoadParam(){  void CaloBragg::LoadParam(){
493    
494    //    //
495    elem[0] = 1.00794; //H  1    elem[0] = 1.00794; //H  1
496    elem[1] = 4.0026;  //He 2    elem[1] = 4.002602;  //He 2
497    elem[2] = 6.941;   //Li 3    elem[2] = 6.941;   //Li 3
498    elem[3] = 9.012182;//Be 4    elem[3] = 9.012182;//Be 4
499    elem[4] = 10.811;  //B  5    elem[4] = 10.811;  //B  5
500    elem[5] = 12.0107; //C  6    elem[5] = 12.0107; //C  6
501    elem[6] = 14.00674;//N  7    elem[6] = 14.00674;//N  7
502    elem[7] = 15.9994; //O  8    elem[7] = 15.9994; //O  8
503    elem[8] = 18.9984; //F  9    elem[8] = 18.9984032; //F  9
504    elem[9] = 20.1797; //Ne 10    elem[9] = 20.1797; //Ne 10
505    elem[10] = 22.98977;//Na 11    elem[10] = 22.98977;//Na 11
506    elem[11] = 24.3050; //Mg 12    elem[11] = 24.3050; //Mg 12
# Line 370  void CaloBragg::LoadParam(){ Line 526  void CaloBragg::LoadParam(){
526    elem[31] = 72.61; //Ge 32    elem[31] = 72.61; //Ge 32
527    
528    
529      //parametri calorimetro
 //parametri calorimetro  
530    NPLA = 22;    NPLA = 22;
531    NCHA = 96;    NCHA = 96;
532    nView = 2;    nView = 2;
# Line 392  void CaloBragg::LoadParam(){ Line 547  void CaloBragg::LoadParam(){
547    pigr = 3.1415;    pigr = 3.1415;
548    Na = 6.02e-23;    Na = 6.02e-23;
549    ZA = 0.49; /*Z/A per Si*/    ZA = 0.49; /*Z/A per Si*/
550    ISi =182e-06; /*MeV*/    //  ISi =182e-06; /*MeV*/
551      ISi = 171e-06; /*MeV*/
552      IW  = 735e-06; /*MeV*/
553      //  ISi =0.0001059994; /*GeV!!*/ no era giusto!!
554    Me = 0.511; /* MeV*/    Me = 0.511; /* MeV*/
555    MassP = 931.27;/*MeV*/    MassP = 931.27;/*MeV*/
556    r2 = 7.95e-26; /*ro*ro in cm */    r2 = 7.95e-26; /*ro*ro in cm */
# Line 420  void CaloBragg::conversione(){ Line 578  void CaloBragg::conversione(){
578    Float_t ytgy = 0;    Float_t ytgy = 0;
579    Float_t a = 0.;    Float_t a = 0.;
580    
   
581    /*silicio*/    /*silicio*/
582    ytgx = ySi * L2->GetCaloLevel2()->tanx[0];    ytgx = ySi * L2->GetCaloLevel2()->tanx[0];
583    ytgy = ySi * L2->GetCaloLevel2()->tany[0];    ytgy = ySi * L2->GetCaloLevel2()->tany[0];
# Line 428  void CaloBragg::conversione(){ Line 585  void CaloBragg::conversione(){
585    //lunghezza effettiva di silicio attraversata (mm)    //lunghezza effettiva di silicio attraversata (mm)
586    SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy));    SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy));
587        
588    spessore[0] = SiCross/10. * rhoSi; //spessore silicio in g/cm2    spessore[0] = (SiCross/10.) * rhoSi; //spessore silicio in g/cm2
     
589        
590    /*tungsteno*/    /*tungsteno*/
591    ytgx = yW * L2->GetCaloLevel2()->tanx[0];    ytgx = yW * L2->GetCaloLevel2()->tanx[0];
# Line 439  void CaloBragg::conversione(){ Line 595  void CaloBragg::conversione(){
595    WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad    WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad
596    //gcm2W = WCross/10. * rhoW;    //gcm2W = WCross/10. * rhoW;
597    
  a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664);  //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664);  
   
598    //       (g/cm2W)/(g/cm2Si)    //       (g/cm2W)/(g/cm2Si)
599      spessore[3] =  (WCross/10.) * rhoW;
600      a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664);  //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664);
601    spessore[1] =  a;    spessore[1] =  a;
   
     
602    //riscala mip allo spessore attraversato    //riscala mip allo spessore attraversato
603    spessore[2] = MIP*(SiCross/ySi);      spessore[2] = MIP*(SiCross/ySi);  
   
604  };//end conversione  };//end conversione
605    
606    
607    
608    
609    
610  void CaloBragg::BetheBloch(Float_t *x, Float_t *z, Float_t *Mass, Float_t *gam, Float_t *Bet, Float_t *out){  void CaloBragg::BetheBloch(Float_t *x, Float_t *z, Float_t *Mass, Float_t *gam, Float_t *Bet, Float_t *out, Float_t II){
611    
612    //rilascio energetico con bethe bloch con correzioni    //rilascio energetico con bethe bloch con correzioni
613    //in:    x: g/cm2    //in:    x: g/cm2
# Line 472  void CaloBragg::BetheBloch(Float_t *x, F Line 625  void CaloBragg::BetheBloch(Float_t *x, F
625    Float_t lg =0.;    Float_t lg =0.;
626    Float_t Energia=0.;    Float_t Energia=0.;
627    Float_t C=0.;    Float_t C=0.;
628      Float_t INo = ISi;
629    
630      if ( usenewBB ) INo = II;
631    
632    eta = (*gam)*(*Bet);    eta = (*gam)*(*Bet);
633    
634    //Bet=3/gam;  SQ(*gam) * SQ(*Bet)    //Bet=3/gam;  SQ(*gam) * SQ(*Bet)
635    Wmax = 2.* Me * SQ(eta) / (1. + 2.*(*gam)*Me/(*Mass) + SQ(Me)/SQ(*Mass));    Wmax = 2.* Me * SQ(eta) / (1. + 2.*(*gam)*Me/(*Mass) + SQ(Me)/SQ(*Mass));
636        
637    lg = 2.* Me * SQ(eta) * Wmax / SQ(ISi);    lg = 2.* Me * SQ(eta) * Wmax / SQ(INo);
638                                                                                                 //  Energia = x* 2 * pigr * Na * r2 * Me * rhoSi *ZA*  SQ(z)/SQ(Bet) * lg;    //  Energia = x* 2 * pigr * Na * r2 * Me * rhoSi *ZA*  SQ(z)/SQ(Bet) * lg;
639    C=(0.42237*pow(eta,-2) + 0.0304*pow(eta,-4) - 0.00038*pow(eta,-6))*pow(10,-6)* pow(ISi,2) +    C=(0.42237*pow(eta,-2.) + 0.0304*pow(eta,-4.) - 0.00038*pow(eta,-6.))*pow(10.,-6.)* pow(INo,2.) +
640      (3.858*pow(eta,-2) - 0.1668*pow(eta,-4) + 0.00158*pow(eta,-6))*pow(10,-9)*pow(ISi,3);      (3.858*pow(eta,-2.) - 0.1668*pow(eta,-4.) + 0.00158*pow(eta,-6.))*pow(10.,-9.)*pow(INo,3.);
641        
642    if(eta <= 0.13) C= C * log(eta/0.0653)/log(0.13/0.0653);          if(eta <= 0.13) C= C * log(eta/0.0653)/log(0.13/0.0653);      
   //   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.);  
     
643    
   /*se ho gia' lo spessore in g/cm2 non se mejo???*/  
644    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.);
645    
646    *out =Energia;//out    *out =Energia;//out
# Line 499  void CaloBragg::BetheBloch(Float_t *x, F Line 650  void CaloBragg::BetheBloch(Float_t *x, F
650    
651    
652    
653  void CaloBragg::ELOSS(Float_t *dx, Int_t *Z, Float_t *Etot, Float_t *out){  void CaloBragg::ELOSS(Float_t *dx, Int_t *Z, Float_t *Etot, Float_t *out, Float_t II){
654        
655    /*perdita di energia per ioni pesanti (come da routine geant)*/    /*perdita di energia per ioni pesanti (come da routine geant)*/
656    //  in : dx    => spessore g/cm2    //  in : dx    => spessore g/cm2
# Line 510  void CaloBragg::ELOSS(Float_t *dx, Int_t Line 661  void CaloBragg::ELOSS(Float_t *dx, Int_t
661    
662    
663    Float_t Q=0.;    Float_t Q=0.;
664    Float_t v=0.;    Float_t v=0.;                                                        
                                                         //  Float_t Mass=0.;  
665    Float_t gam=0.;    Float_t gam=0.;
666    Float_t Bet=0.;    Float_t Bet=0.;
667    Float_t dEP=0.;    Float_t dEP=0.;
# Line 522  void CaloBragg::ELOSS(Float_t *dx, Int_t Line 672  void CaloBragg::ELOSS(Float_t *dx, Int_t
672    
673    Bet = sqrt((SQ(gam) -1.)/SQ(gam));    Bet = sqrt((SQ(gam) -1.)/SQ(gam));
674        
675    v= 121.4139*(Bet/pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet/pow((*Z),(2./3.))));    //  v= 121.4139*(Bet/pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet/pow((*Z),(2./3.))));
676      v= 121.4139*(Bet*pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet*pow((*Z),(2./3.)))); // EMI AAAAGGH!!
677    
678    //carica effettiva    //carica effettiva
679    Q= (*Z)*(1- (1.034 - 0.1777*exp(-0.08114*(*Z)))*exp(-v));    Q= (*Z)*(1- (1.034 - 0.1777*exp(-0.08114*(*Z)))*exp(-v));
680    
681    //perdita energia per un protone    //perdita energia per un protone
682    Float_t protone =1.;    Float_t protone =1.;
683    Float_t Mass=MassP;    //  Float_t Mass=(elem[*Z-1]*MassP); //EMI
684    BetheBloch(dx, &protone, &Mass, &gam, &Bet, &dEP);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP);    //  BetheBloch(dx, &protone, &Mass, &gam, &Bet, &dEP);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP);
685    
686      BetheBloch(dx, &protone, &MassP, &gam, &Bet, &dEP, II);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP); //EMI
687        
688    *out= (SQ(Q)*(dEP));//*dx;    *out= (SQ(Q)*(dEP));//*dx;
689    
 };//end ELOSS  
   
690    
691    };//end ELOSS
692    
693    
694    
695    
696  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[]){
697    
698    //calcola energia rilasciata sulla traccia (usa ELOSS)    //calcola energia rilasciata sulla traccia (usa ELOSS)
699    //  in : Z             =>carica    //  in : Z             =>carica
# Line 560  void CaloBragg::Enetrack(Int_t* Z, Float Line 712  void CaloBragg::Enetrack(Int_t* Z, Float
712    memset(out, 0, 2*NPLA*sizeof(Float_t));    memset(out, 0, 2*NPLA*sizeof(Float_t));
713        
714    Float_t Massa = (elem[(*Z)-1] * MassP);    Float_t Massa = (elem[(*Z)-1] * MassP);
     
715    
716      for( Int_t ipla=((int)(*primo)); ipla<= ((int)(*ultimo)); ipla++){
   //loop piani      (dal primo in cui entra)  
   for( Int_t ipla=((int)(*primo)); ipla< (2*NPLA); ipla++){  
717      dE=0.;      dE=0.;
   
718      //spessore silicio corretto x inclinazione, z, energia, out:rilascio      //spessore silicio corretto x inclinazione, z, energia, out:rilascio
719      ELOSS(&spessore[0], Z, &Ezero, &dE);//spessore in g/cm2!!      ELOSS(&spessore[0], Z, &Ezero, &dE, ISi);//spessore in g/cm2!!
       
720      if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop      if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop
721        out[ipla] = Ezero - Massa; //MeV        out[ipla] = Ezero - Massa; //MeV
722        return;        return;
# Line 577  void CaloBragg::Enetrack(Int_t* Z, Float Line 724  void CaloBragg::Enetrack(Int_t* Z, Float
724      }else{      }else{
725        out[ipla] = dE; //MeV        out[ipla] = dE; //MeV
726        Ezero = Ezero - dE;//energia residua        Ezero = Ezero - dE;//energia residua
727          if ( debug ) printf(" zompa %i out %f dE %f ezero %f \n",ipla,out[ipla],dE,Ezero);
728      };      };
   
729      //se sono su un piano Y (tutti i pari) dopo c'e' il tungsteno      //se sono su un piano Y (tutti i pari) dopo c'e' il tungsteno
730      if(ipla%2 == 0){      if(ipla%2 == 0){
731        /*tungsteno*/        /*tungsteno*/
732        dE=0.;        dE=0.;
733        Float_t sp= spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato  in g/cm2        Float_t sp = 0.;
734        ELOSS(&sp, Z, &Ezero, &dE);        Float_t II = ISi;
735                if ( usenewBB ){
736            sp = spessore[3];        
737            II = IW;
738          } else {
739            sp = spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato  in g/cm2
740          }
741          //      printf(" sp %f II %f \n",sp,II);
742          ELOSS(&sp, Z, &Ezero, &dE,II);
743        if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop        if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop
744          return;          return;
745        }else{        }else{
746          Ezero = Ezero -dE;//energia residua          Ezero = Ezero -dE;//energia residua
         //      cout<<"w calc "<<dE<<endl;  
747        };        };
748      };          };    
749    
750    };//fine loop piani    };//fine loop piani
     
 };//end Enetrack  
751    
752    
753    };//end Enetrack
754    
755    
756    
# Line 613  void CaloBragg::chiquadro(Float_t dE[], Line 766  void CaloBragg::chiquadro(Float_t dE[],
766    
767    Float_t sum = 0.;    Float_t sum = 0.;
768    Float_t PianoPrecedente=0.;    Float_t PianoPrecedente=0.;
   //  Float_t differenza =0.;  
769    Float_t badplane=0.;    Float_t badplane=0.;
770    Float_t badplanetot=0.;    Float_t badplanetot=0.;
771    Float_t w,wi;    Float_t w,wi;
772      //
773    for(Int_t ipla=0; ipla<2*NPLA; ipla++){    if ( newchi2 ){
774          ndf = 0;
775      //tutti i piani attraversati dalla traiettoria      sum = 0.;
776      if(calorimetro[ipla][0] != -1.){      for( Int_t ipla=((int)(estremi[0][0])); ipla<= ((int)(estremi[1][0])); ipla++){
777        w=0.; //normalizzazione;        sum += pow((dE[ipla] - (calorimetro[ipla][1] * spessore[2]))/(0.05*dE[ipla]),2.);
778        wi=1.;//peso        //      printf(" quiqui: dE %f calor %f spessore[2] %f \n",dE[ipla],spessore[2]*calorimetro[ipla][1],spessore[2]);
779          ndf++;
780        }
781        ndf -= 2;
782        if ( ndf > 0 ) sum /= (float)ndf;
783        out[0] = sum;
784        out[1] = 0.;
785        out[2] = (int)(estremi[1][0])-ndf;
786        //    printf(" sum %f ndf %i \n ",sum,ndf);
787      } else {
788        for(Int_t ipla=0; ipla<2*NPLA; ipla++){
789          //tutti i piani attraversati dalla traiettoria
790          if(calorimetro[ipla][0] != -1.){ //
791            w=0.; //normalizzazione;
792            wi=1.;//peso
793    
794        //tolgo piani attraversati dalla traccia ma precedenti il piano individuato come ingresso          //tolgo piani attraversati dalla traccia ma precedenti il piano individuato come ingresso
795        if (ipla<estremi[0][0])  wi=0.;          if (ipla<estremi[0][0])  wi=0.;
796                
797        //tolgo piani attraversati da traccia ma successivi all'ultimo se sono diversi da 0          //tolgo piani attraversati da traccia ma successivi all'ultimo se sono diversi da 0
798        //    if((ipla>estremi[1][0]) && (calorimetro[ipla][1] >0.) ) wi=0.;          //if((ipla>estremi[1][0]) && (calorimetro[ipla][1] >0.) ) wi=0.;
799            if((ipla>estremi[1][0])) wi=0.;
800    
801        //normalizzazione          //normalizzazione      
802        if (calorimetro[ipla][1] != 0.)  w=1./(calorimetro[ipla][1]* MIP);                  if (calorimetro[ipla][1] != 0.)  w=1./(calorimetro[ipla][1]* MIP);    //        
803                
804        //tolgo piani con rilasci inferiori al 30% del precedente          //tolgo piani con rilasci inferiori al 30% del precedente
805        if(calorimetro[ipla][1] <= (0.7*PianoPrecedente)){          if(calorimetro[ipla][1] < (0.7*PianoPrecedente)){ // cosi' i piani senza rilascio non vengono considerati nel calcolo del chi2
806                      wi=0.;
807          wi=0.;            //se sono piani intermedi (non si e' fermta) li considero non buoni
808          //se sono piani intermedi (non si è fermta) li considero non buoni            if( (ipla <= estremi[1][0]) && (calorimetro[ipla][1] !=0.)){//
809          if( (ipla < estremi[1][0]) && (calorimetro[ipla][1] !=0.)){              badplane+=1.;
810            badplane+=1.;              badplanetot+=1.;
811            badplanetot+=1.;            };
812          };          };
813        };    
814                  //meno peso ai piani con rilasci maggiori di 1000 MIP
815            //      if(calorimetro[ipla][1] > 1000) wi=0.5;
816              if(calorimetro[ipla][1] > 1200.) wi=0.5;
       //meno peso ai piani con rilasci maggiori di 1000 MIP  
         if(calorimetro[ipla][1] > 1000) wi=0.5;  
817                            
818          //do peso maggiore alle ultime 6 misure          Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));
         //      if((ipla >= estr2[1][0]-12)) wi=0.;  
   
       Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));  
   
       sum += SQ(arg); // w*wi*(dEpiani[p][v]-(eplane[p][v]*MIP))));//( dEpiani[p][v] - (eplane[p][v]*MIP));  
819    
820        //se trovo piano non buono (tolto quindi wi=0) non modifico il piano precedente          sum += SQ(arg); // w*wi*(dEpiani[p][v]-(eplane[p][v]*MIP))));//( dEpiani[p][v] - (eplane[p][v]*MIP));
821        if(wi != 0.){          if(debug){
822          PianoPrecedente= calorimetro[ipla][1];//tengo piano precedente            printf("dedx  calcolata  %f e reale  %f  \n",dE[ipla],(calorimetro[ipla][1] * MIP));
823          badplane = 0.;//azzero contatore piani scartati consecutivi          }
824            //se trovo piano non buono (tolto quindi wi=0) non modifico il piano precedente
825            if(wi != 0.){//
826              PianoPrecedente= calorimetro[ipla][1];//tengo piano precedente
827              badplane = 0.;//azzero contatore piani scartati consecutivi
828            };
829        };        };
     };  
             
     if(badplane > 2) out[1] =79.;  
830                        
831    };//fine loop piani        //da Emi
832    //chi2,frammentato,pskip          if(badplane > 2){      
833    out[0]=sum;          //      printf(" AAAAAAAAAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG\n");
834            out[1] =79.;
835    out[2]=badplanetot;          break;
836          };
837        
838        };//fine loop piani
839        //chi2,frammentato,pskip  
840        out[0]=sum;
841        out[2]=badplanetot;
842      }  
843  };//end chiquadro  };//end chiquadro
844    
845    
846    
847  //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){
848  void CaloBragg::loopze( Float_t step, Float_t E0,Float_t Zstart, Float_t Zlimite){    //
 //  
849    //loop su z ed energie per trovare miglior z (ed energia)    //loop su z ed energie per trovare miglior z (ed energia)
850    //in:  nloop     => energia massima da provare (nloop x E0)    //in:  nloop     => energia massima da provare (nloop x E0)
851    //     E0        => energia iniziale (intergale)            //     E0        => energia iniziale (intergale)        
852    //     Zstart    => minimo z da cui patire    //     Zstart    => minimo z da cui patire
853    //     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  
854    //    //
855    //out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2    //out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2
856    //    //
857    
858    Float_t dEplan[2*NPLA];//energia rilasciata calcolata  
859    memset(dEplan,0,2*NPLA*sizeof(Float_t));    memset(dEplan,0,2*NPLA*sizeof(Float_t));
860    
861    Int_t Z = 0;// z iniziale    Int_t Z = 0;// z iniziale
862    
863    Float_t Massa = 0.;    Float_t Massa = 0.;
864    
865    Float_t Stepint =(step)/1000.;//passo per il calcolo di energia    Float_t Stepint =(step)/(Float_t)nostep;//passo per il calcolo di energia
866    
867    Float_t energia =0.;//energia del loop    Float_t energia =0.;//energia del loop
868    
869    Float_t chi2[3] = {0,0,0};//out dal calcolo chi2: chi2, piani consecutivi saltati, piani totali saltati    Float_t chi2[3] = {0,0,0};//out dal calcolo chi2: chi2, piani consecutivi saltati, piani totali saltati
870      
871      Int_t zmin = (int)Zstart;  
872    Int_t max=32;//max z di cui so la massa :P    Int_t max=32;//max z di cui so la massa :P
873    if((Zlimite)<=31) max=(int)(Zlimite) + 1;    if((Zlimite)<=31) max=(int)(Zlimite) + 1;
874    
875      if ( fzeta > 0. ){
876        zmin = fzeta;
877        max = fzeta+1;
878      }
879        
880    Int_t colmax=32;    Int_t colmax=32;
881    Int_t rowmax=3000;    Int_t rowmax=3000;
# Line 713  void CaloBragg::loopze( Float_t step, Fl Line 883  void CaloBragg::loopze( Float_t step, Fl
883    Float_t matrixchi2[colmax][rowmax][3];    Float_t matrixchi2[colmax][rowmax][3];
884    memset(matrixchi2, 0, colmax*rowmax*3*sizeof(Float_t));    memset(matrixchi2, 0, colmax*rowmax*3*sizeof(Float_t));
885    
886      Int_t imin = 1-nostep/2;
887      Int_t imax = nostep/2;
888    
889    //loop elementi      //loop elementi  
890    for(Int_t inucl=(int)(Zstart); inucl<max; inucl++){    for(Int_t inucl=zmin; inucl<max; inucl++){
891            
892      Z= inucl;      Z= inucl;
893        
894      Massa = elem[inucl-1]*MassP;      Massa = elem[inucl-1]*MassP;
895            
896      //loop energia      //loop energia
897      for(Int_t iene= 0; iene<3000; iene++){      Int_t iene2 = 0;
898            //    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 <=??
899        energia=  Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica..      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 <=??
       //cout<<"folse"<<estremi[0][0];  
       Enetrack(&Z, &energia, &estremi[0][0], dEplan);//calcola rilascio energetico sui piani  
900                
901          iene2++;
902          energia=  Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica..
903    
904          Enetrack(&Z, &energia, &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani
905        //calcolo chi2        //calcolo chi2
906        chiquadro(dEplan,chi2);        chiquadro(dEplan,chi2);
907    
908          //      printf(" last deplan from: Z = %i iene %i energia %f chi2 %f \n",inucl,iene,energia,chi2[0]);
909                
910        if( (chi2[1] != 79) ){//salto quelli che frammentano        if( (chi2[1] != 79.) ){//salto quelli che frammentano
911          matrixchi2[inucl][iene][0]=chi2[0];//valore chi2 per questo z a questa energia          matrixchi2[inucl][iene2][0]=chi2[0];//valore chi2 per questo z a questa energia
912          matrixchi2[inucl][iene][1]=energia;//energia per questo chi2          matrixchi2[inucl][iene2][1]=energia;//energia per questo chi2
913          matrixchi2[inucl][iene][2]=chi2[2];//piani saltati ne chi2          matrixchi2[inucl][iene2][2]=chi2[2];//piani saltati nel chi2
914        } else {          } else {  
915          matrixchi2[inucl][iene][0]=1000;//valore chi2 per questo z a questa energia          matrixchi2[inucl][iene2][0]=1000.;//valore chi2 per questo z a questa energia
916          matrixchi2[inucl][iene][1]=1000;//energia per questo chi2          matrixchi2[inucl][iene2][1]=1000.;//energia per questo chi2
917          matrixchi2[inucl][iene][2]=1000;//piani saltati ne chi2          matrixchi2[inucl][iene2][2]=1000.;//piani saltati nel chi2
918            break;
919        }        }
920      }//fine loop energia      }//fine loop energia
921    
922    
923   };//fine loop z    };//fine loop z
924    
925    
926    for (Int_t nu=0; nu<colmax; nu++){    //Emi
927      for (Int_t en=0; en<rowmax; en++){    for (Int_t nu=zmin; nu<max; nu++){
928        if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] !=0.)){      for (Int_t en=0; en<nostep; en++){
929          if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] >0.)){
930          bestchi2[0]= matrixchi2[nu][en][0];// chi2          bestchi2[0]= matrixchi2[nu][en][0];// chi2
931          bestchi2[1]= (Float_t)nu;          bestchi2[1]= (Float_t)nu; // z
932          bestchi2[2]= matrixchi2[nu][en][1];//energia;          bestchi2[2]= matrixchi2[nu][en][1];//energia;
933          bestchi2[3]=    matrixchi2[nu][en][2];// totale piani saltati          bestchi2[3]= matrixchi2[nu][en][2];// totale piani saltati
934        }        }
935      }      }
936    }    }
937    
 };//endloopze  
   
938    
939    };//endloopze
940    
941    
942    
 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));  
943    
   for(Int_t ipla=0; ipla< 2*NPLA; ipla++)  ordplane[ipla]=calorimetro[ipla][1]; //energia del piano  
944    
945      // void CaloBragg::mediatroncata(){
946    //ordino tutte le energie dei piani in ordine crescente  //   //calcolo Z con media troncata e utilizzo questo Z per trovare l'energia migliore
947    Int_t ii=0;  //   //in: ordplane[44]   => array con energia dei piani
948  //   while( ii < NPLA-1 ){    //   //    spess[3]    => conversioni spessore di silicio, w,  mip
949        //   //    estr[2][2]  => primo[0][0] e ultimo[1][0] piano attraversati ed energie[][1]
950  //     if(ordplane[ii+1] >= ordplane[ii]){  //   //    calo[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano
951  //       ii++;  //   //    integrale      => energia totale nel calorimetro considerando il W
952  //     }else{  //   //
953  //       SWAP( &(ordplane[ii]), &(ordplane[ii+1]));  //   // out[4] chi2,z,Etot,Pskip
954  //       ii=0;  
955    //   Float_t ordplane[44];//mi serve per la media troncata
956    //   memset(ordplane,0,44*sizeof(Float_t));
957    
958    //   for(Int_t ipla=0; ipla< 2*NPLA; ipla++)  ordplane[ipla]=calorimetro[ipla][1]; //energia del piano
959    
960      
961    //   //ordino tutte le energie dei piani in ordine crescente
962    
963    //     Long64_t work[200];
964    //     Int_t ind = 0;
965    //     //Int_t l = 0;
966    //     Int_t RN = 0;
967    //     Float_t sum4 = 0.;
968    //     Float_t qm = 0.;
969    //     //
970    //     //Float_t qmt = ethr*0.8; // *0.9
971    //     //
972    //     //Int_t uplim = TMath::Max(3,N);
973    //     //
974    //     while ( RN < 4 && ind < 44 ){
975    //       qm = TMath::KOrdStat(44,ordplane,ind,work);
976    //       if (qm >= 0.7 ){      
977    //      if ( RN < 4 ){
978    //        sum4 += qm;
979    //        RN++;
980    //      };
981    // //   l++;
982    // //   if ( debug ) printf(" value no %i qm %f sum4 %f \n",l,qm,sum4);
983    //       };
984    //       ind++;
985  //     };  //     };
986  //   };  //     //
987    //     sum4 /= (Float_t)RN;
988    //     Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2])));//ma non e'/1??
989    //     if(Zmean ==0.) Zmean=1.;
990    //     if ( Zmean < 1. ) Zmean = 1.;
991    
992        
993    // //     Zmean =round(Zmean);
994    // //     if(Zmean <1.) Zmean=1.;
995    
996    //scelgo 4 piani minimo 'sensati'  // //     if(Zmean >0.)Zmean =round(Zmean);
   Float_t sum4=0.;  
   Int_t pi=0;  
   cout<<"sum4="<<sum4<<endl;  
   while(pi<2*NPLA){  
     if ((ordplane[pi] >0.7)&& sum4==0.){// (ordplane[pi] >50)  
       sum4=(ordplane[pi] + ordplane[pi+1] + ordplane[pi+2] + ordplane[pi+3]);  
       pi=2*NPLA;  
     };  
     pi++;  
   };  
   cout<<"sum4="<<sum4<<endl;  
   Float_t Zmean = (sqrt((sum4*MIP)/(4.*spessore[2])));  
   cout<<"Zmean="<<Zmean<<endl;  
   if(Zmean ==0.) Zmean=1.;  
     
   /*z se particella fosse al minimo*/  //energia1piano/mip corretta  
   //  Float_t zmax = (sqrt(estremi[0][1]/spessore[2]));  
997    
998      //     //======== per i nuclei=======
999    //calcolo energia migliore per Z trovato con media troncata  //     if (Zmean >=2.){
1000    //Int_t step = 0;  //     ind = 0;
1001    Float_t zmin=Zmean;  //     RN = 0;
1002    //  Float_t bestchi2[4] = {10000,0,0,0};//chi2,z,Etot,Pskip //era bestchi2[5]  //     sum4 = 0.;
1003    bestchi2[0]=10000.;  //     qm = 0.;
1004    bestchi2[1]=0.;  //     while ( RN < 4 && ind < 44 ){
1005    bestchi2[2]=0.;  //       qm = TMath::KOrdStat(44,ordplane,ind,work);
1006    bestchi2[3]=0.;  //       if (qm >= (Zmean*Zmean)-Zmean*Zmean*0.2 ){    
1007    Float_t zero=0.;  //      if ( RN < 4 ){
1008      //    step  energia   zstart  zstop Si attrav  1 piano  energie piani  out  //        sum4 += qm;
1009    loopze(Integrale,zero,zmin,zmin);  //        RN++;
1010      //      };
1011    qtchi2=bestchi2[0];  //       };
1012    qtz=bestchi2[1];  //       ind++;
1013    qtetot=bestchi2[2];  //     };
1014    qtpskip=bestchi2[3];  //     //
1015  };//end mediatroncata  //     sum4 /= (Float_t)RN;
1016    //  cout<<"z media troncata ok"<<endl;  //     Zmean = (sqrt((sum4*MIP)/(4.*spessore[2])));//ma non e' /1??
1017    //     }
1018    
1019      
1020    //   //calcolo energia migliore per Z trovato con media troncata
1021    //   //  Float_t zmin=Zmean;
1022    //     Float_t zmin=round(Zmean);
1023      
1024    //   bestchi2[0]=10000.;
1025    //   bestchi2[1]=0.;
1026    //   bestchi2[2]=0.;
1027    //   bestchi2[3]=0.;
1028    //   Float_t zero=0.;
1029    
1030    //   //    step   energia zstart zstop  
1031    //   loopze(Integrale,zero,zmin,zmin);
1032      
1033    
1034    //   qtchi2=bestchi2[0];
1035    //   qtz=bestchi2[1];
1036    //   qtetot=bestchi2[2];
1037    //   qtpskip=bestchi2[3];
1038    // };//end mediatroncata
1039    
1040    
1041    
# Line 839  void CaloBragg::Zdaloop(){ Line 1047  void CaloBragg::Zdaloop(){
1047    //    calo1[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano    //    calo1[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano
1048    //    integrale=> energia totale nel calorimetro considerando il W    //    integrale=> energia totale nel calorimetro considerando il W
1049    //    //
1050    // out[6] chi2,z,Etot,Pskip    // out[4] chi2,z,Etot,Pskip
1051    
1052    
1053    /*z se particella fosse al minimo*/  //energia1piano/mip corretta    /*z se particella fosse al minimo*/  //energia1piano/mip corretta
1054    Float_t zmax = (sqrt(estremi[0][1]/spessore[2]));    //  Float_t zmax = round(sqrt(estremi[0][1]/spessore[2]));
1055      //  if(zmax<31)zmax=zmax+1;
1056        
1057    /*calcolo Z ed E con loop sui vari elementi ed energie*/    /*calcolo Z ed E con loop sui vari elementi ed energie*/
1058    //Int_t step = 0;  
1059    Float_t zmin=1.;    Float_t zmin=1.;
1060      Float_t zmax=32.;
1061    Float_t bestchitemp[4] = {0,0,0,0};    Float_t bestchitemp[4] = {0,0,0,0};
1062    //  Float_t bestchi2[4] = {10000,0,0,0};//chi2,z,Etot,Pskip  
1063    bestchi2[0]=10000.;    bestchi2[0]=10000.;
1064    bestchi2[1]=0.;    bestchi2[1]=0.;
1065    bestchi2[2]=0.;    bestchi2[2]=0.;
1066    bestchi2[3]=0.;    bestchi2[3]=0.;
1067    Float_t zero=0.;    Float_t zero=0.;
1068    //primo loop      //------------primo loop   ----------------------
1069    //     energia   ezero, zstart  zstop Si attrav  1 piano  energie piani  out    //     energia   ezero, zstart  zstop
1070    loopze(Integrale,zero,zmin,zmax);    //  loopze(Integrale,zero,zmin,zmax);
1071    
1072      //->  loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,50);
1073      loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,200);
1074    
1075      //  loopze(Integrale*2.,Integrale/100.,zmin,zmax);
1076      if ( debug ) printf(" Integrale %f , outene %f \n",Integrale,bestchi2[2]);
1077        
1078    //secondo loop    //------------secondo loop  ----------------------
1079    for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i];    for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i];
1080    bestchi2[0] = 10000.;    bestchi2[0] = 10000.;
1081    bestchi2[1] = 0.;    bestchi2[1] = 0.;
1082    bestchi2[2] = 0.;    bestchi2[2] = 0.;
1083    bestchi2[3] = 0.;//riazzero    bestchi2[3] = 0.;//riazzero
1084    
1085    Float_t step = Integrale/1000.;    Float_t step = bestchitemp[2];//
1086    zero=bestchitemp[2]-step;//energia da 1 giro - 1step    zero=0.;  // qualsiasi altro valore peggiora le cose
1087    zmin=bestchitemp[1]-2;// zda primo giro    //  zmin=zmax=bestchitemp[1];
1088    if(zmin<1)zmin=1;    zmin=bestchitemp[1]-1;
1089    zmax=bestchitemp[1]+1;    zmax=bestchitemp[1]+1;
1090      //  loopze(step,zero,zmin,zmax); //
1091    
1092      //->  loopze(step,step/2.,zmin,zmax,200); //
1093      loopze(step,step/2.,zmin,zmax,500); //
1094    
1095    loopze(step,zero,zmin,zmax);    if ( debug ) printf(" Integrale2 %f , outene %f step %f \n",Integrale,bestchi2[2],step);
1096    
   //  cout<<"z loop ok"<<endl;  
1097        
1098    //chi2,z,Etot,Pskip    //chi2,z,Etot,Pskip
1099    lpchi2=bestchi2[0];    lpchi2=bestchi2[0];
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.23