/[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.7 by mocchiut, Mon Dec 14 14:37:58 2009 UTC revision 1.8 by mocchiut, Fri Jan 28 10:40:44 2011 UTC
# Line 33  void CaloBragg::Clear(){ Line 33  void CaloBragg::Clear(){
33    //    //
34    tr = 0;    tr = 0;
35    sntr = 0;    sntr = 0;
36  //   qtchi2 = 0.;    //   qtchi2 = 0.;
37  //   qtz = 0.;    //   qtz = 0.;
38  //   qtetot = 0.;    //   qtetot = 0.;
39  //   qtpskip = 0.;    //   qtpskip = 0.;
40    lpchi2 = 0.;    lpchi2 = 0.;
41    lpz = 0.;    lpz = 0.;
42    lpetot = 0.;    lpetot = 0.;
# Line 61  void CaloBragg::Print(){ Line 61  void CaloBragg::Print(){
61    printf(" OBT: %u PKT: %u ATIME: %u Track %i Use track %i \n",OBT,PKT,atime,tr,usetrack);    printf(" OBT: %u PKT: %u ATIME: %u Track %i Use track %i \n",OBT,PKT,atime,tr,usetrack);
62    printf(" first plane: %f \n", estremi[0][0]);    printf(" first plane: %f \n", estremi[0][0]);
63    printf(" last plane: %f \n", estremi[1][0]);    printf(" last plane: %f \n", estremi[1][0]);
64  //   printf(" chi 2 from truncated mean: %f \n", qtchi2);    //   printf(" chi 2 from truncated mean: %f \n", qtchi2);
65  //   printf(" Z from truncated mean %f: \n", qtz);    //   printf(" Z from truncated mean %f: \n", qtz);
66  //   printf(" energy from truncated mean %f: \n", qtetot);    //   printf(" energy from truncated mean %f: \n", qtetot);
67  //   printf(" plane not used for truncated mean %f: \n", qtpskip);    //   printf(" plane not used for truncated mean %f: \n", qtpskip);
68    printf(" chi 2 from loop %f:  \n", lpchi2);    printf(" chi 2 from loop %f:  \n", lpchi2);
69    printf(" Z from loop %f: \n", lpz);    printf(" Z from loop %f: \n", lpz);
70    printf(" energy from loop %f: \n", lpetot);    printf(" energy from loop %f: \n", lpetot);
# Line 83  void CaloBragg::Process(){ Line 83  void CaloBragg::Process(){
83    Process(-1);    Process(-1);
84  };  };
85    
86    
87    void CaloBragg::CleanPlanes(Float_t epiano[22][2]){
88      //  return;
89      Int_t hitplanes = 0;
90      for (Int_t i = 0; i<22; i++){
91        for (Int_t j = 1; j>=0; j--){
92          if ( epiano[i][j] > 0.7 ) hitplanes++;
93        };
94      };
95      Float_t lowlim = 0.85;
96      Float_t dedxone = 0.;
97      Float_t step1 = 0.8*L2->GetCaloLevel2()->qtot/(Float_t)hitplanes;
98      while ( dedxone < step1 ){
99        for (Int_t i = 0; i<22; i++){
100          for (Int_t j = 1; j>=0; j--){
101            if ( epiano[i][j] >= step1 && dedxone < 0.7 ) dedxone = epiano[i][j];
102          };
103        };    
104      }
105      if ( dedxone < 0.7 ){
106        for (Int_t i = 0; i<22; i++){
107          for (Int_t j = 1; j>=0; j--){
108            if ( epiano[i][j] > 0. && dedxone < 0.7 ) dedxone = epiano[i][j];
109          };
110        };    
111      }
112      //
113      //  printf(" dedxone = %f step1 %f  \n",dedxone,step1);
114      Bool_t revulsera = false;
115      Bool_t nullius = false;
116      Int_t nulliferus = 0;
117      for (Int_t i = 0; i<22; i++){
118        for (Int_t j = 1; j>=0; j--){
119          if ( epiano[i][j] < dedxone*lowlim ){
120            //        printf(" %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius);
121            epiano[i][j] = 0.;
122          } else {
123            //x        printf(" %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius);
124            nulliferus = 0;
125            revulsera = true;
126          };
127          if ( epiano[i][j] < 0.7 && revulsera ) nulliferus++;
128          if ( nulliferus > 10 ) nullius = true;
129          if ( nullius ) epiano[i][j] = 0.;
130        };
131      };
132    
133    }
134    
135    
136  void CaloBragg::Process(Int_t ntr){  void CaloBragg::Process(Int_t ntr){
137    //      //  
138    if ( !L2 ){    if ( !L2 ){
# Line 133  void CaloBragg::Process(Int_t ntr){ Line 183  void CaloBragg::Process(Int_t ntr){
183      //      //
184    };    };
185    //    //
186      this->CleanPlanes(*&epiano);
187    //    //
188    PamTrack *ptrack = 0;    PamTrack *ptrack = 0;
189    CaloTrkVar *track = 0;    CaloTrkVar *track = 0;
# Line 164  void CaloBragg::Process(Int_t ntr){ Line 215  void CaloBragg::Process(Int_t ntr){
215    for(Int_t p=0; p<22; p++){    for(Int_t p=0; p<22; p++){
216      for(Int_t v=0; v<2; v++){      for(Int_t v=0; v<2; v++){
217        /*per usare traccia non del calo camboare cibar*/        /*per usare traccia non del calo camboare cibar*/
218        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      
219        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
220      };      };
221    };    };
222    
# Line 181  void CaloBragg::Process(Int_t ntr){ Line 232  void CaloBragg::Process(Int_t ntr){
232        
233    //ordino tutte le energie dei piani in ordine crescente    //ordino tutte le energie dei piani in ordine crescente
234    
235      Long64_t work[200];    Long64_t work[200];
236      Int_t ind = 0;    Int_t ind = 0;
237      //Int_t l = 0;    //Int_t l = 0;
238      Int_t RN = 0;    Int_t RN = 0;
239      Float_t sum4 = 0.;    Float_t sum4 = 0.;
240      Float_t qm = 0.;    Float_t qm = 0.;
241      while ( RN < 4 && ind < 44 ){    while ( RN < 4 && ind < 44 ){
242        qm = TMath::KOrdStat((Long64_t)44,ordplane,(Long64_t)ind,work);      qm = TMath::KOrdStat((Long64_t)44,ordplane,(Long64_t)ind,work);
243        if (qm >= 0.7 ){        if (qm >= 0.7 ){    
244          if ( RN < 4 ){        if ( RN < 4 ){
245            sum4 += qm;          sum4 += qm;
246            RN++;          RN++;
         };  
247        };        };
       ind++;  
248      };      };
249      //      ind++;
250      //sum4 /= (Float_t)RN;    };
251      Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2])));    //
252      if(Zmean ==0.) Zmean=1.;    //sum4 /= (Float_t)RN;
253      if ( Zmean < 1. ) Zmean = 1.;    Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2])));
254      if(Zmean ==0.) Zmean=1.;
255      if ( Zmean < 1. ) Zmean = 1.;
256    
257    
258    /*trova primo e ultimo piano attraversati*/    /*trova primo e ultimo piano attraversati*/
# Line 216  void CaloBragg::Process(Int_t ntr){ Line 267  void CaloBragg::Process(Int_t ntr){
267      };      };
268      p++;      p++;
269    };    };
270   //ultimo parte da 44 e sale  
271      //ultimo parte da 44 e sale
272    p=43;    p=43;
273    while( (estremi[1][1] <= 0.)  &&  (p>(int)estremi[0][0]) ){    while( (estremi[1][1] <= 0.)  &&  (p>(int)estremi[0][0]) ){
274      if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){      if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){
# Line 226  void CaloBragg::Process(Int_t ntr){ Line 278  void CaloBragg::Process(Int_t ntr){
278      p = p-1;      p = p-1;
279    };    };
280    //    //
281    
282      Float_t lastok = 0.;
283      //  Bool_t goback = false;
284      for ( int o = 0; o < estremi[1][0]; o++ ){
285        //
286        if ( calorimetro[o][1] > 0.7 ) lastok = calorimetro[o][1];
287        if ( calorimetro[o][1] < 0.7 && lastok > 0. ) calorimetro[o][1] = lastok;    
288        //    if ( calorimetro[o][1] < 0.7 ) goback = true;
289        //
290      };
291      lastok = 0.;
292      //  if ( goback ){
293      for ( int o = estremi[1][0]; o >= 0;  o-- ){
294        //
295        //    printf(" goback1: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);
296        if ( o < estremi[1][0] && calorimetro[o][1] > calorimetro[o+1][1]*1.2 && lastok > 0. ) calorimetro[o][1] = lastok;  
297        if ( calorimetro[o][1] > 0.7 ) lastok = calorimetro[o][1];
298        if ( calorimetro[o][1] < 0.7 && lastok > 0. ) calorimetro[o][1] = lastok;    
299        //    printf(" goback2: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);
300        //
301      };
302      //  };
303    
304    
305        
306    /*integrale: energia totale rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/    /*integrale: energia totale rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/
307    for(Int_t pl=0; pl<(2*NPLA); pl++){    for(Int_t pl=0; pl<(2*NPLA); pl++){
308        //    printf(" integrale: calorimetro %f  \n",calorimetro[pl][1]);
309      //calcolo intergale in unita di spessori di silicio      //calcolo intergale in unita di spessori di silicio
310      Integrale += calorimetro[pl][1] * MIP;//piano di silicio      Integrale += calorimetro[pl][1] * MIP;//piano di silicio
311      // 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          
# Line 236  void CaloBragg::Process(Int_t ntr){ Line 313  void CaloBragg::Process(Int_t ntr){
313        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]);
314      };      };
315    };    };
316    Integrale=24000;//Integrale*1000;    //Integrale=24000;//Integrale*1000;
317      Integrale *= 1000.;
318        
319    /*z ed energia con media troncata*/    /*z ed energia con media troncata*/
320    //  mediatroncata();  // out:  1)chi2, 2)z, 3)Etot, 4)Pskip    //  mediatroncata();  // out:  1)chi2, 2)z, 3)Etot, 4)Pskip
321        
322    /*z ed energia con loop*/    /*z ed energia con loop*/
323      Zdaloop(); // out:  1)chi2, 2)z, 3)Etot, 4)Pskip    Zdaloop(); // out:  1)chi2, 2)z, 3)Etot, 4)Pskip
324    
325        
326    if ( debug ) this->Print();    if ( debug ) this->Print();
# Line 255  void CaloBragg::Draw(){ Line 333  void CaloBragg::Draw(){
333    
334    Process();    Process();
335    
336  //  Float_t dEpianimean[44];    //  Float_t dEpianimean[44];
337   Float_t dEpianiloop[44];    Float_t dEpianiloop[44];
338   Float_t Depth[44];    Float_t Depth[44];
339  //  Int_t tz=(Int_t)qtz;    //  Int_t tz=(Int_t)qtz;
340   Int_t tz1=(Int_t)lpz;    Int_t tz1=(Int_t)lpz;
341  //  Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata    //  Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata
342   Enetrack(&tz1, &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop    Enetrack(&tz1, &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop
343    
344   Float_t sp= spessore[0]*spessore[1];    Float_t sp= spessore[0]*spessore[1];
345   for(Int_t i=0;i<44;i++)Depth[i]=i*sp;    for(Int_t i=0;i<44;i++)Depth[i]=i*sp;
346    //    //
347    gStyle->SetLabelSize(0.04);    gStyle->SetLabelSize(0.04);
348    gStyle->SetNdivisions(510,"XY");    gStyle->SetNdivisions(510,"XY");
349    //    //
350          TString hid = Form("cCaloBragg");          TString hid = Form("cCaloBragg");    
351          TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));    TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
352          if ( tc ){    if ( tc ){
353  //       tc->Clear();      //   tc->Clear();
354          } else {    } else {
355           tc = new TCanvas(hid,hid);      tc = new TCanvas(hid,hid);
356  //       tc->Divide(1,2);      //   tc->Divide(1,2);
357          };    };
358          //    //
359  //      TString thid = Form("hCaloBragg");          //    TString thid = Form("hCaloBragg");      
360  //         TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));    //         TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
361  //      if ( th ) th->Delete();    //    if ( th ) th->Delete();
362  //       th->Clear();    //     th->Clear();
363  //       th->Reset();    //     th->Reset();
364  //      } else {    //    } else {
365  //      th = new TH2F(thid,thid,300,-0.5,300.,1000,0.,150.);    //    th = new TH2F(thid,thid,300,-0.5,300.,1000,0.,150.);
366  //      th->SetMarkerStyle(20);    //    th->SetMarkerStyle(20);
367  //      };    //    };
368          //    //
369          TString thid2 = Form("hCaloBragg2");        tc->cd();
370          TH2F *th2  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid2));    TString thid2 = Form("hCaloBragg2");  
371          if ( th2 ) th2->Delete();    TH2F *th2  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid2));
372          th2 = new TH2F(thid2,thid2,300,-0.5,300.,1000,0.,150.);    if ( th2 ) th2->Delete();
373          th2->SetMarkerStyle(20);    th2 = new TH2F(thid2,thid2,300,-0.5,300.,1000,0.,150.);
374          th2->SetMarkerColor(kRed);    th2->SetMarkerStyle(20);
375          //    th2->SetMarkerColor(kRed);
376          TString thid3 = Form("hCaloBragg3");        //
377          TH2F *th3  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid3));    TString thid3 = Form("hCaloBragg3");  
378          if ( th3 ) th3->Delete();    TH2F *th3  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid3));
379          th3 = new TH2F(thid3,thid3,300,-0.5,300.,1000,0.,150.);    if ( th3 ) th3->Delete();
380          th3->SetMarkerStyle(20);    th3 = new TH2F(thid3,thid3,300,-0.5,300.,1000,0.,150.);
381          th3->SetMarkerColor(kBlue);    th3->SetMarkerStyle(20);
382      th3->SetMarkerColor(kBlue);
383    
384          tc->cd(1);  
385          //    //  tc->cd(1);
386  //      for(Int_t i=0;i<=estremi[1][0];i++)th->Fill(Depth[i],dEpianimean[i]);    //
387          for(Int_t i=0;i<=estremi[1][0];i++)th2->Fill(Depth[i],calorimetro[i][1]*MIP);    //    for(Int_t i=0;i<=estremi[1][0];i++)th->Fill(Depth[i],dEpianimean[i]);
388  //      th->Draw();    for(Int_t i=0;i<=estremi[1][0];i++)th2->Fill(Depth[i],calorimetro[i][1]*MIP);
389          th2->Draw("same");    //    th->Draw();
390      th2->Draw("same");
391          tc->cd(2);  
392          //    //  tc->cd(2);
393          for(Int_t i=0;i<=estremi[1][0];i++)th3->Fill(Depth[i],dEpianiloop[i]);    tc->cd();
394          th3->Draw();    //
395          th2->Draw("same");    for(Int_t i=0;i<=estremi[1][0];i++)th3->Fill(Depth[i],dEpianiloop[i]);
396      th3->Draw();
397      th2->Draw("same");
398    
399          tc->Modified();    tc->Modified();
400          tc->Update();    tc->Update();
401    
402    //    //
403    gStyle->SetLabelSize(0);    gStyle->SetLabelSize(0);
# Line 364  void CaloBragg::LoadParam(){ Line 444  void CaloBragg::LoadParam(){
444    elem[31] = 72.61; //Ge 32    elem[31] = 72.61; //Ge 32
445    
446    
447  //parametri calorimetro    //parametri calorimetro
448    NPLA = 22;    NPLA = 22;
449    NCHA = 96;    NCHA = 96;
450    nView = 2;    nView = 2;
# Line 430  void CaloBragg::conversione(){ Line 510  void CaloBragg::conversione(){
510    WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad    WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad
511    //gcm2W = WCross/10. * rhoW;    //gcm2W = WCross/10. * rhoW;
512    
513   a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664);  //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664);    a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664);  //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664);
514    
515    //       (g/cm2W)/(g/cm2Si)    //       (g/cm2W)/(g/cm2Si)
516    spessore[1] =  a;    spessore[1] =  a;
# Line 469  void CaloBragg::BetheBloch(Float_t *x, F Line 549  void CaloBragg::BetheBloch(Float_t *x, F
549    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));
550        
551    lg = 2.* Me * SQ(eta) * Wmax / SQ(ISi);    lg = 2.* Me * SQ(eta) * Wmax / SQ(ISi);
552                                                                                                 //  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;
553    C=(0.42237*pow(eta,-2.) + 0.0304*pow(eta,-4.) - 0.00038*pow(eta,-6.))*pow(10.,-6.)* pow(ISi,2.) +    C=(0.42237*pow(eta,-2.) + 0.0304*pow(eta,-4.) - 0.00038*pow(eta,-6.))*pow(10.,-6.)* pow(ISi,2.) +
554      (3.858*pow(eta,-2.) - 0.1668*pow(eta,-4.) + 0.00158*pow(eta,-6.))*pow(10.,-9.)*pow(ISi,3.);      (3.858*pow(eta,-2.) - 0.1668*pow(eta,-4.) + 0.00158*pow(eta,-6.))*pow(10.,-9.)*pow(ISi,3.);
555        
556    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);      
557    
558     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.);
559    
560    *out =Energia;//out    *out =Energia;//out
561    
# Line 619  void CaloBragg::chiquadro(Float_t dE[], Line 699  void CaloBragg::chiquadro(Float_t dE[],
699        };        };
700        
701        //meno peso ai piani con rilasci maggiori di 1000 MIP        //meno peso ai piani con rilasci maggiori di 1000 MIP
702          if(calorimetro[ipla][1] > 1000) wi=0.5;        //      if(calorimetro[ipla][1] > 1000) wi=0.5;
703          if(calorimetro[ipla][1] > 1200.) wi=0.5;
704                            
705        Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));        Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));
706    
# Line 636  void CaloBragg::chiquadro(Float_t dE[], Line 717  void CaloBragg::chiquadro(Float_t dE[],
717                        
718      //da Emi      //da Emi
719      if(badplane > 2){            if(badplane > 2){      
720          //      printf(" AAAAAAAAAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG\n");
721        out[1] =79.;        out[1] =79.;
722        break;        break;
723      };      };
# Line 649  void CaloBragg::chiquadro(Float_t dE[], Line 731  void CaloBragg::chiquadro(Float_t dE[],
731    
732    
733    
734  void CaloBragg::loopze( Float_t step, Float_t E0,Float_t Zstart, Float_t Zlimite){  void CaloBragg::loopze( Float_t step, Float_t E0,Float_t Zstart, Float_t Zlimite, Int_t nostep = 1000){
735  //    //
736    //loop su z ed energie per trovare miglior z (ed energia)    //loop su z ed energie per trovare miglior z (ed energia)
737    //in:  nloop     => energia massima da provare (nloop x E0)    //in:  nloop     => energia massima da provare (nloop x E0)
738    //     E0        => energia iniziale (intergale)            //     E0        => energia iniziale (intergale)        
# Line 668  void CaloBragg::loopze( Float_t step, Fl Line 750  void CaloBragg::loopze( Float_t step, Fl
750    
751    Float_t Massa = 0.;    Float_t Massa = 0.;
752    
753    Float_t Stepint =(step)/1000.;//passo per il calcolo di energia    Float_t Stepint =(step)/(Float_t)nostep;//passo per il calcolo di energia
754    
755    Float_t energia =0.;//energia del loop    Float_t energia =0.;//energia del loop
756    
# Line 683  void CaloBragg::loopze( Float_t step, Fl Line 765  void CaloBragg::loopze( Float_t step, Fl
765    Float_t matrixchi2[colmax][rowmax][3];    Float_t matrixchi2[colmax][rowmax][3];
766    memset(matrixchi2, 0, colmax*rowmax*3*sizeof(Float_t));    memset(matrixchi2, 0, colmax*rowmax*3*sizeof(Float_t));
767    
768      Int_t imin = 1-nostep/2;
769      Int_t imax = nostep/2;
770    
771    //loop elementi      //loop elementi  
772    for(Int_t inucl=(int)(Zstart); inucl<max; inucl++){    for(Int_t inucl=(int)(Zstart); inucl<max; inucl++){
# Line 692  void CaloBragg::loopze( Float_t step, Fl Line 776  void CaloBragg::loopze( Float_t step, Fl
776      Massa = elem[inucl-1]*MassP;      Massa = elem[inucl-1]*MassP;
777            
778      //loop energia      //loop energia
779      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 <=??      Int_t iene2 = 0;
780            //    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 <=??
781        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 <=??
782          
783          iene2++;
784        energia=  Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica..        energia=  Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica..
785    
786        Enetrack(&Z, &energia, &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani        Enetrack(&Z, &energia, &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani
# Line 702  void CaloBragg::loopze( Float_t step, Fl Line 789  void CaloBragg::loopze( Float_t step, Fl
789        chiquadro(dEplan,chi2);        chiquadro(dEplan,chi2);
790                
791        if( (chi2[1] != 79.) ){//salto quelli che frammentano        if( (chi2[1] != 79.) ){//salto quelli che frammentano
792          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
793          matrixchi2[inucl][iene][1]=energia;//energia per questo chi2          matrixchi2[inucl][iene2][1]=energia;//energia per questo chi2
794          matrixchi2[inucl][iene][2]=chi2[2];//piani saltati nel chi2          matrixchi2[inucl][iene2][2]=chi2[2];//piani saltati nel chi2
795        } else {          } else {  
796          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
797          matrixchi2[inucl][iene][1]=1000;//energia per questo chi2          matrixchi2[inucl][iene2][1]=1000.;//energia per questo chi2
798          matrixchi2[inucl][iene][2]=1000;//piani saltati nel chi2          matrixchi2[inucl][iene2][2]=1000.;//piani saltati nel chi2
799          break;          break;
800        }        }
801      }//fine loop energia      }//fine loop energia
802    
803    
804   };//fine loop z    };//fine loop z
805    
806    
807    //Emi    //Emi
808    for (Int_t nu=(int)(Zstart); nu<max; nu++){    for (Int_t nu=(int)(Zstart); nu<max; nu++){
809      for (Int_t en=0; en<1000; en++){      for (Int_t en=0; en<nostep; en++){
810        if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] >0.)){        if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] >0.)){
811          bestchi2[0]= matrixchi2[nu][en][0];// chi2          bestchi2[0]= matrixchi2[nu][en][0];// chi2
812          bestchi2[1]= (Float_t)nu; // z          bestchi2[1]= (Float_t)nu; // z
# Line 845  void CaloBragg::Zdaloop(){ Line 932  void CaloBragg::Zdaloop(){
932    
933    
934    /*z se particella fosse al minimo*/  //energia1piano/mip corretta    /*z se particella fosse al minimo*/  //energia1piano/mip corretta
935    Float_t zmax = round(sqrt(estremi[0][1]/spessore[2]));    //  Float_t zmax = round(sqrt(estremi[0][1]/spessore[2]));
936    if(zmax<31)zmax=zmax+1;    //  if(zmax<31)zmax=zmax+1;
937        
938    /*calcolo Z ed E con loop sui vari elementi ed energie*/    /*calcolo Z ed E con loop sui vari elementi ed energie*/
939    
940    Float_t zmin=1.;    Float_t zmin=1.;
941      Float_t zmax=32.;
942    Float_t bestchitemp[4] = {0,0,0,0};    Float_t bestchitemp[4] = {0,0,0,0};
943    
944    bestchi2[0]=10000.;    bestchi2[0]=10000.;
# Line 860  void CaloBragg::Zdaloop(){ Line 948  void CaloBragg::Zdaloop(){
948    Float_t zero=0.;    Float_t zero=0.;
949    //------------primo loop   ----------------------    //------------primo loop   ----------------------
950    //     energia   ezero, zstart  zstop    //     energia   ezero, zstart  zstop
951    loopze(Integrale,zero,zmin,zmax);    //  loopze(Integrale,zero,zmin,zmax);
952      loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,50);
953      //  loopze(Integrale*2.,Integrale/100.,zmin,zmax);
954      //  printf(" Integrale %f , outene %f \n",Integrale,bestchi2[2]);
955        
956    //------------secondo loop  ----------------------    //------------secondo loop  ----------------------
957    for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i];    for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i];
# Line 871  void CaloBragg::Zdaloop(){ Line 961  void CaloBragg::Zdaloop(){
961    bestchi2[3] = 0.;//riazzero    bestchi2[3] = 0.;//riazzero
962    
963    Float_t step = bestchitemp[2];//    Float_t step = bestchitemp[2];//
964    zero=0;  // qualsiasi altro valore peggiora le cose    zero=0.;  // qualsiasi altro valore peggiora le cose
965    zmin=zmax=bestchitemp[1];    //  zmin=zmax=bestchitemp[1];
966    loopze(step,zero,zmin,zmax); //    zmin=bestchitemp[1]-1;
967      zmax=bestchitemp[1]+1;
968      //  loopze(step,zero,zmin,zmax); //
969      loopze(step,step/2.,zmin,zmax,200); //
970    
971        
972    //chi2,z,Etot,Pskip    //chi2,z,Etot,Pskip
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.23