/[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.13 by mocchiut, Thu Oct 6 09:27:39 2011 UTC revision 1.15 by biancoa, Wed Nov 23 14:33:32 2011 UTC
# Line 92  void CaloBragg::Process(){ Line 92  void CaloBragg::Process(){
92  };  };
93    
94    
95  void CaloBragg::CleanPlanes(Float_t epiano[22][2]){  void CaloBragg::CleanPlanes(Float_t epiano[22][2], Bool_t zpiano[22][2]){
96    //  return;    //  return;
97    Int_t hitplanes = 0;    Int_t hitplanes = 0;
98      Float_t f5 = 0.;
99    for (Int_t i = 0; i<22; i++){    for (Int_t i = 0; i<22; i++){
100      for (Int_t j = 1; j>=0; j--){      for (Int_t j = 1; j>=0; j--){
101        if ( epiano[i][j] > 0.7 ) hitplanes++;        zpiano[i][j] = false;
102          if ( epiano[i][j] > 0.7 ){
103            if ( hitplanes < 100 ) f5 += epiano[i][j];
104            hitplanes++;
105          };
106      };      };
107    };    };
108      Int_t atl5 = TMath::Min(hitplanes,100);
109      atl5 = TMath::Max(atl5,1);
110    Float_t lowlim = 0.85;    Float_t lowlim = 0.85;
111      //Float_t lowlim = 1.;
112    Float_t dedxone = 0.;    Float_t dedxone = 0.;
113    Float_t step1 = 0.8*L2->GetCaloLevel2()->qtot/(Float_t)hitplanes;    //  Float_t step1 = 0.8*L2->GetCaloLevel2()->qtot/(Float_t)hitplanes;
114    while ( dedxone < step1 ){    Float_t step1 = 0.8*f5/atl5;
115      //  while ( dedxone < step1 ){
116      for (Int_t i = 0; i<22; i++){      for (Int_t i = 0; i<22; i++){
117        for (Int_t j = 1; j>=0; j--){        for (Int_t j = 1; j>=0; j--){
118            if (debug) printf("Acleanplanes: i %i j %i step1 %f dedxone %f epiano[i][j] %f \n",i,j,step1,dedxone,epiano[i][j]);
119          if ( epiano[i][j] >= step1 && dedxone < 0.7 ) dedxone = epiano[i][j];          if ( epiano[i][j] >= step1 && dedxone < 0.7 ) dedxone = epiano[i][j];
120            if ( dedxone >= step1 ) break; // new
121        };        };
122          if ( dedxone >= step1 ) break; // new
123      };          };    
124    }      //  }
125    if ( dedxone < 0.7 ){    if ( dedxone < 0.7 ){   // here we could have instead while dedxone == 0. ... perhaps better...
126      for (Int_t i = 0; i<22; i++){      for (Int_t i = 0; i<22; i++){
127        for (Int_t j = 1; j>=0; j--){        for (Int_t j = 1; j>=0; j--){
128            if (debug) printf("Bcleanplanes dedxone < 0.7: i %i j %i step1 %f dedxone %f epiano[i][j] %f \n",i,j,step1,dedxone,epiano[i][j]);
129          if ( epiano[i][j] > 0. && dedxone < 0.7 ) dedxone = epiano[i][j];          if ( epiano[i][j] > 0. && dedxone < 0.7 ) dedxone = epiano[i][j];
130            if ( dedxone >= 0.7 ) break; // new
131        };        };
132          if ( dedxone >= 0.7 ) break; // new
133      };          };    
134    }    }
135    //    //
# Line 125  void CaloBragg::CleanPlanes(Float_t epia Line 140  void CaloBragg::CleanPlanes(Float_t epia
140    for (Int_t i = 0; i<22; i++){    for (Int_t i = 0; i<22; i++){
141      for (Int_t j = 1; j>=0; j--){      for (Int_t j = 1; j>=0; j--){
142        if ( epiano[i][j] < dedxone*lowlim ){        if ( epiano[i][j] < dedxone*lowlim ){
143          //        printf(" %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius);          if ( debug ) printf("Ccleanplanes: %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius);
144          epiano[i][j] = 0.;          //        epiano[i][j] = 0.;
145            zpiano[i][j] = true;
146            if ( epiano[i][j] < dedxone*0.05 ) epiano[i][j] = 0.;
147        } else {        } else {
148          //x        printf(" %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius);          if ( debug ) printf("Dcleanplanes else: %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius);
149          nulliferus = 0;          nulliferus = 0;
150          revulsera = true;          revulsera = true;
151        };        };
152        if ( epiano[i][j] < 0.7 && revulsera ) nulliferus++;        //      if ( epiano[i][j] < 0.7 && revulsera ) nulliferus++;
153          if ( (zpiano[i][j] || epiano[i][j] < 0.7 ) && revulsera ) nulliferus++;
154        if ( nulliferus > 10 ) nullius = true;        if ( nulliferus > 10 ) nullius = true;
155        if ( nullius ) epiano[i][j] = 0.;        //      if ( nullius ) epiano[i][j] = 0.;
156          if ( nullius ) zpiano[i][j] = true;
157      };      };
158    };    };
159    
# Line 191  void CaloBragg::Process(Int_t ntr){ Line 210  void CaloBragg::Process(Int_t ntr){
210      //      //
211    };    };
212    //    //
213    this->CleanPlanes(*&epiano);    Bool_t zpiano[22][2];
214      this->CleanPlanes(*&epiano, *&zpiano);
215    //    //
216    PamTrack *ptrack = 0;    PamTrack *ptrack = 0;
217    CaloTrkVar *track = 0;    CaloTrkVar *track = 0;
# Line 220  void CaloBragg::Process(Int_t ntr){ Line 240  void CaloBragg::Process(Int_t ntr){
240    if(L2->GetCaloLevel2()->npcfit[0]==0 && L2->GetCaloLevel2()->npcfit[1]==0 && L2->GetCaloLevel2()->npcfit[2]==0 && L2->GetCaloLevel2()->npcfit[3]==0) return;// controllo sulla traccia nel calorimetro    if(L2->GetCaloLevel2()->npcfit[0]==0 && L2->GetCaloLevel2()->npcfit[1]==0 && L2->GetCaloLevel2()->npcfit[2]==0 && L2->GetCaloLevel2()->npcfit[3]==0) return;// controllo sulla traccia nel calorimetro
241    
242    //    //
243      Bool_t zcalo[44];
244    for(Int_t p=0; p<22; p++){    for(Int_t p=0; p<22; p++){
245      for(Int_t v=0; v<2; v++){      for(Int_t v=0; v<2; v++){
246        /*per usare traccia non del calo camboare cibar*/        /*per usare traccia non del calo camboare cibar*/
247        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      
248        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
249          zcalo[(2*p)+1-v] = zpiano[p][v];
250          if ( debug ) printf(" idx %i %f %i \n",(2*p)+1-v,epiano[p][v], zpiano[p][v]);
251      };      };
252    };    };
253    
# Line 288  void CaloBragg::Process(Int_t ntr){ Line 311  void CaloBragg::Process(Int_t ntr){
311    //    //
312    
313    Float_t lastok = 0.;    Float_t lastok = 0.;
314      //  if ( false ){
315    //  Bool_t goback = false;    //  Bool_t goback = false;
316    for ( int o = 0; o < estremi[1][0]; o++ ){    for ( int o = 0; o < estremi[1][0]; o++ ){
317      //      //
318      if ( calorimetro[o][1] > 0.7 ) lastok = calorimetro[o][1];      if (debug)    printf(" goforth1: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);
319      if ( calorimetro[o][1] < 0.7 && lastok > 0. ) calorimetro[o][1] = lastok;          if ( calorimetro[o][1] > 0.7 && !zcalo[o] ) lastok = calorimetro[o][1];
320        if ( (zcalo[o] || calorimetro[o][1] < 0.7) && lastok > 0. ){
321          if ( fabs(calorimetro[o][1]-lastok)/calorimetro[o][1] > 0.5 ) {
322            if (debug) printf(" goforthchange %f %f \n",calorimetro[o][1],lastok);
323            calorimetro[o][1] = lastok;    
324            if (debug)        printf(" goforthchang+ %f %f \n",calorimetro[o][1],lastok);
325          }
326        }
327        if (debug)    printf(" goforth2: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);
328      //    if ( calorimetro[o][1] < 0.7 ) goback = true;      //    if ( calorimetro[o][1] < 0.7 ) goback = true;
329      //      //
330    };    };
# Line 300  void CaloBragg::Process(Int_t ntr){ Line 332  void CaloBragg::Process(Int_t ntr){
332    //  if ( goback ){    //  if ( goback ){
333    for ( int o = estremi[1][0]; o >= 0;  o-- ){    for ( int o = estremi[1][0]; o >= 0;  o-- ){
334      //      //
335      //    printf(" goback1: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);      if (debug)    printf(" goback1: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);
336      if ( o < estremi[1][0] && calorimetro[o][1] > calorimetro[o+1][1]*1.2 && lastok > 0. ) calorimetro[o][1] = lastok;        if ( o < estremi[1][0] && calorimetro[o][1] > calorimetro[o+1][1]*1.2 && lastok > 0. ) calorimetro[o][1] = lastok;  
337      if ( calorimetro[o][1] > 0.7 ) lastok = calorimetro[o][1];      if ( calorimetro[o][1] > 0.7 && !zcalo[o] ) lastok = calorimetro[o][1];
338      if ( calorimetro[o][1] < 0.7 && lastok > 0. ) calorimetro[o][1] = lastok;          if ( (zcalo[o] || calorimetro[o][1] < 0.7) && lastok > 0. ){
339      //    printf(" goback2: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);        if ( fabs(calorimetro[o][1]-lastok)/calorimetro[o][1] > 0.5 ) {
340            if ( debug ) printf(" gobackchange \n");
341            calorimetro[o][1] = lastok;    
342          }
343        }
344        if (debug)    printf(" goback2: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);
345      //      //
346    };    };
347    //  };    //  };
348      //}
349    
350    if ( startZero ) {    if ( startZero ) {
351      estremi[0][0] = 0.;      estremi[0][0] = 0.;
# Line 456  void CaloBragg::Draw(Int_t Z, Int_t isot Line 494  void CaloBragg::Draw(Int_t Z, Int_t isot
494    TString thid2 = Form("hCaloBragg2");      TString thid2 = Form("hCaloBragg2");  
495    TH2F *th2  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid2));    TH2F *th2  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid2));
496    if ( th2 ) th2->Delete();    if ( th2 ) th2->Delete();
497    th2 = new TH2F(thid2,thid2,300,-0.5,300.,1000,0.,150.);    th2 = new TH2F(thid2,thid2,300,-0.5,300.,1000,0.,25.); //150
498    th2->SetMarkerStyle(20);    th2->SetMarkerStyle(20);
499    th2->SetMarkerColor(kRed);    th2->SetMarkerColor(kRed);
500    //    //
501    TString thid3 = Form("hCaloBragg3");      TString thid3 = Form("hCaloBragg3");  
502    TH2F *th3  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid3));    TH2F *th3  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid3));
503    if ( th3 ) th3->Delete();    if ( th3 ) th3->Delete();
504    th3 = new TH2F(thid3,thid3,300,-0.5,300.,1000,0.,150.);    th3 = new TH2F(thid3,thid3,300,-0.5,300.,1000,0.,25.);//150.
505    th3->SetMarkerStyle(20);    th3->SetMarkerStyle(20);
506    th3->SetMarkerColor(kBlue);    th3->SetMarkerColor(kBlue);
507    
# Line 527  void CaloBragg::LoadParam(){ Line 565  void CaloBragg::LoadParam(){
565    
566    elem[3][0] = 9.012182;  //Be 4    elem[3][0] = 9.012182;  //Be 4
567    elem[3][1] = 10.01353;  //10Be (Isotope) (most stable)    elem[3][1] = 10.01353;  //10Be (Isotope) (most stable)
568    elem[3][2] = -1.;    elem[3][2] = 7.01693;   //9Be no EC in space?
569    elem[3][3] = -1.;    elem[3][3] = -1.;
570    elem[3][4] = -1.;    elem[3][4] = -1.;
571    elem[3][5] = -1.;    elem[3][5] = -1.;
572    elem[3][6] = -1.;    elem[3][6] = -1.;
573    
574    elem[4][0] = 11.0093;   //B 5    elem[4][0] = 11.00930;   //B 5
575    elem[4][1] = 10.01294;  //10B  (Isotope)    elem[4][1] = 10.01294;  //10B  (Isotope)
576    elem[4][2] = -1.;    elem[4][2] = -1.;
577    elem[4][3] = -1.;    elem[4][3] = -1.;
# Line 958  void CaloBragg::Enetrack(Int_t* Z, Int_t Line 996  void CaloBragg::Enetrack(Int_t* Z, Int_t
996      }else{      }else{
997        out[ipla] = dE; //MeV        out[ipla] = dE; //MeV
998        Ezero = Ezero - dE;//energia residua        Ezero = Ezero - dE;//energia residua
999        if ( debug ) printf(" zompa %i out %f dE %f ezero %f \n",ipla,out[ipla],dE,Ezero);        //      if ( debug ) printf(" zompa %i out %f dE %f ezero %f \n",ipla,out[ipla],dE,Ezero);
1000      };      };
1001      //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
1002      if(ipla%2 == 0){      if(ipla%2 == 0){
# Line 1131  if(debug)  printf("loopze inizio max  %d Line 1169  if(debug)  printf("loopze inizio max  %d
1169    
1170    //loop isotopi    //loop isotopi
1171    while ( elem[inucl-1][isotope] > 0. ){      while ( elem[inucl-1][isotope] > 0. ){  
1172    
1173              if( fiso != -1 ){
1174                      isotope=fiso;
1175                      if(debug) printf("In Loopze - Isotope N %d",isotope);
1176              }
1177      Massa = elem[inucl-1][isotope]*MassP;      Massa = elem[inucl-1][isotope]*MassP;
1178    
1179      //loop energia      //loop energia
# Line 1158  if(debug)  printf("loopze inizio max  %d Line 1201  if(debug)  printf("loopze inizio max  %d
1201                          matrixchi2[inucl][isotope][iene2][2]=chi2[2];//piani saltati nel chi2                          matrixchi2[inucl][isotope][iene2][2]=chi2[2];//piani saltati nel chi2
1202              if( fene > 0. ) break;              if( fene > 0. ) break;
1203           } else {                 } else {      
1204                          matrixchi2[inucl][isotope][iene2][0]=1000.;//valore chi2 per questo z a questa energia                          matrixchi2[inucl][isotope][iene2][0]=numeric_limits<Float_t>::max();//valore chi2 per questo z a questa energia
1205                          matrixchi2[inucl][isotope][iene2][1]=1000.;//energia per questo chi2                          matrixchi2[inucl][isotope][iene2][1]=numeric_limits<Float_t>::max();//energia per questo chi2
1206                          matrixchi2[inucl][isotope][iene2][2]=1000.;//piani saltati nel chi2                          matrixchi2[inucl][isotope][iene2][2]=numeric_limits<Float_t>::max();//piani saltati nel chi2
1207                  break;                  break;
1208        }        }
1209    
1210      }//fine loop energia      }//fine loop energia
1211    
1212      if( fiso != -1 ){
1213              if(debug) printf("exited form isotopes loop");
1214              break;
1215      }
1216    
1217    isotope++; //incremento il contatore isotopi    isotope++; //incremento il contatore isotopi
1218    }//fine loop isotopi    }//fine loop isotopi
1219    isotope=0; //riazzero il contatore isotopi    isotope=0; //riazzero il contatore isotopi
# Line 1175  if(debug)  printf("loopze inizio max  %d Line 1223  if(debug)  printf("loopze inizio max  %d
1223    isotope=0;//non dovrebbe servire    isotope=0;//non dovrebbe servire
1224    
1225    //Emi    //Emi
1226    for (Int_t nu=zmin; nu<max; nu++){   for (Int_t nu=zmin; nu<max; nu++){
1227      while(elem[nu-1][isotope]> 0.){     if( fiso != -1 ){
1228              isotope=fiso;
1229              if(debug) printf("In Loopze EMI - Isotope N %d",isotope);
1230       }
1231         while(elem[nu-1][isotope]> 0.){
1232          for (Int_t en=0; en<nostep; en++){          for (Int_t en=0; en<nostep; en++){
1233          if((matrixchi2[nu][isotope][en][0]<bestchi2[0]) && (matrixchi2[nu][isotope][en][0] >0.)){          if((matrixchi2[nu][isotope][en][0]<bestchi2[0]) && (matrixchi2[nu][isotope][en][0] >0.)){
1234            bestchi2[0]= matrixchi2[nu][isotope][en][0];// chi2            bestchi2[0]= matrixchi2[nu][isotope][en][0];// chi2
# Line 1185  if(debug)  printf("loopze inizio max  %d Line 1237  if(debug)  printf("loopze inizio max  %d
1237            bestchi2[3]= matrixchi2[nu][isotope][en][2];// totale piani saltati            bestchi2[3]= matrixchi2[nu][isotope][en][2];// totale piani saltati
1238        bestchi2[4]= (Float_t)isotope; //isotopo        bestchi2[4]= (Float_t)isotope; //isotopo
1239          }          }
1240        }       }
1241    
1242      if( fiso != -1 ){
1243        if(debug) printf("exited form isotopes loop");
1244        break;
1245      }
1246    
1247      isotope++;      isotope++;
1248      }      }
1249     isotope=0;     isotope=0;
1250    }   }
1251    
1252  };//endloopze  };//endloopze
1253    
Legend:
Removed from v.1.13  
changed lines
  Added in v.1.15
  ViewVC Help
Powered by ViewVC 1.1.23