/[PAMELA software]/DarthVader/ToFLevel2/src/ToFLevel2.cpp
ViewVC logotype

Diff of /DarthVader/ToFLevel2/src/ToFLevel2.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.23 by pamelats, Thu Dec 4 11:41:08 2008 UTC revision 1.26 by carbone, Fri Nov 20 11:05:21 2009 UTC
# Line 5  Line 5 
5   * WM dec 2008: Description of "GetdEdx" changed   * WM dec 2008: Description of "GetdEdx" changed
6   * WM dec 2008: "GetdEdxPaddle" modified: Now includes saturation limit   * WM dec 2008: "GetdEdxPaddle" modified: Now includes saturation limit
7   *              PMTs higher than the saturation limit are not used for dEdx   *              PMTs higher than the saturation limit are not used for dEdx
8     * WM apr 2009: bug found by Nicola in method "GetPaddlePlane"
9   */   */
10    
 #include <TObject.h>  
11  #include <ToFLevel2.h>  #include <ToFLevel2.h>
 #include <iostream>  
12  using namespace std;  using namespace std;
13  ClassImp(ToFPMT);  ClassImp(ToFPMT);
14    ClassImp(ToFdEdx);
15    ClassImp(ToFGeom);
16  ClassImp(ToFTrkVar);  ClassImp(ToFTrkVar);
17  ClassImp(ToFLevel2);  ClassImp(ToFLevel2);
18    
# Line 736  void ToFLevel2::GetPMTPaddle(Int_t pmt_i Line 737  void ToFLevel2::GetPMTPaddle(Int_t pmt_i
737  // gf Apr 07  // gf Apr 07
738    
739  void ToFLevel2::GetPaddlePMT(Int_t paddle, Int_t &pmtleft, Int_t &pmtright){  void ToFLevel2::GetPaddlePMT(Int_t paddle, Int_t &pmtleft, Int_t &pmtright){
740      pmtleft=paddle*2;
741    if(paddle==0){    pmtright= pmtleft+1;  
     pmtleft=0;  
     pmtright=1;  
   }  
   
   if(paddle==1){  
     pmtleft=2;  
     pmtright=3;  
   }  
   
   if(paddle==2){  
     pmtleft=4;  
     pmtright=5;  
   }  
   
   if(paddle==3){  
     pmtleft=6;  
     pmtright=7;  
   }  
   
   if(paddle==4){  
     pmtleft=8;  
     pmtright=9;  
   }  
   
   if(paddle==5){  
     pmtleft=10;  
     pmtright=11;  
   }  
   
   if(paddle==6){  
     pmtleft=12;  
     pmtright=13;  
   }  
   
   if(paddle==7){  
     pmtleft=14;  
     pmtright=15;  
   }  
   
   if(paddle==8){  
     pmtleft=16;  
     pmtright=17;  
   }  
   
   if(paddle==9){  
     pmtleft=18;  
     pmtright=19;  
   }  
   
   if(paddle==10){  
     pmtleft=20;  
     pmtright=21;  
   }  
   
   if(paddle==11){  
     pmtleft=22;  
     pmtright=23;  
   }  
   
   if(paddle==12){  
     pmtleft=24;  
     pmtright=25;  
   }  
   
   if(paddle==13){  
     pmtleft=26;  
     pmtright=27;  
   }  
   
   if(paddle==14){  
     pmtleft=28;  
     pmtright=29;  
   }  
   
   if(paddle==15){  
     pmtleft=30;  
     pmtright=31;  
   }  
   
   if(paddle==16){  
     pmtleft=32;  
     pmtright=33;  
   }  
   
   if(paddle==17){  
     pmtleft=34;  
     pmtright=35;  
   }  
   
   if(paddle==18){  
     pmtleft=36;  
     pmtright=37;  
   }  
   
   if(paddle==19){  
     pmtleft=38;  
     pmtright=39;  
   }  
   
   if(paddle==20){  
     pmtleft=40;  
     pmtright=41;  
   }  
   
   if(paddle==21){  
     pmtleft=42;  
     pmtright=43;  
   }  
   
   if(paddle==22){  
     pmtleft=44;  
     pmtright=45;  
   }  
   
   if(paddle==23){  
     pmtleft=46;  
     pmtright=47;  
   }  
     
742    return;    return;
743  }  }
744    
# Line 970  void ToFLevel2::GetPaddleGeometry(Int_t Line 852  void ToFLevel2::GetPaddleGeometry(Int_t
852   */   */
853  Int_t ToFLevel2::GetPaddleid(Int_t plane, Int_t paddle)  Int_t ToFLevel2::GetPaddleid(Int_t plane, Int_t paddle)
854  {  {
   
855    Int_t padid=-1;    Int_t padid=-1;
856    Int_t pads11=8;    Int_t pads[6]={8,6,2,2,3,3};
   Int_t pads12=6;  
   Int_t pads21=2;  
   Int_t pads22=2;  
   Int_t pads31=3;  
   //  Int_t pads32=3;  
   
   
   if(plane == 0){  
     padid=paddle;  
   }  
857    
858    if(plane == 1){    int somma=0;
859      padid=pads11+paddle;    int np=plane;
860    }    for(Int_t j=0; j<np; j++){
861        somma+=pads[j];
   if(plane == 2){  
     padid=pads11+pads12+paddle;  
862    }    }
863      padid=paddle+somma;
   if(plane == 3){  
     padid=pads11+pads12+pads21+paddle;  
   }  
   
   if(plane == 4){  
     padid=pads11+pads12+pads21+pads22+paddle;  
   }  
   
   if(plane == 5){  
     padid=pads11+pads12+pads21+pads22+pads31+paddle;  
   }  
   
864    return padid;    return padid;
865    
866  }  }
# Line 1034  void ToFLevel2::GetPaddlePlane(Int_t pad Line 891  void ToFLevel2::GetPaddlePlane(Int_t pad
891      return;      return;
892    }    }
893    
894    if(7<pad<14){    if((7<pad)&&(pad<14)){
895      plane=1;      plane=1;
896      paddle=pad-pads11;      paddle=pad-pads11;
897      return;      return;
898    }    }
899        
900    if(13<pad<16){    if((13<pad)&&(pad<16)){
901      plane=2;      plane=2;
902      paddle=pad-pads11-pads12;      paddle=pad-pads11-pads12;
903      return;      return;
904    }    }
905    
906    if(15<pad<18){    if((15<pad)&&(pad<18)){
907      plane=3;      plane=3;
908      paddle=pad-pads11-pads12-pads21;      paddle=pad-pads11-pads12-pads21;
909      return;      return;
910    }    }
911    
912    if(17<pad<21){    if((17<pad)&&(pad<21)){
913      plane=4;      plane=4;
914      paddle=pad-pads11-pads12-pads21-pads22;      paddle=pad-pads11-pads12-pads21-pads22;
915      return;      return;
916    }    }
917    
918    if(20<pad<24){    if((20<pad)&&(pad<24)){
919      plane=5;      plane=5;
920      paddle=pad-pads11-pads12-pads21-pads22-pads31;      paddle=pad-pads11-pads12-pads21-pads22-pads31;
921      return;      return;
# Line 1283  void ToFLevel2::GetLevel2Struct(cToFLeve Line 1140  void ToFLevel2::GetLevel2Struct(cToFLeve
1140        }        }
1141    } //ELENA    } //ELENA
1142  }  }
1143    
1144    
1145    //
1146    // Reprocessing tool // Emiliano 08/04/07
1147    //
1148    Int_t ToFLevel2::Process(TrkLevel2 *trk, TrigLevel2 *trg, GL_RUN *run, OrbitalInfo *orb, Bool_t force){
1149      //
1150      // Copiare qui qualcosa di simile a calonuclei per evitare di riprocessare sempre tutto
1151      //
1152    
1153    
1154    
1155    
1156      //
1157      // structures to communicate with F77
1158      //
1159      extern struct ToFInput  tofinput_;
1160      extern struct ToFOutput tofoutput_;
1161      //
1162      // DB connection
1163      //
1164      TString host;
1165      TString user;
1166      TString psw;
1167      const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
1168      const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
1169      const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
1170      if ( !pamdbhost ) pamdbhost = "";
1171      if ( !pamdbuser ) pamdbuser = "";
1172      if ( !pamdbpsw ) pamdbpsw = "";
1173      if ( strcmp(pamdbhost,"") ) host = pamdbhost;
1174      if ( strcmp(pamdbuser,"") ) user = pamdbuser;
1175      if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
1176      //
1177      //
1178      TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
1179      if ( !dbc->IsConnected() ) return 1;
1180      stringstream myquery;
1181      myquery.str("");
1182      myquery << "SET time_zone='+0:00'";
1183      dbc->Query(myquery.str().c_str());
1184      GL_PARAM *glparam = new GL_PARAM();
1185      glparam->Query_GL_PARAM(1,1,dbc); // parameters stored in DB in GL_PRAM table
1186      trk->LoadField(glparam->PATH+glparam->NAME);
1187      //
1188      Bool_t defcal = true;
1189      Int_t error=glparam->Query_GL_PARAM(run->RUNHEADER_TIME,201,dbc); // parameters stored in DB in GL_PRAM table
1190      if ( error<0 ) {
1191        return(1);
1192      };
1193      printf(" Reading ToF parameter file: %s \n",(glparam->PATH+glparam->NAME).Data());
1194      if ( (UInt_t)glparam->TO_TIME != (UInt_t)4294967295UL ) defcal = false;
1195      //
1196      Int_t nlen = (Int_t)(glparam->PATH+glparam->NAME).Length();
1197      rdtofcal((char *)(glparam->PATH+glparam->NAME).Data(),&nlen);
1198      //
1199      Int_t adc[4][12];
1200      Int_t tdc[4][12];
1201      Float_t tdcc[4][12];
1202      //
1203      // process tof data
1204      //
1205      for (Int_t hh=0; hh<12;hh++){
1206        for (Int_t kk=0; kk<4;kk++){
1207               adc[kk][hh] = 4095;
1208               tdc[kk][hh] = 4095;
1209               tdcc[kk][hh] = 4095.;
1210               tofinput_.adc[hh][kk] = 4095;
1211               tofinput_.tdc[hh][kk] = 4095;
1212        };
1213      };
1214      Int_t ntrkentry = 0;
1215      Int_t npmtentry = 0;
1216      Int_t gg = 0;
1217      Int_t hh = 0;
1218      Int_t adcf[48];
1219      memset(adcf, 0, 48*sizeof(Int_t));
1220      Int_t tdcf[48];
1221      memset(tdcf, 0, 48*sizeof(Int_t));
1222      for (Int_t pm=0; pm < this->ntrk() ; pm++){
1223         ToFTrkVar *ttf = this->GetToFTrkVar(pm);
1224         for ( Int_t nc=0; nc < ttf->npmttdc; nc++){
1225                if ( (ttf->tdcflag).At(nc) != 0 ) tdcf[(ttf->pmttdc).At(nc)] = 1;
1226         };
1227         for ( Int_t nc=0; nc < ttf->npmtadc; nc++){
1228                if ( (ttf->adcflag).At(nc) != 0 ) adcf[(ttf->pmtadc).At(nc)] = 1;
1229         };
1230      };
1231      //
1232      for (Int_t pm=0; pm < this->npmt() ; pm++){
1233         ToFPMT *pmt = this->GetToFPMT(pm);
1234         this->GetPMTIndex(pmt->pmt_id, gg, hh);
1235         if ( adcf[pmt->pmt_id] == 0 ){
1236                 tofinput_.adc[gg][hh] = (int)pmt->adc;
1237                 adc[hh][gg] = (int)pmt->adc;
1238         };
1239         if ( tdcf[pmt->pmt_id] == 0 ){
1240                 tofinput_.tdc[gg][hh] = (int)pmt->tdc;
1241                 tdc[hh][gg] = (int)pmt->tdc;
1242         };
1243         tdcc[hh][gg] = (float)pmt->tdc_tw;
1244         // Int_t pppid = this->GetPMTid(hh,gg);
1245         //      printf(" pm %i pmt_id %i pppid %i hh %i gg %i tdcc %f tdc %f adc %f \n",pm,pmt->pmt_id,pppid,hh,gg,pmt->tdc_tw,pmt->tdc,pmt->adc);
1246      };
1247      //
1248      Int_t unpackError = this->unpackError;
1249      //
1250      for (Int_t hh=0; hh<5;hh++){
1251         tofinput_.patterntrig[hh]=trg->patterntrig[hh];
1252      };
1253      //
1254      this->Clear();
1255      //
1256          Int_t pmt_id = 0;
1257          ToFPMT *t_pmt = new ToFPMT();
1258          if(!(this->PMT)) this->PMT = new TClonesArray("ToFPMT",12); //ELENA
1259          TClonesArray &tpmt = *this->PMT;
1260          ToFTrkVar *t_tof = new ToFTrkVar();
1261          if(!(this->ToFTrk)) this->ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA
1262          TClonesArray &t = *this->ToFTrk;
1263          //
1264          //
1265          // Here we have calibrated data, ready to be passed to the FORTRAN routine which will extract common and track-related  variables.
1266          //
1267          npmtentry = 0;
1268          //
1269          ntrkentry = 0;
1270          //
1271          // Calculate tracks informations from ToF alone
1272          //
1273          tofl2com();
1274          //
1275          memcpy(this->tof_j_flag,tofoutput_.tof_j_flag,6*sizeof(Int_t));
1276          //
1277          t_tof->trkseqno = -1;
1278          //
1279          // and now we must copy from the output structure to the level2 class:
1280          //
1281          t_tof->npmttdc = 0;
1282          //
1283          for (Int_t hh=0; hh<12;hh++){
1284            for (Int_t kk=0; kk<4;kk++){
1285              if ( tofoutput_.tofmask[hh][kk] != 0 ){
1286                pmt_id = this->GetPMTid(kk,hh);
1287                t_tof->pmttdc.AddAt(pmt_id,t_tof->npmttdc);
1288                t_tof->tdcflag.AddAt(tofoutput_.tdcflagtof[hh][kk],t_tof->npmttdc); // gf: Jan 09/07
1289                t_tof->npmttdc++;
1290              };
1291            };
1292          };
1293          for (Int_t kk=0; kk<13;kk++){
1294            t_tof->beta[kk] = tofoutput_.betatof_a[kk];
1295          }
1296          //
1297          t_tof->npmtadc = 0;
1298          for (Int_t hh=0; hh<12;hh++){
1299            for (Int_t kk=0; kk<4;kk++){
1300              if ( tofoutput_.adctof_c[hh][kk] < 1000 ){
1301                t_tof->dedx.AddAt(tofoutput_.adctof_c[hh][kk],t_tof->npmtadc);
1302                pmt_id = this->GetPMTid(kk,hh);
1303                t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
1304                t_tof->adcflag.AddAt(tofoutput_.adcflagtof[hh][kk],t_tof->npmtadc); // gf: Jan 09/07
1305                t_tof->npmtadc++;
1306              };
1307            };
1308          };
1309          //
1310          memcpy(t_tof->xtofpos,tofoutput_.xtofpos,sizeof(t_tof->xtofpos));
1311          memcpy(t_tof->ytofpos,tofoutput_.ytofpos,sizeof(t_tof->ytofpos));
1312          memcpy(t_tof->xtr_tof,tofoutput_.xtr_tof,sizeof(t_tof->xtr_tof));
1313          memcpy(t_tof->ytr_tof,tofoutput_.ytr_tof,sizeof(t_tof->ytr_tof));
1314          //
1315          new(t[ntrkentry]) ToFTrkVar(*t_tof);
1316          ntrkentry++;
1317          t_tof->Clear();
1318          //
1319          //
1320          //
1321          t_pmt->Clear();
1322          //
1323          for (Int_t hh=0; hh<12;hh++){
1324            for (Int_t kk=0; kk<4;kk++){
1325             // new WM
1326              if ( tofoutput_.tdc_c[hh][kk] < 4095 || adc[kk][hh] < 4095  || tdc[kk][hh] < 4095 ){
1327    //          if ( tdcc[kk][hh] < 4095. || adc[kk][hh] < 4095  || tdc[kk][hh] < 4095 ){
1328                //
1329                t_pmt->pmt_id = this->GetPMTid(kk,hh);
1330                t_pmt->tdc_tw = tofoutput_.tdc_c[hh][kk];
1331                t_pmt->adc = (Float_t)adc[kk][hh];
1332                t_pmt->tdc = (Float_t)tdc[kk][hh];
1333                //
1334                new(tpmt[npmtentry]) ToFPMT(*t_pmt);
1335                npmtentry++;
1336                t_pmt->Clear();
1337              };
1338            };
1339          };
1340          //
1341          // Calculate track-related variables
1342          //
1343          if ( trk->ntrk() > 0 ){
1344            //
1345            // We have at least one track
1346            //
1347            //
1348            // Run over tracks
1349            //
1350            for(Int_t nt=0; nt < trk->ntrk(); nt++){
1351              //
1352              TrkTrack *ptt = trk->GetStoredTrack(nt);
1353              //
1354              // Copy the alpha vector in the input structure
1355              //
1356              for (Int_t e = 0; e < 5 ; e++){
1357                tofinput_.al_pp[e] = ptt->al[e];
1358              };
1359              //
1360              // Get tracker related variables for this track
1361              //
1362              toftrk();
1363              //
1364              // Copy values in the class from the structure (we need to use a temporary class to store variables).
1365              //
1366              t_tof->npmttdc = 0;
1367              for (Int_t hh=0; hh<12;hh++){
1368                for (Int_t kk=0; kk<4;kk++){
1369                  if ( tofoutput_.tofmask[hh][kk] != 0 ){
1370                    pmt_id = this->GetPMTid(kk,hh);
1371                    t_tof->pmttdc.AddAt(pmt_id,t_tof->npmttdc);
1372                    t_tof->tdcflag.AddAt(tofoutput_.tdcflag[hh][kk],t_tof->npmttdc); // gf: Jan 09/07
1373                    t_tof->npmttdc++;
1374                  };
1375                };
1376              };
1377              for (Int_t kk=0; kk<13;kk++){
1378                t_tof->beta[kk] = tofoutput_.beta_a[kk];
1379              };
1380              //
1381              t_tof->npmtadc = 0;
1382              for (Int_t hh=0; hh<12;hh++){
1383                for (Int_t kk=0; kk<4;kk++){
1384                  if ( tofoutput_.adc_c[hh][kk] < 1000 ){
1385                    t_tof->dedx.AddAt(tofoutput_.adc_c[hh][kk],t_tof->npmtadc);
1386                    pmt_id = this->GetPMTid(kk,hh);
1387                    t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
1388                    t_tof->adcflag.AddAt(tofoutput_.adcflag[hh][kk],t_tof->npmtadc); // gf: Jan 09/07
1389                    t_tof->npmtadc++;
1390                  };
1391                };
1392              };
1393              //
1394              memcpy(t_tof->xtofpos,tofoutput_.xtofpos,sizeof(t_tof->xtofpos));
1395              memcpy(t_tof->ytofpos,tofoutput_.ytofpos,sizeof(t_tof->ytofpos));
1396              memcpy(t_tof->xtr_tof,tofoutput_.xtr_tof,sizeof(t_tof->xtr_tof));
1397              memcpy(t_tof->ytr_tof,tofoutput_.ytr_tof,sizeof(t_tof->ytr_tof));
1398              //
1399              // Store the tracker track number in order to be sure to have shyncronized data during analysis
1400              //
1401              t_tof->trkseqno = nt;
1402              //
1403              // create a new object for this event with track-related variables
1404              //
1405              new(t[ntrkentry]) ToFTrkVar(*t_tof);
1406              ntrkentry++;
1407              t_tof->Clear();
1408              //
1409            }; // loop on all the tracks
1410          //
1411          this->unpackError = unpackError;
1412          if ( defcal ){
1413            this->default_calib = 1;
1414          } else {
1415            this->default_calib = 0;
1416          };
1417     };
1418    
1419    
1420    
1421      return(0);
1422    }
1423    
1424    
1425    
1426    
1427    
1428    
1429    
1430    
1431    
1432    
1433    
1434    
1435    
1436    
1437    
1438    
1439    
1440    
1441    
1442    
1443    
1444    
1445    
1446    
1447    ToFdEdx::ToFdEdx()
1448    {
1449      memset(conn,0,12*sizeof(Bool_t));
1450      memset(ts,0,12*sizeof(UInt_t));
1451      memset(te,0,12*sizeof(UInt_t));
1452      Define_PMTsat();
1453      Clear();
1454    }
1455    //------------------------------------------------------------------------
1456    void ToFdEdx::CheckConnectors(UInt_t atime, GL_PARAM *glparam, TSQLServer *dbc)
1457    {
1458      for(int i=0; i<12; i++){
1459        if(atime<=ts[i] || atime>te[i]){
1460          Int_t error=glparam->Query_GL_PARAM(atime,210+i,dbc); // parameters stored in DB in GL_PRAM table
1461          if ( error<0 ) {
1462            conn[i]=false;
1463            ts[i]=0;
1464            te[i]=numeric_limits<UInt_t>::max();
1465          };
1466          if ( !error ){
1467            conn[i]=true;
1468            ts[i]=glparam->FROM_TIME;
1469            te[i]=glparam->TO_TIME;
1470          }
1471          if ( error>0 ){
1472            conn[i]=false;
1473            ts[i]=glparam->TO_TIME;
1474            TSQLResult *pResult;
1475            TSQLRow *row;
1476            TString query= Form("SELECT FROM_TIME FROM GL_PARAM WHERE TYPE=%i AND FROM_TIME>=%i ORDER BY FROM_TIME ASC LIMIT 1;",210+i,atime);
1477            pResult=dbc->Query(query.Data());
1478            if(!pResult->GetRowCount()){
1479              te[i]=numeric_limits<UInt_t>::max();
1480            }else{
1481              row=pResult->Next();
1482              te[i]=(UInt_t)atoll(row->GetField(0));
1483            }
1484          }
1485          //
1486          
1487        }
1488      }
1489    
1490    }
1491    //------------------------------------------------------------------------
1492    void ToFdEdx::Clear(Option_t *option)
1493    {
1494      //
1495      // Set arrays and initialize structure
1496    
1497      eDEDXpmt.Set(48);    eDEDXpmt.Reset(-1);   // Set array size  and reset structure
1498      eZpmt.Set(48);       eZpmt.Reset(-1);
1499      eDEDXpad.Set(24);    eDEDXpad.Reset(-1);
1500      eZpad.Set(24);       eZpad.Reset(-1);
1501      eDEDXlayer.Set(6);   eDEDXlayer.Reset(-1);
1502      eZlayer.Set(6);      eZlayer.Reset(-1);
1503      eDEDXplane.Set(3);   eDEDXplane.Reset(-1);
1504      eZplane.Set(3);      eZplane.Reset(-1);
1505      INFOpmt.Set(48);     INFOpmt.Reset(0);
1506      INFOlayer.Set(6);    INFOlayer.Reset(0);
1507      //
1508    };
1509    
1510    //------------------------------------------------------------------------
1511    void ToFdEdx::Print(Option_t *option)
1512    {
1513      //
1514      printf("========================================================================\n");
1515    
1516    };
1517    
1518    
1519    //------------------------------------------------------------------------
1520    // void ToFdEdx::InitPar(TString parname, TString parfile)
1521    // {
1522    //   // expensive function - call it once/run
1523    
1524    
1525    
1526    //   ReadParAtt(            Form("%s/attenuation.txt"              , pardir) );
1527    //   ReadParPos(            Form("%s/desaturation_position.txt"    , pardir) );
1528    //   ReadParBBneg(          Form("%s/BetheBloch.txt"               , pardir) );
1529    //   ReadParBBpos(          Form("%s/BetheBloch_betagt1.txt"       , pardir) );
1530    //   ReadParDesatBB(        Form("%s/desaturation_beta.txt"        , pardir) );
1531    
1532    // };
1533    
1534    
1535    //------------------------------------------------------------------------
1536    void ToFdEdx::Process(UInt_t atime, Float_t betamean, Float_t *xtr_tof, Float_t *ytr_tof, pamela::tof::TofEvent *tofl0 )
1537    {
1538      // the parameters should be already initialised by InitPar()
1539    
1540    
1541      Clear();
1542    
1543    
1544    
1545      //  Float_t betamean = fabs(trackTRK->GetToFTrack()->beta[12]);
1546    
1547      if(betamean<0.05 || betamean>2){
1548        for(int i=0;i<48;i++)INFOpmt[i]=1;
1549      }
1550    
1551     // define angle:  
1552      double dx   = xtr_tof[1] - xtr_tof[5];
1553      double dy   = ytr_tof[0] - ytr_tof[4];
1554      double dr   = sqrt(dx*dx+dy*dy);
1555      double theta=atan(dr/76.81);
1556    
1557    
1558    
1559      //  TArrayF adc;
1560      Float_t adc[48];
1561    
1562      ToFLevel2 tf;
1563    
1564      for (Int_t gg=0; gg<4;gg++){
1565        for (Int_t hh=0; hh<12;hh++){
1566          //          tofinput_.tdc[hh][gg]=tofEvent->tdc[gg][hh];          
1567          int mm = tf.GetPMTid(gg,hh);        
1568          adc[mm]=tofl0->adc[gg][hh];
1569        };      
1570      };
1571      
1572      
1573      
1574      
1575      for( int ii=0; ii<48; ii++ ) {
1576        if( adc[ii] >= PMTsat[ii]-5 )  continue;
1577        if( adc[ii] <= 0. )            continue;
1578      
1579        double adcpC   = f_adcPC( adc[ii] );    // - adc conversion in pC
1580        double adccorr = adcpC*fabs(cos(theta));
1581    
1582             if(adccorr<=0.)           continue;
1583    
1584    
1585    
1586        //--------------------- TABLE OF PERIODS WITH HV PROBLEMS ----------------------------
1587    
1588        int Aconn=conn[0];    // PMT 0,20,22,24
1589        int Bconn=conn[1];    // PMT 6,12,26,34
1590        int Cconn=conn[2];    // PMT 4,14,28,32
1591        int Dconn=conn[3];    // PMT 2,8,10,30
1592        int Econn=conn[4];    // PMT 42,43,44,47
1593        int Fconn=conn[5];    // PMT 7,19,23,27
1594        int Gconn=conn[6];    // PMT 3,11,25,33
1595        int Hconn=conn[7];    // PMT 1,9,13,21
1596        int Iconn=conn[8];    // PMT 5,29,31,35
1597        int Lconn=conn[9];    // PMT 37,40,45,46
1598        int Mconn=conn[10];    // PMT 15,16,17,18
1599        int Nconn=conn[11];    // PMT 36,38,39,41
1600    
1601        //    int standard=0;
1602        if( false ) cout << Gconn << Iconn << Lconn <<endl;
1603        int S115B_ok=0;
1604        int S115B_break=0;
1605    
1606    //   if(atime>=1153660001 && atime<=1154375000)Dconn=1;
1607    //     else if(atime>=1155850001 && atime<=1156280000){
1608    //       Hconn=1;
1609    //       Nconn=1;
1610    //     }
1611    
1612    //  else if(atime>=1168490001 && atime<=1168940000)Dconn=1;
1613    //     else if(atime>=1168940001 && atime<=1169580000){
1614    //       Fconn=1;
1615    //       Mconn=1;
1616    //     }
1617    
1618    //  else if(atime>=1174665001 && atime<=1175000000)Bconn=1;
1619    //     else if(atime>=1176120001 && atime<=1176800000)Hconn=1;
1620    //     else if(atime>=1176800001 && atime<=1178330000)Econn=1;
1621    //     else if(atime>=1178330001 && atime<=1181322000)Hconn=1;
1622    //     else if(atime>=1182100001 && atime<=1183030000)Aconn=1;
1623    //     else if(atime>=1184000001 && atime<=1184570000)Hconn=1;
1624    //     else if(atime>=1185090001 && atime<=1185212000)Dconn=1;
1625    //     else if(atime>=1191100001 && atime<=1191940000)Dconn=1;
1626    //     else if(atime>=1196230001 && atime<=1196280000)Hconn=1;
1627    //     else if(atime>=1206100001 && atime<=1206375600)Cconn=1;
1628    //     else if(atime>=1217989201 && atime<=1218547800)Econn=1;
1629    //     else if(atime>=1225789201 && atime<=1226566800)Econn=1;
1630    //     else if(atime>=1229400901 && atime<=1229700000)Econn=1;
1631    //     else if(atime>=1230318001 && atime<=1230415200)Econn=1;
1632    //     else {
1633    //       standard=1;
1634    //     }
1635        if(atime<1158720000)S115B_ok=1;
1636        else S115B_break=1;
1637    
1638    
1639     //------------------------------------------------------------------------
1640    
1641    //---------------------------------------------------- Z reconstruction
1642    
1643    double adcHe, adcnorm, adclin, dEdx, Zeta;
1644    
1645     adcHe=-2;
1646     adcnorm=-2;
1647     adclin=-2;
1648     dEdx=-2;
1649     Zeta=-2;
1650    
1651    
1652    //  float ZetaH=-2;
1653    //  float dEdxH=-2;
1654    
1655    //  double day = (atime-1150000000)/84600;
1656    
1657        if(Aconn==1 && (ii==0 || ii==20 || ii==22 || ii==24)){
1658           adcHe   = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.675;
1659        }
1660        else if(Bconn==1 && (ii==6 || ii==12 || ii==26 || ii==34)){
1661           adcHe   = (Get_adc_he(ii, xtr_tof, ytr_tof))/2.482;
1662        }
1663        else if(Cconn==1 && (ii==4 || ii==14 || ii==28 || ii==32)){
1664          adcHe   = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.464;
1665        }
1666        else if(Dconn==1 && (ii==2 || ii==8 || ii==10 || ii==30)){
1667           adcHe   = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.995;
1668        }
1669        else if(Econn==1 && (ii==42 || ii==43 || ii==44 || ii==47)){
1670           adcHe   = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.273;
1671        }
1672        else if(Fconn==1 && (ii==7 || ii==19 || ii==23 || ii==27)){
1673           adcHe   = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.565;
1674        }
1675        else if(Mconn==1 && (ii==15 || ii==16 || ii==17 || ii==18)){
1676           adcHe   = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.565;
1677        }
1678        else if(Nconn==1 && (ii==36 || ii==38 || ii==39 || ii==41)){
1679           adcHe   = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.018;
1680        }
1681        else if(Hconn==1 && (ii==1 || ii==13 || ii==21 || (ii==9&&S115B_ok==1))){
1682          adcHe   = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.84;
1683        }
1684        else if(S115B_break==1 && ii==9 && Hconn==0){
1685           adcHe   = f_att5B( ytr_tof[0] );   //N.B.: this function refers to the Carbon!!!
1686        }
1687        else if(S115B_break==1 && ii==9 && Hconn==1){
1688           adcHe   = (f_att5B( ytr_tof[0] ))/1.64;
1689        }
1690        else  adcHe   = Get_adc_he(ii, xtr_tof, ytr_tof);
1691    
1692        if(adcHe<=0)   continue;
1693    
1694        if(ii==9 && S115B_break==1)  adcnorm = f_pos5B(adccorr);
1695        else adcnorm = f_pos( (parPos[ii]), adccorr);
1696    
1697        if(adcnorm<=0) continue;
1698    
1699        if(ii==9 && S115B_break==1)  adclin  = 36.*adcnorm/adcHe;
1700        else  adclin  = 4.*adcnorm/adcHe;
1701    
1702        if(adclin<=0)  continue;
1703    
1704        double dEdxHe=-2;
1705        if(ii==9 && S115B_break==1){
1706          if( betamean <1. ) dEdxHe = f_BB5B( betamean );
1707          else                       dEdxHe = 33;
1708        } else {
1709          if( betamean <1. ) dEdxHe = f_BB( (parBBneg[ii]), betamean );
1710          else                       dEdxHe = parBBpos[ii];
1711        }
1712    
1713        if(dEdxHe<=0)  continue;
1714    
1715        if(ii==9 && S115B_break==1)  dEdx = f_desatBB5B( adclin );
1716        else  dEdx = f_desatBB((parDesatBB[ii]), adclin );
1717    
1718        if(dEdx<=0)    continue;
1719    
1720        if(ii==9 && S115B_break==1)  Zeta = sqrt(36.*(dEdx/dEdxHe));
1721        else  Zeta = sqrt(4.*(dEdx/dEdxHe));
1722    
1723        if(Zeta<=0)    continue;
1724    
1725    //--------------------- TIME DEPENDENCE ----------------------------------
1726    
1727    
1728    //      TArrayF &binx = TDx[ii];
1729    //      TArrayF &biny = TDy[ii];
1730    //      for(int k=0; k<200; k++) {
1731    //        if (day > binx[k]-2.5 && day<=binx[k]+2.5 && biny[k]>0)  {
1732    //       ZetaH=Zeta/biny[k]*6;
1733    //       dEdxH=dEdx/(pow(biny[k],2))*36;
1734    //        }
1735    //      }
1736    
1737    //      if(ZetaH!=-2)eZpmt[ii]=(Float_t)ZetaH;
1738    //      else eZpmt[ii]=(Float_t)Zeta;
1739    
1740    //      if(dEdxH!=-2)eDEDXpmt[ii]=(Float_t)dEdxH;
1741    //      else eDEDXpmt[ii]=(Float_t)dEdx;
1742    
1743    //      printf("%5d %8.2f %8.2f %8.2f  %8.2f %8.2f  %8.2f %5.4f \n",               ii, adcpC,  adccorr, adcHe, dEdxHe, dEdx, Zeta, betamean );
1744    
1745        eZpmt[ii]=(Float_t)Zeta;
1746        eDEDXpmt[ii]=(Float_t)dEdx;
1747    
1748    
1749     }  //end loop on 48 PMT
1750    
1751    //---------------------------------------------------  paddle + layer --------------------
1752    
1753      for(int j=0;j<48;j++){
1754        int k=100;
1755        if(j%2==0 || j==0)k=j/2;
1756        
1757        double zpdl=-1;
1758        
1759        if((j%2==0 || j==0) && eZpmt[j]!=-1 && eZpmt[j+1]!=-1){
1760          zpdl=0.5*(eZpmt[j]+eZpmt[j+1]);
1761        }else if((j%2==0 || j==0) && eZpmt[j]!=-1 && eZpmt[j+1]==-1){
1762          zpdl=eZpmt[j];
1763        }else if((j%2==0 || j==0) && eZpmt[j]==-1 && eZpmt[j+1]!=-1){
1764          zpdl=eZpmt[j+1];
1765        }
1766        
1767        if(j%2==0 || j==0)eZpad[k]= (Float_t)zpdl;
1768        
1769        if((j%2==0 || j==0)&&eZpad[k]!=-1){
1770          if(k>=0&&k<8)eZlayer[0]=eZpad[k];
1771          if(k>=8&&k<14)eZlayer[1]=eZpad[k];
1772          if(k>=14&&k<16)eZlayer[2]=eZpad[k];
1773          if(k>=16&&k<18)eZlayer[3]=eZpad[k];
1774          if(k>=18&&k<21)eZlayer[4]=eZpad[k];
1775          if(k>=21)eZlayer[5]=eZpad[k];
1776        }
1777    
1778        if(eZlayer[0]!=-1&&eZlayer[1]!=-1&&fabs(eZlayer[0]-eZlayer[1])<1.5)eZplane[0]=0.5*(eZlayer[0]+eZlayer[1]);
1779        else if(eZlayer[0]!=-1&&eZlayer[1]==-1)eZplane[0]=eZlayer[0];
1780        else if(eZlayer[1]!=-1&&eZlayer[0]==-1)eZplane[0]=eZlayer[1];
1781    
1782        if(eZlayer[2]!=-1&&eZlayer[3]!=-1&&fabs(eZlayer[2]-eZlayer[3])<1.5)eZplane[1]=0.5*(eZlayer[2]+eZlayer[3]);
1783        else if(eZlayer[2]!=-1&&eZlayer[3]==-1)eZplane[1]=eZlayer[2];
1784        else if(eZlayer[3]!=-1&&eZlayer[2]==-1)eZplane[1]=eZlayer[3];
1785    
1786        if(eZlayer[4]!=-1&&eZlayer[5]!=-1&&fabs(eZlayer[4]-eZlayer[5])<1.5)eZplane[2]=0.5*(eZlayer[4]+eZlayer[5]);
1787        else if(eZlayer[4]!=-1&&eZlayer[5]==-1)eZplane[2]=eZlayer[4];
1788        else if(eZlayer[5]!=-1&&eZlayer[4]==-1)eZplane[2]=eZlayer[5];
1789    
1790      }
1791    
1792      for(int jj=0;jj<48;jj++){
1793        int k=100;
1794        if(jj%2==0 || jj==0)k=jj/2;
1795        
1796        double dedxpdl=-1;
1797        
1798        if((jj%2==0 || jj==0) && eDEDXpmt[jj]!=-1 && eDEDXpmt[jj+1]!=-1){
1799          dedxpdl=0.5*(eDEDXpmt[jj]+eDEDXpmt[jj+1]);
1800        }else if((jj%2==0 || jj==0) && eDEDXpmt[jj]!=-1 && eDEDXpmt[jj+1]==-1){
1801          dedxpdl=eDEDXpmt[jj];
1802        }else if((jj%2==0 || jj==0) && eDEDXpmt[jj]==-1 && eDEDXpmt[jj+1]!=-1){
1803          dedxpdl=eDEDXpmt[jj+1];
1804        }
1805        
1806        if(jj%2==0 || jj==0)eDEDXpad[k]= (Float_t)dedxpdl;
1807        
1808        if((jj%2==0 || jj==0)&&eDEDXpad[k]!=-1){
1809          if(k>=0&&k<8)eDEDXlayer[0]=eDEDXpad[k];
1810          if(k>=8&&k<14)eDEDXlayer[1]=eDEDXpad[k];
1811          if(k>=14&&k<16)eDEDXlayer[2]=eDEDXpad[k];
1812          if(k>=16&&k<18)eDEDXlayer[3]=eDEDXpad[k];
1813          if(k>=18&&k<21)eDEDXlayer[4]=eDEDXpad[k];
1814          if(k>=21)eDEDXlayer[5]=eDEDXpad[k];
1815        }
1816    
1817        if(eDEDXlayer[0]!=-1&&eDEDXlayer[1]!=-1&&fabs(eDEDXlayer[0]-eDEDXlayer[1])<10)eDEDXplane[0]=0.5*(eDEDXlayer[0]+eDEDXlayer[1]);
1818        else if(eDEDXlayer[0]!=-1&&eDEDXlayer[1]==-1)eDEDXplane[0]=eDEDXlayer[0];
1819        else if(eDEDXlayer[1]!=-1&&eDEDXlayer[0]==-1)eDEDXplane[0]=eDEDXlayer[1];
1820    
1821        if(eDEDXlayer[2]!=-1&&eDEDXlayer[3]!=-1&&fabs(eDEDXlayer[2]-eDEDXlayer[3])<10)eDEDXplane[1]=0.5*(eDEDXlayer[2]+eDEDXlayer[3]);
1822        else if(eDEDXlayer[2]!=-1&&eDEDXlayer[3]==-1)eDEDXplane[1]=eDEDXlayer[2];
1823        else if(eDEDXlayer[3]!=-1&&eDEDXlayer[2]==-1)eDEDXplane[1]=eDEDXlayer[3];
1824    
1825        if(eDEDXlayer[4]!=-1&&eDEDXlayer[5]!=-1&&fabs(eDEDXlayer[4]-eDEDXlayer[5])<10)eDEDXplane[2]=0.5*(eDEDXlayer[4]+eDEDXlayer[5]);
1826        else if(eDEDXlayer[4]!=-1&&eDEDXlayer[5]==-1)eDEDXplane[2]=eDEDXlayer[4];
1827        else if(eDEDXlayer[5]!=-1&&eDEDXlayer[4]==-1)eDEDXplane[2]=eDEDXlayer[5];
1828    
1829      }
1830      
1831    
1832    
1833    };
1834    
1835    
1836    //------------------------------------------------------------------------
1837    void ToFdEdx::PrintTD()
1838    {
1839      for(int i=0; i<48; i++) {  
1840        TArrayF &binx = TDx[i];
1841        TArrayF &biny = TDy[i];
1842        for(int k=0; k<200; k++) {  // bin temporali
1843          printf("%d %d %f %f", i,k, binx[k], biny[k]);
1844          
1845        }
1846      }
1847    }
1848    
1849    
1850    //------------------------------------------------------------------------
1851    void ToFdEdx::Define_PMTsat()
1852    {
1853      Float_t  sat[48] = {
1854        3176.35,3178.19,3167.38,3099.73,3117.00,3126.29,3111.44,3092.27,
1855        3146.48,3094.41,3132.13,3115.37,3099.32,3110.97,3111.80,3143.14,
1856        3106.72,3153.44,3136.00,3188.96,3104.73,3140.45,3073.18,3106.62,
1857        3112.48,3146.92,3127.24,3136.52,3109.59,3112.89,3045.15,3147.26,
1858        3095.92,3121.05,3083.25,3123.62,3150.92,3125.30,3067.60,3160.18,
1859        3119.36,3108.92,3164.77,3133.64,3111.47,3131.98,3128.87,3135.56 };
1860      PMTsat.Set(48,sat);
1861    }
1862    
1863    //------------------------------------------------------------------------
1864    // void ToFdEdx::ReadParTD( Int_t ipmt, const char *fname )
1865    // {
1866    //   printf("read %s\n",fname);
1867    //   if(ipmt<0)  return;
1868    //   if(ipmt>47) return;
1869    //   FILE *fattin = fopen( fname , "r" );
1870    //   Float_t yTD[200],xTD[200];
1871    //   for(int j=0;j<200;j++){
1872    //     float x,y,ym,e;
1873    //     if(fscanf(fattin,"%f %f %f %f",
1874    //            &x, &y, &ym, &e )!=4) break;
1875    //     xTD[j]=x;
1876    //     if(ym>0&&fabs(y-ym)>1)  yTD[j]=ym;
1877    //     else                    yTD[j]=y;
1878    //   }
1879    //   TDx[ipmt].Set(200,xTD);
1880    //   TDy[ipmt].Set(200,yTD);
1881    //   fclose(fattin);
1882    // }
1883    
1884    //------------------------------------------------------------------------
1885    void ToFdEdx::ReadParBBpos( const char *fname )
1886    {
1887      printf("read %s\n",fname);
1888      parBBpos.Set(48);
1889      FILE *fattin = fopen( fname , "r" );
1890      for (int i=0; i<48; i++) {
1891        int   tid=0;
1892        float  tp;
1893        if(fscanf(fattin,"%d %f",
1894                  &tid, &tp )!=2) break;
1895        parBBpos[i]=tp;
1896      }
1897      fclose(fattin);
1898    }
1899    
1900    //------------------------------------------------------------------------
1901    void ToFdEdx::ReadParDesatBB( const char *fname )
1902    {
1903      printf("read %s\n",fname);
1904      FILE *fattin = fopen( fname , "r" );
1905      for (int i=0; i<48; i++) {
1906        int   tid=0;
1907        float  tp[3];
1908        if(fscanf(fattin,"%d %f %f %f",
1909                  &tid, &tp[0], &tp[1], &tp[2] )!=4) break;
1910        parDesatBB[i].Set(3,tp);
1911      }
1912      fclose(fattin);
1913    }
1914    
1915    
1916    //------------------------------------------------------------------------
1917    void ToFdEdx::ReadParBBneg( const char *fname )
1918    
1919    {
1920      printf("read %s\n",fname);
1921      FILE *fattin = fopen( fname , "r" );
1922      for (int i=0; i<48; i++) {
1923        int   tid=0;
1924        float  tp[3];
1925        if(fscanf(fattin,"%d %f %f %f",
1926                  &tid, &tp[0], &tp[1], &tp[2] )!=4) break;
1927        parBBneg[i].Set(3,tp);
1928      }
1929      fclose(fattin);
1930    }
1931    
1932    //------------------------------------------------------------------------
1933    void ToFdEdx::ReadParPos( const char *fname )
1934    {
1935      printf("read %s\n",fname);
1936      FILE *fattin = fopen( fname , "r" );
1937      for (int i=0; i<48; i++) {
1938        int   tid=0;
1939        float  tp[4];
1940        if(fscanf(fattin,"%d %f %f %f %f",
1941                  &tid, &tp[0], &tp[1], &tp[2], &tp[3])!=5) break;
1942        parPos[i].Set(4,tp);
1943      }
1944      fclose(fattin);
1945    }
1946    
1947    //------------------------------------------------------------------------
1948    void ToFdEdx::ReadParAtt( const char *fname )
1949    {
1950      printf("read %s\n",fname);
1951      FILE *fattin = fopen( fname , "r" );
1952      for (int i=0; i<48; i++) {
1953        int   tid=0;
1954        float  tp[6];
1955        if(fscanf(fattin,"%d %f %f %f %f %f %f",
1956                  &tid, &tp[0], &tp[1], &tp[2], &tp[3], &tp[4], &tp[5] )!=7) break;
1957        parAtt[i].Set(6,tp);
1958      }
1959      fclose(fattin);
1960    }
1961    
1962    
1963    
1964    
1965    
1966    
1967    double ToFdEdx::f_att( TArrayF &p, float x )
1968    {
1969      return
1970        p[0] +
1971        p[1]*x +
1972        p[2]*x*x +
1973        p[3]*x*x*x +
1974        p[4]*x*x*x*x +
1975        p[5]*x*x*x*x*x;
1976    }
1977    //------------------------------------------------------------------------
1978    double ToFdEdx::f_att5B( float x )
1979    {
1980      return
1981        101.9409 +
1982        6.643781*x +
1983        0.2765518*x*x +
1984        0.004617647*x*x*x +
1985        0.0006195132*x*x*x*x +
1986        0.00002813734*x*x*x*x*x;
1987    }
1988    
1989    
1990    double ToFdEdx::f_pos( TArrayF &p, float x )
1991    {
1992      return
1993        p[0] +
1994        p[1]*x +
1995        p[2]*x*x +
1996        p[3]*x*x*x;
1997    }
1998    
1999    double ToFdEdx::f_pos5B( float x )
2000    {
2001      return
2002        15.45132 +
2003        0.8369721*x +
2004        0.0005*x*x;
2005    }
2006    
2007    
2008    
2009    double ToFdEdx::f_adcPC( float x )
2010    {
2011      return 28.12+0.6312*x-5.647e-05*x*x+3.064e-08*x*x*x;
2012    }
2013    
2014    
2015    float ToFdEdx::Get_adc_he( int id, float pl_x[6], float pl_y[6])
2016    {
2017    
2018      //
2019      // input: id - pmt [0:47}
2020      //             pl_x - coord x of the tof plane
2021      //             pl_y - coord y
2022    
2023       adc_he = 0;
2024      if( eGeom.GetXY(id)==1 )  adc_he = f_att( (parAtt[id]), pl_x[eGeom.GetPlane(id)] );
2025      if( eGeom.GetXY(id)==2 )  adc_he = f_att( (parAtt[id]), pl_y[eGeom.GetPlane(id)] );
2026      return adc_he;
2027    }
2028    
2029    //------------------------------------------------------------------------
2030    double ToFdEdx::f_BB( TArrayF &p, float x )
2031    {
2032      return  p[0]/(x*x)*(log(x*x/(1-x*x)) - p[1]*x*x - p[2]);
2033    }
2034    
2035    //------------------------------------------------------------------------
2036    double ToFdEdx::f_BB5B( float x )
2037    {
2038      return  0.165797/(x*x)*(log(x*x/(1-x*x)) + 140.481*x*x + 52.9258);
2039    }
2040    //------------------------------------------------------------------------
2041    double ToFdEdx::f_desatBB( TArrayF &p, float x )
2042    {
2043      return
2044        p[0] +
2045        p[1]*x +
2046        p[2]*x*x;
2047    }
2048    
2049    //------------------------------------------------------------------------
2050    double ToFdEdx::f_desatBB5B( float x )
2051    {
2052      return
2053        -2.4 +
2054        0.75*x +
2055        0.009*x*x;
2056    }
2057    
2058    
2059    
2060    
2061    
2062    
2063    
2064    
2065    
2066    
2067    
2068    
2069    
2070    
2071    
2072    
Legend:
Removed from v.1.23  
changed lines
  Added in v.1.26
  ViewVC Help
Powered by ViewVC 1.1.23