--- DarthVader/CalorimeterLevel2/src/CaloLevel0.cpp 2008/02/07 20:02:09 1.16 +++ DarthVader/CalorimeterLevel2/src/CaloLevel0.cpp 2014/01/16 15:29:08 1.34 @@ -56,6 +56,31 @@ clevel1 = &clevel1_; clevel2 = &clevel2_; // +// extern struct FlEventi eventi_; +// extern struct FlGruppo gruppo_; +// extern struct FlGruppo2 gruppo2_; +// extern struct FlGruppo4 gruppo4_; +// extern struct FlTaglioen taglioen_; +// extern struct FlAngolo angolo_; +// extern struct FlWhere where_; +// extern struct FlGeneral general_; +// extern struct FlCh ch_; +// extern struct FlCalofit calofit_; +// extern struct FlPawcd pawcd_; +// extern struct FlQuestd questd_; +// eventi = &eventi_; +// gruppo = &gruppo_; +// gruppo2 = &gruppo2_; +// gruppo4 = &gruppo4_; +// taglioen = &taglioen_; +// angolo = &angolo_; +// where = &where_; +// general = &general_; +// ch = &ch_; +// calofit = &calofit_; +// pawcd = &pawcd_; +// questd = &questd_; + // trkseqno = 0; ClearStructs(); // @@ -67,17 +92,25 @@ memset(obadmask, 0, 2*22*96*sizeof(Int_t)); memset(obadpulsemask, 0, 2*22*6*sizeof(Int_t)); memset(ctprecor, 0, 2*22*6*sizeof(Float_t)); + memset(ctsicor, 0, 2*22*9*sizeof(Float_t)); memset(ctneigcor, 0, 2*22*6*sizeof(Float_t)); calopar1 = true; calopar2 = true; calopar3 = true; + calopar4 = true; + calopar5 = true; crosst = true; + mask18 = false; ftcalopar1 = 0; ttcalopar1 = 0; ftcalopar2 = 0; ttcalopar2 = 0; ftcalopar3 = 0; ttcalopar3 = 0; + ftcalopar4 = 0; + ttcalopar4 = 0; + ftcalopar5 = 0; + ttcalopar5 = 0; } void CaloLevel0::SetCrossTalk(Bool_t ct){ @@ -88,6 +121,18 @@ ctground = ct; } +void CaloLevel0::SetCrossTalkType(Int_t ct){ + if ( ct == 0 ) ctground = true; + if ( ct == 1 ){ + ctground = false; + noselfct = false; + }; + if ( ct == 2 ){ + ctground = false; + noselfct = true; + }; +} + void CaloLevel0::SetVerbose(Bool_t ct){ verbose = ct; } @@ -95,6 +140,14 @@ /** * Initialize CaloLevel0 object **/ +void CaloLevel0::ProcessingInit(TSQLServer *dbc, UInt_t hs, Int_t &sgnl, TTree *l0tree, Bool_t isdeb, Bool_t isverb){ + if ( !dbc->IsConnected() ) throw -116; + this->InitDo(dbc,hs,sgnl,l0tree,isdeb,isverb); +} + +/** + * Initialize CaloLevel0 object +**/ void CaloLevel0::ProcessingInit(GL_TABLES *glt, UInt_t hs, Int_t &sgnl, TTree *l0tree, Bool_t isdeb, Bool_t isverb){ // const TString host = glt->CGetHost(); @@ -102,6 +155,14 @@ const TString psw = glt->CGetPsw(); TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data()); if ( !dbc->IsConnected() ) throw -116; + this->InitDo(dbc,hs,sgnl,l0tree,isdeb,isverb); + dbc->Close(); + delete dbc; + dbc = 0; +} + + +void CaloLevel0::InitDo(TSQLServer *dbc, UInt_t hs, Int_t &sgnl, TTree *l0tree, Bool_t isdeb, Bool_t isverb){ stringstream myquery; myquery.str(""); myquery << "SET time_zone='+0:00'"; @@ -182,8 +243,6 @@ // delete glcalo; delete glroot; - dbc->Close(); - delete dbc; // return; // @@ -200,6 +259,11 @@ return(sgnl); } +Int_t CaloLevel0::ChkParam(TSQLServer *dbc, UInt_t runheader, Bool_t mechal){ + Int_t sig = this->ChkParamDo(dbc,runheader,mechal); + return(sig); +} + Int_t CaloLevel0::ChkParam(GL_TABLES *glt, UInt_t runheader, Bool_t mechal){ const TString host = glt->CGetHost(); const TString user = glt->CGetUser(); @@ -211,6 +275,15 @@ myquery << "SET time_zone='+0:00'"; dbc->Query(myquery.str().c_str()); // + Int_t sig = this->ChkParamDo(dbc,runheader,mechal); + dbc->Close(); + delete dbc; + dbc = 0; + return(sig); +} + +Int_t CaloLevel0::ChkParamDo(TSQLServer *dbc, UInt_t runheader, Bool_t mechal){ + // stringstream calfile; stringstream bmfile; stringstream aligfile; @@ -221,18 +294,48 @@ // if ( calopar1 || ( ttcalopar1 != 0 && ttcalopar1 < runheader ) ){ // - // - // if ( debug ) printf(" calopar1 %i ftcalopar1 %u ttcalopar1 %u runheader %u \n",calopar1,ftcalopar1,ttcalopar1,runheader); // + if ( calopar1 ){ + // + // determine where I can find calorimeter ADC to MIP conversion file + // + if ( verbose ) printf(" Querying DB for calorimeter parameters files...\n"); + // + error = 0; + error = glparam->Query_GL_PARAM(runheader,101,dbc); + if ( error < 0 ) return(error); + // + calfile.str(""); + calfile << glparam->PATH.Data() << "/"; + calfile << glparam->NAME.Data(); + // + if ( verbose ) printf("\n Using ADC to MIP conversion file: \n %s \n",calfile.str().c_str()); + f = fopen(calfile.str().c_str(),"rb"); + if ( !f ){ + if ( verbose ) printf(" CALORIMETER - ERROR: no ADC to MIP file!\n"); + return(-105); + }; + // + for (Int_t m = 0; m < 2 ; m++ ){ + for (Int_t k = 0; k < 22; k++ ){ + for (Int_t l = 0; l < 96; l++ ){ + fread(&mip[m][k][l],sizeof(mip[m][k][l]),1,f); + if ( debug ) printf(" %f \n",mip[m][k][l]); + }; + }; + }; + fclose(f); + }; + // calopar1 = false; // - // determine where I can find calorimeter ADC to MIP conversion file + // flight extra corrections: // - if ( verbose ) printf(" Querying DB for calorimeter parameters files...\n"); + if ( verbose ) printf(" Querying DB for calorimeter flight ADC to MIP files...\n"); // error = 0; - error = glparam->Query_GL_PARAM(runheader,101,dbc); + error = glparam->Query_GL_PARAM(runheader,110,dbc); if ( error < 0 ) return(error); // calfile.str(""); @@ -241,25 +344,28 @@ ftcalopar1 = glparam->FROM_TIME; ttcalopar1 = glparam->TO_TIME; // - if ( verbose ) printf("\n Using ADC to MIP conversion file: \n %s \n",calfile.str().c_str()); - f = fopen(calfile.str().c_str(),"rb"); - if ( !f ){ - if ( verbose ) printf(" CALORIMETER - ERROR: no ADC to MIP file!\n"); - return(-105); + if ( verbose ) printf("\n Using ADC to MIP special conversion file: \n %s \n",calfile.str().c_str()); + ifstream spfile; + spfile.open(calfile.str().c_str()); + if ( !spfile ){ + if ( verbose ) printf(" CALORIMETER - ERROR: no special calibration file!\n"); + return(-123); }; + // + Int_t vview = 0; + Int_t vplane = 0; + Int_t vstrip = 0; + Float_t vval = 0.; + while ( spfile >> vview && spfile >> vplane && spfile >> vstrip && spfile >> vval){ + if ( debug ) printf(" Setting ADC to MIP conversion factor: view %i plane %i strip %i mip %f \n",vview,vplane,vstrip,vval); + mip[vview][vplane][vstrip] = vval; + }; // - for (Int_t m = 0; m < 2 ; m++ ){ - for (Int_t k = 0; k < 22; k++ ){ - for (Int_t l = 0; l < 96; l++ ){ - fread(&mip[m][k][l],sizeof(mip[m][k][l]),1,f); - if ( debug ) printf(" %f \n",mip[m][k][l]); - }; - }; - }; - fclose(f); }; // + // if ( calopar2 || ( ttcalopar2 != 0 && ttcalopar2 < runheader ) ){ + // if ( debug ) printf(" calopar2 %i ftcalopar2 %u ttcalopar2 %u runheader %u \n",calopar2,ftcalopar2,ttcalopar2,runheader); calopar2 = false; // @@ -366,17 +472,145 @@ badfile.close(); }; // + // calopar4 + // + if ( calopar4 || ( ttcalopar4 != 0 && ttcalopar4 < runheader ) ){ + // + if ( debug ) printf(" calopar4 %i ftcalopar4 %u ttcalopar4 %u runheader %u \n",calopar4,ftcalopar4,ttcalopar4,runheader); + // + calopar4 = false; + // + // flight extra corrections: + // + if ( verbose ) printf(" Querying DB for calorimeter max rms file...\n"); + // + error = 0; + error = glparam->Query_GL_PARAM(runheader,109,dbc); + if ( error < 0 ) return(error); + // + calfile.str(""); + calfile << glparam->PATH.Data() << "/"; + calfile << glparam->NAME.Data(); + ftcalopar4 = glparam->FROM_TIME; + ttcalopar4 = glparam->TO_TIME; + // + if ( verbose ) printf("\n Using calorimeter max rms file: \n %s \n",calfile.str().c_str()); + ifstream spfile; + spfile.open(calfile.str().c_str()); + if ( !spfile ){ + if ( verbose ) printf(" CALORIMETER - ERROR: no max rms file!\n"); + return(-124); + }; + // + Int_t vview = 0; + Int_t vplane = 0; + Int_t vval = 0; + for (Int_t l=0; l<2; l++){ + for (Int_t m=0; m<22; m++){ + maxrms[l][m] = 26; + }; + }; + while ( spfile >> vview && spfile >> vplane && spfile >> vval){ + if ( debug ) printf(" Setting view %i plane %i max rms %i \n",vview,vplane,vval); + maxrms[vview][vplane] = vval; + }; + spfile.close(); + // + }; + // + // calopar5 + // + if ( calopar5 || ( ttcalopar5 != 0 && ttcalopar5 < runheader ) ){ + // + if ( debug ) printf(" calopar5 %i ftcalopar5 %u ttcalopar5 %u runheader %u \n",calopar5,ftcalopar5,ttcalopar5,runheader); + // + calopar5 = false; + // + // flight extra corrections: + // + if ( verbose ) printf(" Querying DB for calorimeter noise to signal threshold file...\n"); + // + error = 0; + error = glparam->Query_GL_PARAM(runheader,111,dbc); + if ( error < 0 ) return(error); + // + calfile.str(""); + calfile << glparam->PATH.Data() << "/"; + calfile << glparam->NAME.Data(); + ftcalopar5 = glparam->FROM_TIME; + ttcalopar5 = glparam->TO_TIME; + // + if ( verbose ) printf("\n Using calorimeter noise to signal threshold file: \n %s \n",calfile.str().c_str()); + ifstream spfile; + spfile.open(calfile.str().c_str()); + if ( !spfile ){ + if ( verbose ) printf(" CALORIMETER - ERROR: no noise to signal threshold file!\n"); + return(-125); + }; + // + Int_t vview = 0; + Int_t vplane = 0; + Int_t vstrip = 0; + Float_t vval = 0.; + for (Int_t l=0; l<2; l++){ + for (Int_t m=0; m<22; m++){ + for (Int_t n=0; n<96; n++){ + memin[l][m][n] = 0.7; + }; + }; + }; + while ( spfile >> vview && spfile >> vplane && spfile >> vstrip && spfile >> vval){ + if ( vstrip == -1 ){ + for (Int_t ll=0; ll<96; ll++){ + if ( debug ) printf(" Setting view %i plane %i strip %i noise to signal ratio %f \n",vview,vplane,ll,vval); + memin[vview][vplane][ll] = vval; + }; + } else { + if ( debug ) printf(" Setting view %i plane %i strip %i noise to signal ratio %f \n",vview,vplane,vstrip,vval); + memin[vview][vplane][vstrip] = vval; + }; + }; + spfile.close(); + // + }; + // + // delete glparam; - dbc->Close(); - delete dbc; // return(0); } -Int_t CaloLevel0::CalcCrossTalkCorr(GL_TABLES *glt, UInt_t runheader){ +Int_t CaloLevel0::CalcCrossTalkCorr(TSQLServer *dbc, UInt_t runheader, Bool_t ctusetable){ + Int_t sig = CalcCrossTalkCorrDo(dbc,runheader,ctusetable); + return(sig); +}; + +Int_t CaloLevel0::CalcCrossTalkCorr(TSQLServer *dbc, UInt_t runheader){ + Int_t sig = CalcCrossTalkCorrDo(dbc,runheader,true); + return(sig); +} + +Int_t CaloLevel0::CalcCrossTalkCorr(GL_TABLES *glt, UInt_t runheader, Bool_t usetable){ + const TString host = glt->CGetHost(); + const TString user = glt->CGetUser(); + const TString psw = glt->CGetPsw(); + TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data()); + if ( !dbc->IsConnected() ) throw -116; + stringstream myquery; + myquery.str(""); + myquery << "SET time_zone='+0:00'"; + dbc->Query(myquery.str().c_str()); // - if ( ctground ) return(0); + Int_t sig = CalcCrossTalkCorrDo(dbc,runheader,usetable); + dbc->Close(); + delete dbc; + dbc = 0; // + return(sig); + // +}; + +Int_t CaloLevel0::CalcCrossTalkCorr(GL_TABLES *glt, UInt_t runheader){ const TString host = glt->CGetHost(); const TString user = glt->CGetUser(); const TString psw = glt->CGetPsw(); @@ -387,13 +621,25 @@ myquery << "SET time_zone='+0:00'"; dbc->Query(myquery.str().c_str()); // - stringstream bmfile; + Int_t sig = CalcCrossTalkCorrDo(dbc,runheader,true); + dbc->Close(); + delete dbc; + dbc = 0; + // + return(sig); + // +} + +Int_t CaloLevel0::CalcCrossTalkCorrDo(TSQLServer *dbc, UInt_t runheader, Bool_t usetable){ + // + if ( ctground ) return(0); + // Int_t error = 0; - ifstream badfile; - GL_PARAM *glparam = new GL_PARAM(); + GL_PARAM *glparam = new GL_PARAM(); // // determine where I can find file with offline bad pulser mask // + stringstream bmfile; error = 0; error = glparam->Query_GL_PARAM(runheader,105,dbc); if ( error < 0 ) return(error); @@ -402,6 +648,7 @@ bmfile << glparam->PATH.Data() << "/"; bmfile << glparam->NAME.Data(); // + ifstream badfile; if ( verbose ) printf("\n Using bad pulser offline mask file: \n %s \n\n",bmfile.str().c_str()); badfile.open(bmfile.str().c_str()); if ( !badfile ){ @@ -430,259 +677,350 @@ }; }; // - delete glparam; badfile.close(); - // - // Let's start with cross-talk correction calculation - // - GL_CALOPULSE_CALIB *glp = new GL_CALOPULSE_CALIB(); - Float_t adcp[2][22][96]; - Float_t adcpcal[2][22][96]; - memset(adcp , 0, 2*22*96*sizeof(Float_t)); - memset(adcpcal , 0, 2*22*96*sizeof(Float_t)); - // - UInt_t pampli = 0; - for (Int_t s=0; s<4; s++){ + if ( !usetable ){ // - // Save into matrix adcp the values of the highest pulse calibration (pulse amplitude = 2) + // Let's start with cross-talk correction calculation // - pampli = 2; - error = 0; - error = glp->Query_GL_CALOPULSE_CALIB(runheader,s,pampli,dbc); - if ( error < 0 ){ - if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n"); - return(error); - }; - // - UInt_t idcalib = glp->ID_ROOT_L0; - UInt_t fromtime = glp->FROM_TIME; - UInt_t calibno = glp->EV_ROOT; - // - // determine path and name and entry of the calibration file + GL_CALOPULSE_CALIB *glp = new GL_CALOPULSE_CALIB(); + Float_t adcp[2][22][96]; + Float_t adcpcal[2][22][96]; + memset(adcp , 0, 2*22*96*sizeof(Float_t)); + memset(adcpcal , 0, 2*22*96*sizeof(Float_t)); // - GL_ROOT *glroot = new GL_ROOT(); - if ( verbose ) printf("\n"); - if ( verbose ) printf(" ** SECTION %i **\n",s); - // - error = 0; - error = glroot->Query_GL_ROOT(idcalib,dbc); - if ( error < 0 ){ - if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n"); - return(error); - }; - // - stringstream name; - name.str(""); - name << glroot->PATH.Data() << "/"; - name << glroot->NAME.Data(); - // - TString fcalname = (TString)name.str().c_str(); - ifstream myfile; - myfile.open(fcalname.Data()); - if ( !myfile ){ - return(-107); - }; - myfile.close(); - // - TFile *File = new TFile(fcalname.Data()); - if ( !File ) return(-108); - TTree *tr = (TTree*)File->Get("CalibCalPulse2"); - if ( !tr ) return(-119); - // - TBranch *calo = tr->GetBranch("CalibCalPulse2"); - // - pamela::CalibCalPulse2Event *ce = 0; - tr->SetBranchAddress("CalibCalPulse2", &ce); - // - Long64_t ncalibs = calo->GetEntries(); + UInt_t pampli = 0; + for (Int_t s=0; s<4; s++){ + // + // Save into matrix adcp the values of the highest pulse calibration (pulse amplitude = 2) + // + pampli = 2; + error = 0; + error = glp->Query_GL_CALOPULSE_CALIB(runheader,s,pampli,dbc); + if ( error < 0 ){ + if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n"); + return(error); + }; + // + UInt_t idcalib = glp->ID_ROOT_L0; + UInt_t fromtime = glp->FROM_TIME; + UInt_t calibno = glp->EV_ROOT; + // + // determine path and name and entry of the calibration file + // + GL_ROOT *glroot = new GL_ROOT(); + if ( verbose ) printf("\n"); + if ( verbose ) printf(" ** SECTION %i **\n",s); + // + error = 0; + error = glroot->Query_GL_ROOT(idcalib,dbc); + if ( error < 0 ){ + if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n"); + return(error); + }; + // + stringstream name; + name.str(""); + name << glroot->PATH.Data() << "/"; + name << glroot->NAME.Data(); + // + TString fcalname = (TString)name.str().c_str(); + ifstream myfile; + myfile.open(fcalname.Data()); + if ( !myfile ){ + return(-107); + }; + myfile.close(); + // + TFile *File = new TFile(fcalname.Data()); + if ( !File ) return(-108); + TTree *tr = (TTree*)File->Get("CalibCalPulse2"); + if ( !tr ) return(-119); + // + TBranch *calo = tr->GetBranch("CalibCalPulse2"); + // + pamela::CalibCalPulse2Event *ce = 0; + tr->SetBranchAddress("CalibCalPulse2", &ce); + // + Long64_t ncalibs = calo->GetEntries(); + // + if ( !ncalibs ) return(-110); + // + if ( calo->GetEntry(calibno) <= 0) throw -36; + if ( verbose ) printf(" PULSE2 using entry %u from file %s",calibno,fcalname.Data()); + // + // retrieve calibration table + // + if ( ce->pstwerr[s] && ce->pperror[s] == 0 && ce->unpackError == 0 ){ + for ( Int_t d=0 ; d<11 ;d++ ){ + for ( Int_t j=0; j<96 ;j++){ + if ( s == 2 ){ + adcp[0][2*d+1][j] = ce->calpuls[3][d][j]; + }; + if ( s == 3 ){ + adcp[0][2*d][j] = ce->calpuls[1][d][j]; + }; + if ( s == 0 ){ + adcp[1][2*d][j] = ce->calpuls[0][d][j]; + }; + if ( s == 1 ){ + adcp[1][2*d+1][j] = ce->calpuls[2][d][j]; + }; + }; + }; + } else { + if ( verbose ) printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n "); + return(-111); + }; + // + File->Close(); + delete glroot; + // + // Save into matrix adcpcal the calibrated values of the pulse calibration (subtraction of pulse amplitude = 0 relative to the pulse2 calibration used) + // + pampli = 0; + error = 0; + error = glp->Query_GL_CALOPULSE_CALIB(fromtime,s,pampli,dbc); + if ( error < 0 ){ + if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n"); + return(error); + }; + // + idcalib = glp->ID_ROOT_L0; + calibno = glp->EV_ROOT; + // + // determine path and name and entry of the calibration file + // + glroot = new GL_ROOT(); + if ( verbose ) printf("\n"); + if ( verbose ) printf(" ** SECTION %i **\n",s); + // + error = 0; + error = glroot->Query_GL_ROOT(idcalib,dbc); + if ( error < 0 ){ + if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n"); + return(error); + }; + // + name.str(""); + name << glroot->PATH.Data() << "/"; + name << glroot->NAME.Data(); + // + fcalname = (TString)name.str().c_str(); + myfile.open(fcalname.Data()); + if ( !myfile ){ + return(-107); + }; + myfile.close(); + // + TFile *File1 = new TFile(fcalname.Data()); + if ( !File1 ) return(-108); + TTree *tr1 = (TTree*)File1->Get("CalibCalPulse1"); + if ( !tr1 ) return(-120); + // + TBranch *calo1 = tr1->GetBranch("CalibCalPulse1"); + // + pamela::CalibCalPulse1Event *ce1 = 0; + tr1->SetBranchAddress("CalibCalPulse1", &ce1); + // + ncalibs = calo1->GetEntries(); + // + if ( !ncalibs ) return(-110); + // + if ( calo1->GetEntry(calibno) <= 0 ) throw -36; + if ( verbose ) printf(" PULSE1 using entry %u from file %s",calibno,fcalname.Data()); + // + // retrieve calibration table + // + if ( ce1->pstwerr[s] && ce1->pperror[s] == 0 && ce1->unpackError == 0 ){ + for ( Int_t d=0 ; d<11 ;d++ ){ + for ( Int_t j=0; j<96 ;j++){ + if ( s == 2 ){ + adcpcal[0][2*d+1][j] = adcp[0][2*d+1][j] - ce1->calpuls[3][d][j]; + }; + if ( s == 3 ){ + adcpcal[0][2*d][j] = adcp[0][2*d][j] - ce1->calpuls[1][d][j]; + }; + if ( s == 0 ){ + adcpcal[1][2*d][j] = adcp[1][2*d][j] - ce1->calpuls[0][d][j]; + }; + if ( s == 1 ){ + adcpcal[1][2*d+1][j] = adcp[1][2*d+1][j] - ce1->calpuls[2][d][j]; + }; + }; + }; + } else { + if ( verbose ) printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n "); + return(-111); + }; + // + File1->Close(); + // + delete glroot; + // + };// loop on the four sections // - if ( !ncalibs ) return(-110); // - calo->GetEntry(calibno); - if ( verbose ) printf(" PULSE2 using entry %u from file %s",calibno,fcalname.Data()); + delete glp; // - // retrieve calibration table + // Ok, now we can try to calculate the cross-talk correction for each pre-amplifier // - if ( ce->pstwerr[s] && ce->pperror[s] == 0 && ce->unpackError == 0 ){ - for ( Int_t d=0 ; d<11 ;d++ ){ - for ( Int_t j=0; j<96 ;j++){ - if ( s == 2 ){ - adcp[0][2*d+1][j] = ce->calpuls[3][d][j]; - }; - if ( s == 3 ){ - adcp[0][2*d][j] = ce->calpuls[1][d][j]; - }; - if ( s == 0 ){ - adcp[1][2*d][j] = ce->calpuls[0][d][j]; + for ( Int_t v=0; v<2; v++){ + if ( debug ) printf(" \n\n NEW VIEW \n"); + for ( Int_t p=0; p<22; p++){ + for ( Int_t npre=0; npre<6; npre++){ + ctprecor[v][p][npre] = 1000.; + ctneigcor[v][p][npre] = 1000.; + Int_t str0=npre*16; + Int_t str16= -1 + (1+npre)*16; + // + UInt_t neigc = 0; + UInt_t farc = 0; + UInt_t pulsc = 0; + Float_t sigpulsed = 0.; + Float_t neigbase = 0.; + Float_t farbase = 0.; + // + // Loop over the strip of the pre and sum all signal far away from pulsed strip, signal in the neighbour(s) strip(s) and save the pulsed signal + // moreover count the number of strips in each case + // + for (Int_t s=str0; s<=str16; s++){ + if ( adcpcal[v][p][s] > 10000.){ + sigpulsed = adcpcal[v][p][s]; + pulsc++; + if ( s > str0 ){ + neigbase += adcpcal[v][p][s-1]; + neigc++; + farbase -= adcpcal[v][p][s-1]; + farc--; + }; + if ( s < str16 ){ + neigbase += adcpcal[v][p][s+1]; + neigc++; + farbase -= adcpcal[v][p][s+1]; + farc--; + }; + } else { + farc++; + farbase += adcpcal[v][p][s]; + }; }; - if ( s == 1 ){ - adcp[1][2*d+1][j] = ce->calpuls[2][d][j]; + // + // Now calculate the corrections + // + Float_t avefarbase = 0.; + if ( farc ) avefarbase = farbase/(Float_t)farc; + Float_t aveneigbase = 0.; + if ( neigc ) aveneigbase = neigbase/(Float_t)neigc; + // + if ( pulsc == 1 && farc && neigc ){ + ctprecor[v][p][npre] = -avefarbase/(sigpulsed+fabs(avefarbase)); + ctneigcor[v][p][npre] = fabs(aveneigbase-avefarbase)/(sigpulsed+fabs(avefarbase)); + if ( debug ) printf(" Cross-talk correction View %i Plane %i Pre %i : pre-correction: %f neighbour strips correction %f \n",v,p,npre,ctprecor[v][p][npre],ctneigcor[v][p][npre]); + } else { + // + // did not find the pulsed strip or more than one pulsed strip found! + // + if ( debug ) printf(" Problems finding the cross-talk corrections: \n View %i Plane %i Pre %i number of pulsed strip %i \n Average faraway baseline %f number of strips %i Average neighbour baseline %f number of neighbour strips %i \n",v,p,npre,pulsc,avefarbase,farc,aveneigbase,neigc); + // }; }; + if ( debug ) printf(" \n ==================== \n"); }; - } else { - if ( verbose ) printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n "); - return(-111); }; + } else { // - File->Close(); - delete glroot; + // use pre-amply table + // // - // Save into matrix adcpcal the calibrated values of the pulse calibration (subtraction of pulse amplitude = 0 relative to the pulse2 calibration used) + // determine where I can find file with offline neighbour correction table // - pampli = 0; + stringstream bmfile2; error = 0; - error = glp->Query_GL_CALOPULSE_CALIB(fromtime,s,pampli,dbc); - if ( error < 0 ){ - if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n"); - return(error); - }; - // - idcalib = glp->ID_ROOT_L0; - calibno = glp->EV_ROOT; + error = glparam->Query_GL_PARAM(runheader,106,dbc); + if ( error < 0 ) return(error); // - // determine path and name and entry of the calibration file + bmfile2.str(""); + bmfile2 << glparam->PATH.Data() << "/"; + bmfile2 << glparam->NAME.Data(); + // + ifstream badfile2; + if ( verbose ) printf("\n Using pre-amply neighbour crosstalk table file: \n %s \n\n",bmfile2.str().c_str()); + badfile2.open(bmfile2.str().c_str()); + if ( !badfile2 ){ + if ( verbose ) printf(" CALORIMETER - ERROR: no pre-amply neighbour crosstalk table file!\n"); + return(-121); + }; + // + Int_t vview = 0; + Int_t vplane = 0; + Int_t vpre = 0; + Float_t vcorr = 0.; + while ( badfile2 >> vview && badfile2 >> vplane && badfile2 >> vpre && badfile2 >> vcorr){ + if ( debug ) printf(" Pre-amply neighbour correction: view %i plane %i pre %i correction %f \n",vview,vplane,vpre,vcorr); + ctneigcor[vview][vplane][vpre] = vcorr; + }; // - glroot = new GL_ROOT(); - if ( verbose ) printf("\n"); - if ( verbose ) printf(" ** SECTION %i **\n",s); + // determine where I can find file with offline SECOND neighbour correction table // + stringstream bmfile3; error = 0; - error = glroot->Query_GL_ROOT(idcalib,dbc); - if ( error < 0 ){ - if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n"); - return(error); - }; - // - name.str(""); - name << glroot->PATH.Data() << "/"; - name << glroot->NAME.Data(); + error = glparam->Query_GL_PARAM(runheader,107,dbc); + if ( error < 0 ) return(error); // - fcalname = (TString)name.str().c_str(); - myfile.open(fcalname.Data()); - if ( !myfile ){ - return(-107); + bmfile3.str(""); + bmfile3 << glparam->PATH.Data() << "/"; + bmfile3 << glparam->NAME.Data(); + // + ifstream badfile3; + if ( verbose ) printf("\n Using pre-amply second neighbour crosstalk table file: \n %s \n\n",bmfile3.str().c_str()); + badfile3.open(bmfile3.str().c_str()); + if ( !badfile3 ){ + if ( verbose ) printf(" CALORIMETER - ERROR: no pre-amply second neighbour crosstalk table file!\n"); + return(-122); }; - myfile.close(); - // - TFile *File1 = new TFile(fcalname.Data()); - if ( !File1 ) return(-108); - TTree *tr1 = (TTree*)File1->Get("CalibCalPulse1"); - if ( !tr1 ) return(-120); - // - TBranch *calo1 = tr1->GetBranch("CalibCalPulse1"); - // - pamela::CalibCalPulse1Event *ce1 = 0; - tr1->SetBranchAddress("CalibCalPulse1", &ce1); - // - ncalibs = calo1->GetEntries(); - // - if ( !ncalibs ) return(-110); + // + Int_t pview = 0; + Int_t pplane = 0; + Int_t ppre = 0; + Float_t pcorr = 0.; + while ( badfile3 >> pview && badfile3 >> pplane && badfile3 >> ppre && badfile3 >> pcorr){ + if ( debug ) printf(" Pre-amply second neighbour correction: view %i plane %i pre %i correction %f \n",pview,pplane,ppre,-pcorr); + ctprecor[pview][pplane][ppre] = -pcorr; // data are saved as negatives in the file + }; // - calo1->GetEntry(calibno); - if ( verbose ) printf(" PULSE1 using entry %u from file %s",calibno,fcalname.Data()); + // determine where to find the file containing the Silicon crosstalk correction table // - // retrieve calibration table + stringstream bmfile4; + error = 0; + error = glparam->Query_GL_PARAM(runheader,108,dbc); + if ( error < 0 ) return(error); // - if ( ce1->pstwerr[s] && ce1->pperror[s] == 0 && ce1->unpackError == 0 ){ - for ( Int_t d=0 ; d<11 ;d++ ){ - for ( Int_t j=0; j<96 ;j++){ - if ( s == 2 ){ - adcpcal[0][2*d+1][j] = adcp[0][2*d+1][j] - ce1->calpuls[3][d][j]; - }; - if ( s == 3 ){ - adcpcal[0][2*d][j] = adcp[0][2*d][j] - ce1->calpuls[1][d][j]; - }; - if ( s == 0 ){ - adcpcal[1][2*d][j] = adcp[1][2*d][j] - ce1->calpuls[0][d][j]; - }; - if ( s == 1 ){ - adcpcal[1][2*d+1][j] = adcp[1][2*d+1][j] - ce1->calpuls[2][d][j]; - }; - }; - }; - } else { - if ( verbose ) printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n "); - return(-111); + bmfile4.str(""); + bmfile4 << glparam->PATH.Data() << "/"; + bmfile4 << glparam->NAME.Data(); + // + ifstream badfile4; + if ( verbose ) printf("\n Using Silicon crosstalk table file: \n %s \n\n",bmfile4.str().c_str()); + badfile4.open(bmfile4.str().c_str()); + if ( !badfile4 ){ + if ( verbose ) printf(" CALORIMETER - ERROR: no Silicon crosstalk table file!\n"); + return(-125); }; + // + Int_t spview = 0; + Int_t spplane = 0; + Int_t psil = 0; + Float_t spcorr = 0.; + memset(ctsicor, 0, 2*22*9*sizeof(Float_t)); + while ( badfile4 >> spview && badfile4 >> spplane && badfile4 >> psil && badfile4 >> spcorr){ + if ( debug ) printf(" Silicon correction: view %i plane %i silicon %i correction %f \n",spview,spplane,psil,-spcorr); + ctsicor[spview][spplane][psil] = -spcorr; // data are saved as negatives in the file + }; // - File1->Close(); - // - delete glroot; - // - };// loop on the four sections - // - // - delete glp; - dbc->Close(); - delete dbc; - // - // Ok, now we can try to calculate the cross-talk correction for each pre-amplifier - // - for ( Int_t v=0; v<2; v++){ - if ( debug ) printf(" \n\n NEW VIEW \n"); - for ( Int_t p=0; p<22; p++){ - for ( Int_t npre=0; npre<6; npre++){ - ctprecor[v][p][npre] = 1000.; - ctneigcor[v][p][npre] = 1000.; - Int_t str0=npre*16; - Int_t str16= -1 + (1+npre)*16; - // - UInt_t neigc = 0; - UInt_t farc = 0; - UInt_t pulsc = 0; - Float_t sigpulsed = 0.; - Float_t neigbase = 0.; - Float_t farbase = 0.; - // - // Loop over the strip of the pre and sum all signal far away from pulsed strip, signal in the neighbour(s) strip(s) and save the pulsed signal - // moreover count the number of strips in each case - // - for (Int_t s=str0; s<=str16; s++){ - if ( adcpcal[v][p][s] > 10000.){ - sigpulsed = adcpcal[v][p][s]; - pulsc++; - if ( s > str0 ){ - neigbase += adcpcal[v][p][s-1]; - neigc++; - farbase -= adcpcal[v][p][s-1]; - farc--; - }; - if ( s < str16 ){ - neigbase += adcpcal[v][p][s+1]; - neigc++; - farbase -= adcpcal[v][p][s+1]; - farc--; - }; - } else { - farc++; - farbase += adcpcal[v][p][s]; - }; - }; - // - // Now calculate the corrections - // - Float_t avefarbase = 0.; - if ( farc ) avefarbase = farbase/(Float_t)farc; - Float_t aveneigbase = 0.; - if ( neigc ) aveneigbase = neigbase/(Float_t)neigc; - // - if ( pulsc == 1 && farc && neigc ){ - ctprecor[v][p][npre] = -avefarbase/(sigpulsed+fabs(avefarbase)); - ctneigcor[v][p][npre] = fabs(aveneigbase-avefarbase)/(sigpulsed+fabs(avefarbase)); - if ( debug ) printf(" Cross-talk correction View %i Plane %i Pre %i : pre-correction: %f neighbour strips correction %f \n",v,p,npre,ctprecor[v][p][npre],ctneigcor[v][p][npre]); - } else { - // - // did not find the pulsed strip or more than one pulsed strip found! - // - if ( debug ) printf(" Problems finding the cross-talk corrections: \n View %i Plane %i Pre %i number of pulsed strip %i \n Average faraway baseline %f number of strips %i Average neighbour baseline %f number of neighbour strips %i \n",v,p,npre,pulsc,avefarbase,farc,aveneigbase,neigc); - // - }; - }; - if ( debug ) printf(" \n ==================== \n"); - }; }; // + delete glparam; + // // Check the calculated corrections // Int_t opre=0; @@ -723,82 +1061,149 @@ }; }; }; - }; + }; // return(0); } void CaloLevel0::FindBaseCompress(Int_t l, Int_t m, Int_t pre){ + Int_t n = 0; + Float_t q = 0; + this->FindBaseCompress(l,m,pre,n,q); +} + +void CaloLevel0::FindBaseCompress(Int_t l, Int_t m, Int_t pre, Int_t &nst, Float_t &qp){ for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ dexy[l][m][e] = dexyc[l][m][e]; }; - this->FindBaseRaw(l,m,pre); + this->FindBaseRaw(l,m,pre,nst,qp); } void CaloLevel0::FindBaseRaw(Int_t l, Int_t m, Int_t pre){ - Float_t minstrip = 100000.; - Float_t rms = 0.; + Int_t n = 0; + Float_t q = 0; + this->FindBaseRaw(l,m,pre,n,q); +} + +void CaloLevel0::FindBaseRaw(Int_t l, Int_t m, Int_t pre, Int_t &nst, Float_t &qp){ + // + Float_t minstrip = 100000.; + Float_t rms = 0.; + Int_t process = 0; + Int_t onlmask[16]; + memset(onlmask, 0, 16*sizeof(Int_t)); + // + while ( process < 2 ){ + // + minstrip = 100000.; + rms = 0.; base[l][m][pre] = 0.; + qp = 0.; + // + Int_t spos = -1; + Int_t ee = 0; for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ - if ( calgood[l][m][e] == 0. && obadmask[l][m][e] == 0 && dexy[l][m][e]-calped[l][m][e] < minstrip && dexy[l][m][e] > 0.) { - minstrip = dexy[l][m][e]-calped[l][m][e]; - rms = calthr[l][m][pre]; - }; + if ( calgood[l][m][e] == 0. && obadmask[l][m][e] == 0 && dexy[l][m][e]-calped[l][m][e] < minstrip && dexy[l][m][e] > 0. && onlmask[ee] == 0 ) { + minstrip = dexy[l][m][e]-calped[l][m][e]; + rms = calthr[l][m][pre]; + spos = ee; + }; + ee++; + qp += (dexy[l][m][e]-calped[l][m][e]-sbase[l][m][e]); }; - if ( debug && l==1 ){ - printf("\n BASELINE CALCULATION for view %i pl %i pre %i: \n => minstrip %f rms %f \n",l,m,pre,minstrip,rms); + // + if ( debug && l==0 ){ + printf("\n BASELINE CALCULATION for view %i pl %i pre %i: \n => minstrip %f rms %f \n => qp = %f \n",l,m,pre,minstrip,rms,qp); }; if ( minstrip != 100000. ) { - Float_t strip6s = 0.; - for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ - if ( (dexy[l][m][e]-calped[l][m][e]) >= minstrip && (dexy[l][m][e]-calped[l][m][e]) <= (minstrip+rms) ) { - strip6s += 1.; - base[l][m][pre] += (dexy[l][m][e] - calped[l][m][e]); - }; - // - // compression - // - if ( abs((int)(dexy[l][m][e]-calped[l][m][e])) <= (minstrip+rms) ) { - dexyc[l][m][e] = 0.; - } else { - dexyc[l][m][e] = dexy[l][m][e]; - }; - }; - if ( debug && l==1 ){ - printf(" strip6s %f \n",strip6s); + Float_t strip6s = 0.; + for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ + if ( (dexy[l][m][e]-calped[l][m][e]) >= minstrip && (dexy[l][m][e]-calped[l][m][e]) <= (minstrip+rms) ) { + strip6s += 1.; + base[l][m][pre] += (dexy[l][m][e] - calped[l][m][e]); }; - // if ( strip6s >= 9. ){ - if ( strip6s >= 2. ){ - Double_t arro = base[l][m][pre]/strip6s; - Float_t deci = 1000.*((float)arro - float(int(arro))); - if ( deci < 500. ) { - arro = double(int(arro)); - } else { - arro = 1. + double(int(arro)); - }; - base[l][m][pre] = arro; + // + // compression + // + // if ( abs((int)(dexy[l][m][e]-calped[l][m][e])) <= (minstrip+rms) ) { + // dexyc[l][m][e] = 0.; + // } else { + dexyc[l][m][e] = dexy[l][m][e]; + // }; + }; + // + if ( strip6s == 1. && process < 1 ){ + onlmask[spos] = 1; + process++; + if ( debug ) printf(" Warning, only one strip to calculate baseline: minstrip %f rms %f spos %i l %i m %i pre %i \n",minstrip,rms,spos,l,m,pre); + continue; + }; + process += 2; + nst = (Int_t)strip6s; + // + if ( debug ){ + printf(" strip6s %f \n",strip6s); + }; + // if ( strip6s >= 9. ){ + if ( (strip6s >= 2. && process == 2) || (strip6s >= 9. && process > 2) ){ + //if ( (strip6s >= 4. && process == 2) || (strip6s >= 9. && process > 2) ){ + Double_t arro = base[l][m][pre]/strip6s; + Float_t deci = 1000.*((float)arro - float(int(arro))); + if ( deci < 500. ) { + arro = double(int(arro)); } else { - base[l][m][pre] = 31000.; - for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ - dexyc[l][m][e] = dexy[l][m][e]; - }; + arro = 1. + double(int(arro)); }; - } else { + base[l][m][pre] = arro; + // + // if too few strips were used to determine the baseline check if it is comparable with the previous event, if not mark it as bad + // + if ( debug && process > 2 ) printf(" AGH low strip value was discarded process %i strip6s %f minstrip %f rms %f spos %i\n",process,strip6s,minstrip,rms,spos); + if ( debug ) printf(" Calculated baseline: base %f sbase-0.02*qp %f \n",base[l][m][pre],(-qp*0.02+sbase[l][m][pre])); + // + if ( strip6s < 4 && base[l][m][pre] > (-0.015*qp+sbase[l][m][pre]) && sbase[l][m][pre] > 0. ){ + if ( debug ) printf(" Suspicious calculated baseline: base %f sbase-0.02*qp %f strip6s %i \n",base[l][m][pre],(-qp*0.02+sbase[l][m][pre]),(Int_t)strip6s); + base[l][m][pre] = 31000.; + nst = 0; // 9RED BUG + qp = 0.; // 9RED BUG + for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ + dexyc[l][m][e] = dexy[l][m][e]; + }; + }; + } else { + if ( debug ) printf(" reset baseline here if ! ( (strip6s >=2 && process == 2) || (strip6s >= 9 and process > 2) ) \n"); base[l][m][pre] = 31000.; + nst = 0; // 9RED BUG + qp = 0.; // 9RED BUG + for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ + dexyc[l][m][e] = dexy[l][m][e]; + }; + }; + } else { + if ( debug ) printf(" reset baseline here if no minimum find\n"); + nst = 0; // 9RED BUG + qp = 0.; // 9RED BUG + process += 2; + base[l][m][pre] = 31000.; + for (Int_t e = pre*16; e < (pre+1)*16 ; e++){ + dexyc[l][m][e] = dexy[l][m][e]; + }; }; + }; + if ( debug ) printf(" Baseline calculation: baseline for view %i plane %i pre %i is %f nst %i qp %f \n",l,m,pre,base[l][m][pre],nst,qp); } Int_t CaloLevel0::Calibrate(Int_t ei){ // // get entry ei // - l0calo->GetEntry(ei); + if ( l0calo->GetEntry(ei) <= 0 ) throw -36; // // if it was not a selftrigger event, could it ever been a selftrigger event? if so trigty = 3. // clevel2->nsatstrip = 0.; Int_t val = 0; - Int_t del = 1100; + Int_t del = 1000; for (Int_t sec = 0; sec < 4; sec++){ for (Int_t dsec = 0; dsec < 7; dsec++){ val = (Int_t)de->calselftrig[sec][dsec]; @@ -807,26 +1212,26 @@ }; }; val = 0; - del = 1100; - if ( clevel2->trigty != 2. ){ - Bool_t ck = false; + del = 1000; + if ( clevel2->trigty < 2. ){ + // Bool_t ck = false; for (Int_t sec = 0; sec < 4; sec++){ val = (Int_t)de->calselftrig[sec][6]; del = delay(val); - if ( del < 1100 ){ + if ( del < 1000 ){ clevel2->wartrig = 0.; clevel2->trigty = 3.; - ck = true; + // ck = true; break; }; }; - if ( !ck ) clevel2->wartrig = 100.; + // if ( !ck ) clevel2->wartrig = 100.; } else { Bool_t ck = false; for (Int_t sec = 0; sec < 4; sec++){ val = (Int_t)de->calselftrig[sec][6]; del = delay(val); - if ( del < 1100 ){ + if ( del < 1000 ){ clevel2->wartrig = 0.; ck = true; }; @@ -837,27 +1242,34 @@ Int_t se = 5; Int_t done = 0; Int_t pre = -1; - Bool_t isCOMP = false; - Bool_t isFULL = false; + // Bool_t isCOMP = false; + // Bool_t isFULL = false; Bool_t isRAW = false; Float_t ener; Int_t doneb = 0; Int_t donec = 0; - Int_t ck = 0; + Int_t ck[2][22][6]; + memset(ck, 0, 2*22*6*sizeof(Int_t)); Int_t ipre = 0; - Int_t ip[3] = {0}; + // Int_t ip[3] = {0}; + Int_t ip[3] = {0,0,0}; + Int_t ipp = 0; Float_t base0, base1, base2; base0 = 0.; base1 = 0.; base2 = 0.; - Float_t qpre[6] = {0.,0.,0.,0.,0.,0.}; + Float_t qpre[2][22][6]; + memset(qpre, 0, 2*22*6*sizeof(Float_t)); Float_t ene[96]; Int_t chdone[4] = {0,0,0,0}; Int_t pe = 0; // Float_t ener0 = 0.; Float_t cbase0 = 0.; + Float_t totbase = 0.; + Float_t totped = 0.; Bool_t pproblem = false; + Bool_t negbase = false; // Float_t tim = 0.; Int_t plo = 0; @@ -871,6 +1283,7 @@ // // determine the section number // + negbase = false; se = 5; if (l == 0 && m%2 == 0) se = 3; if (l == 0 && m%2 != 0) se = 2; @@ -879,11 +1292,11 @@ // // determine what kind of event we are going to analyze // - isCOMP = false; - isFULL = false; + // isCOMP = false; + // isFULL = false; isRAW = false; - if ( de->stwerr[se] & (1 << 16) ) isCOMP = true; - if ( de->stwerr[se] & (1 << 17) ) isFULL = true; + // if ( de->stwerr[se] & (1 << 16) ) isCOMP = true; + // if ( de->stwerr[se] & (1 << 17) ) isFULL = true; if ( de->stwerr[se] & (1 << 3) ) isRAW = true; if ( !chdone[se] ){ // @@ -896,7 +1309,7 @@ }; clevel2->perr[se] = 0; if ( de->perror[se] != 0 ){ - clevel2->perr[se] = 1; + clevel2->perr[se] = (Int_t)de->perror[se]; pe++; }; clevel2->swerr[se] = 0; @@ -912,7 +1325,7 @@ pre = -1; // for (Int_t nn = 0; nn < 96; nn++){ - ene[nn] = 0.; + // ene[nn] = 0.; dexy[l][m][nn] = de->dexy[l][m][nn] ; dexyc[l][m][nn] = de->dexyc[l][m][nn] ; }; @@ -921,44 +1334,66 @@ // pre = -1; cbase0 = 0.; + Int_t nstt[2]; + Float_t rqp[2]; for (Int_t i = 0; i < 3; i++){ + nstt[0] = 1000; + nstt[1] = 1000; + rqp[0] = 0.; + rqp[1] = 0.; for (Int_t j = 0; j < 2; j++){ pre = j + i*2; // // baseline check and calculation // if ( !isRAW ){ + // + // if it is a compress event with fully transmitted pre try to calculate the baseline + // if ( de->base[l][m][pre] != 0. && de->base[l][m][pre]<31000. ) { base[l][m][pre] = de->base[l][m][pre] ; } else { - FindBaseCompress(l,m,pre); + FindBaseCompress(l,m,pre,nstt[j],rqp[j]); }; cbase0 += base[l][m][pre]; } else { // - // if it is a raw event and we haven't checked - // yet, calculate the baseline. + // if it is a raw event calculate the baseline. // - FindBaseRaw(l,m,pre); + FindBaseRaw(l,m,pre,nstt[j],rqp[j]); cbase0 += base[l][m][pre]; }; }; + // + // if we are able to calculate the baseline with more than 3 strips on one pre and not in the other one choose the pre with more calculated strips + // + if ( nstt[0] < 4 && nstt[1] >= 4 && nstt[0] != 1000 && nstt[1] != 1000 ) base[l][m][pre-1] = 31000.; + if ( nstt[0] >= 4 && nstt[1] < 4 && nstt[0] != 1000 && nstt[1] != 1000 ) base[l][m][pre] = 31000.; + // // + // // if we are NOT able to calculate the baseline with more than 3 strips on both pres take the baseline (if any) of the one which has less energy + // // + // if ( nstt[0] < 4 && nstt[1] < 4 ){ + // if ( rqp[0] >= rqp[1] ) base[l][m][pre-1] = 31000.; + // if ( rqp[0] < rqp[1] ) base[l][m][pre] = 31000.; + // }; }; // // run over strips // pre = -1; ener0 = 0.; + totbase = 0.; + totped = 0.; for (Int_t i = 0 ; i < 3 ; i++){ ip[i] = 0; for (Int_t n = i*32 ; n < (i+1)*32 ; n++){ if (n%16 == 0) { - ck = 0; done = 0; doneb = 0; donec = 0; pre++; - qpre[pre] = 0.; + ck[l][m][pre] = 0; + qpre[l][m][pre] = 0.; }; // // baseline check and calculation @@ -966,8 +1401,9 @@ // no suitable new baseline, use old ones! // if ( !done ){ + if ( debug ) printf(" l %i m %i pre %i ip[i] %i base %f base ip[i] %f sbase %f \n",l,m,pre,ip[i],base[l][m][pre],base[l][m][ip[i]],sbase[l][m][pre]); if ( (base[l][m][pre] == 31000. || base[l][m][pre] == 0.) ){ - ck = 1; + ck[l][m][pre] = 1; if (pre%2 == 0) { ip[i] = pre + 1; } else { @@ -975,20 +1411,20 @@ }; if ( (base[l][m][ip[i]] == 31000. || base[l][m][ip[i]] == 0. || !crosst ) ){ // - ck = 2; + ck[l][m][pre] = 2; if ( sbase[l][m][pre] == 31000. || sbase[l][m][pre] == 0. ) { - ck = 3; + ck[l][m][pre] = 3; }; }; - done = 1; }; + done = 1; }; // // CALIBRATION ALGORITHM // if ( !doneb ){ - if ( debug ) printf(" ck is %i \n",ck); - switch (ck) { + if ( debug ) printf(" ck[l][m][pre] is %i \n",ck[l][m][pre]); + switch (ck[l][m][pre]) { case 0: base0 = base[l][m][pre]; base2 = calbase[l][m][pre]; @@ -1016,11 +1452,15 @@ ener = dexyc[l][m][n]; ener0 += ener; clevel1->estrip[n][m][l] = 0.; + totbase += de->base[l][m][pre]/96.; + totped += fabs(calped[l][m][n]); + if ( de->base[l][m][pre] < 0 ) negbase = true; if ( base0>0 && base0 < 30000. ){ - // if ( !donec && (base0 - base1 + base2) != 0. ){ - // sbase[l][m][pre] = base0 - base1 + base2; - if ( !donec && (base0 + base1 - base2) != 0. ){ - sbase[l][m][pre] = base0 + base1 - base2; + // + // save the baseline only if the energy release is "small" + // + if ( !donec && (base0 + base1 - base2) != 0. && (n+1)%16==0 ){ + if ( qpre[l][m][pre] < 200. ) sbase[l][m][pre] = base0 + base1 - base2; donec = 1; }; if ( ener > 0. ){ @@ -1028,157 +1468,410 @@ // // OK, now in estrip we have the energy deposit in MIP of all the strips for this event (at the end of loops of course) // - qpre[pre] += clevel1->estrip[n][m][l]; + if ( debug && l==0 && (m==17 || m==18) ) printf(" view %i plane %i strip %i ener %f calped %f base0 %f base1 %f base2 %f mip %f ENERGIA %f \n",l,m,n,ener,calped[l][m][n],base0,base1,base2,mip[l][m][n],clevel1->estrip[n][m][l]); + if ( clevel1->estrip[n][m][l] > 0. ) qpre[l][m][pre] += clevel1->estrip[n][m][l]; // // }; }; }; - if ( crosst ){ - if (ck == 1){ - if (ip[i]%2 == 0) { - ipre = ip[i] + 1; - } else { - ipre = ip[i] - 1; - }; - for (Int_t j = ipre*16 ; j < (ipre+1)*16 ; j++){ - if ( !ctground ){ - clevel1->estrip[j][m][l] += (qpre[ipre] - qpre[ip[i]]) * ctprecor[l][m][ipre]; - } else { - clevel1->estrip[j][m][l] += (qpre[ipre] - qpre[ip[i]]) * 0.00478; - }; - }; - }; - if (ck == 2){ - for (Int_t j = i*32 ; j < (i+1)*32 ; j++){ - ipre = j/16 + 1; - if ( !ctground ){ - clevel1->estrip[j][m][l] += qpre[ipre] * ctprecor[l][m][ipre]; - } else { - clevel1->estrip[j][m][l] += qpre[ipre] * 0.00478; - }; - }; - }; - }; }; // - if ( ener0 == 0. && cbase0 == 0. && !pproblem && clevel2->perr[se] == 0){ - if ( verbose ) printf(" L0 entry %i : calorimeter power problems! event marked as bad perr %f swerr %X view %i plane %i \n",ei,de->perror[se],de->stwerr[se],l,m); + // check if there were problems with 5.7 or glitches in the power supply + // + // if ( ((ener0 == 0. && cbase0 == 0.) || negbase || totbase > 196600. || totped < 1. ) && !pproblem && clevel2->perr[se] == 0){ // check pedestal and baseline values for one plane, if all zeros calibration is not valid (calorimeter power problems) [8th data reduction bug, fixed on 25/11/2009 by E.M.] + if ( ((ener0 == 0. && cbase0 == 0.) || negbase || totbase > 32700. || totped < 1. ) && !pproblem && clevel2->perr[se] == 0){ // check pedestal and baseline values for one plane, if all zeros calibration is not valid (calorimeter power problems) [8th data reduction bug, fixed on 25/11/2009 by E.M.] + if ( verbose ) printf(" L0 entry %i : calorimeter power problems! event marked as bad perr %f swerr %X view %i plane %i negbase %i totbase %f totped %f\n",ei,de->perror[se],de->stwerr[se],l,m, negbase, totbase, totped); pproblem = true; pe++; }; // - Int_t j4 = -4; - Int_t jjj = -3; - Int_t jj = -2; - Int_t jjpre = -1; - Int_t jjjpre = -1; - for (Int_t j = 0 ; j < 100 ; j++){ - jj++; - jjj++; - j4++; - if ( j < 96 ) ene[j] = clevel1->estrip[j][m][l]; - if ( crosst ){ - if ( jj >= 0 && jj < 96 ){ - if ( !ctground ){ - if ( jj%16 == 0 ) jjpre++; - if ( jj != 0 && jj != 32 && jj != 64 ) ene[jj-1] += -clevel1->estrip[jj][m][l] * ctneigcor[l][m][jjpre]; - if ( jj != 31 && jj != 63 && jj != 95 ) ene[jj+1] += -clevel1->estrip[jj][m][l] * ctneigcor[l][m][jjpre]; - } else { - if ( jj != 0 && jj != 32 && jj != 64 ) ene[jj-1] += -clevel1->estrip[jj][m][l] * 0.01581; - if ( jj != 31 && jj != 63 && jj != 95 ) ene[jj+1] += -clevel1->estrip[jj][m][l] * 0.01581; + } else { + for (Int_t nn = 0; nn < 96; nn++){ + clevel1->estrip[nn][m][l] = 0.; + }; + }; + }; + }; + // + // run over views and planes to apply crosstalk corrections + // + for (Int_t l = 0; l < 2; l++){ + for (Int_t m = 0; m < 22; m++){ + // + // determine the section number + // + se = 5; + if (l == 0 && m%2 == 0) se = 3; + if (l == 0 && m%2 != 0) se = 2; + if (l == 1 && m%2 != 0) se = 1; + if (l == 1 && m%2 == 0) se = 0; + // + // check for any error in the event + // + if ( clevel2->crc[se] == 0 && (clevel1->good2 == 1 || clevel2->trigty >= 2) ){ + // + // Cross-talk corrections + // + if ( crosst ){ + // + // energy on silicon ladders + // + Float_t qsi[3]; + qsi[0] = qpre[l][m][0]+qpre[l][m][1]; + qsi[1] = qpre[l][m][2]+qpre[l][m][3]; + qsi[2] = qpre[l][m][4]+qpre[l][m][5]; + // + for ( pre = 1; pre < 6; pre += 2 ){ + Int_t ladder = (pre - 1)/2; + // + // If the noselfct flag is set the strip doesn't suffer the self crosstalk due to electronics so we must subtract some energy + // + if ( noselfct ){ + for (Int_t j = ladder*32 ; j < (ladder+1)*32 ; j++){ + ipre = j/16 ; + if ( debug ) printf(" CT STEP1 %i %i %i estrip %f ctprecor %f \n",j,m,l,clevel1->estrip[j][m][l],ctprecor[l][m][ipre]); + if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] -= clevel1->estrip[j][m][l] * ctprecor[l][m][ipre]; + if ( debug ) printf(" CT STEP2 %i %i %i estrip %f ctprecor %f \n",j,m,l,clevel1->estrip[j][m][l],ctprecor[l][m][ipre]); }; }; - if ( jjj >= 0 && jjj < 96 ){ - if ( !ctground ){ - if ( jjj%16 == 0 ) jjjpre++; - if ( jjj != 0 && jjj != 32 && jjj != 64 ) clevel1->estrip[jjj-1][m][l] += -ene[jjj] * ctneigcor[l][m][jjjpre]; - if ( jjj != 31 && jjj != 63 && jjj != 95 ) clevel1->estrip[jjj+1][m][l] += -ene[jjj] * ctneigcor[l][m][jjjpre]; + // + // Using the neighbour pre baseline + // + if (ck[l][m][pre] == 1 || ck[l][m][pre-1] == 1){ + // + // pre-amplifier effect on baseline when using the neighbour pre (ck=1) + // + if (ck[l][m][pre] == 1){ + ipre = pre; + ipp = pre - 1; } else { - if ( jjj != 0 && jjj != 32 && jjj != 64 ) clevel1->estrip[jjj-1][m][l] += -ene[jjj] * 0.01581; - if ( jjj != 31 && jjj != 63 && jjj != 95 ) clevel1->estrip[jjj+1][m][l] += -ene[jjj] * 0.01581; + ipre = pre - 1; + ipp = pre; + }; + Int_t it = 0; + Float_t nqpre = 0.; + // + if ( debug ) printf(" CK1 Limit for while: 0.07 \n"); + for (Int_t j = ipre*16 ; j < (ipre+1)*16 ; j++){ + if ( debug ) printf(" CT STEP3 %i %i %i estrip %f ctprecor %f \n",j,m,l,clevel1->estrip[j][m][l],ctprecor[l][m][ipp]); + if ( !ctground ){ + if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += - qpre[l][m][ipp] * ctprecor[l][m][ipp]; + } else { + if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += - qpre[l][m][ipp] * 0.00478; + }; + if ( clevel1->estrip[j][m][l] > 0. ) nqpre += clevel1->estrip[j][m][l] ; + if ( debug ) printf(" CT STEP4 %i %i %i estrip %f ctprecor %f \n",j,m,l,clevel1->estrip[j][m][l],ctprecor[l][m][ipp]); }; + qpre[l][m][ipre] = nqpre; + nqpre = 0.; + Float_t deltaqpre = qpre[l][m][ipre]; + // + // these values are empirically determined, usually the routine converge due to deltaqsi and the latest applied correction is based on less than 1 mip + // + while ( it < 10 && deltaqpre > 0.07 ){ + nqpre = 0.; + for (Int_t j = ipre*16 ; j < (ipre+1)*16 ; j++){ + if ( debug ) printf(" CT STEP5 %i %i %i estrip %f ctprecor %f \n",j,m,l,clevel1->estrip[j][m][l],ctprecor[l][m][ipre]); + if ( !ctground ){ + if ( debug ) printf(" CK1 pre correction: iteration %i deltaqpre %f ctprecor %f TOTAL CORRECTION %f \n",it,deltaqpre,ctprecor[l][m][ipre],deltaqpre * ctprecor[l][m][ipre]); + if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += deltaqpre * ctprecor[l][m][ipre]; + } else { + if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += deltaqpre * 0.00478; + }; + if ( clevel1->estrip[j][m][l] > 0. ) nqpre += clevel1->estrip[j][m][l] ; + if ( debug ) printf(" CT STEP6 %i %i %i estrip %f ctprecor %f \n",j,m,l,clevel1->estrip[j][m][l],ctprecor[l][m][ipre]); + }; + if ( ctground ) it = 100; + it++; + deltaqpre = nqpre - qpre[l][m][ipre]; + if ( debug ) printf(" CK1 BEFORE: qpre %f \n",qpre[l][m][ipre]); + qpre[l][m][ipre] = nqpre; + if ( debug ) printf(" CK1 AFTER: qpre %f \n",qpre[l][m][ipre]); + }; + // }; - }; - if ( j4 >= 0 && j4 < 96 ){ // - // NOTICE: THE FOLLOWING LINE EXCLUDE ALL STRIPS FOR WHICH THE RMS*4 IS GREATER THAN 26 !!! <=============== IMPORTANT! =================> + // No baseline calculation due to high energy release // - if ( obadmask[l][m][j4] == 1 || clevel1->estrip[j4][m][l] <= clevel1->emin || calrms[l][m][j4] > 26 ){ - clevel1->estrip[j4][m][l] = 0.; - }; - // - // code and save the energy for each strip in svstrip - // - if ( clevel1->estrip[j4][m][l] > clevel1->emin ){ + if (ck[l][m][pre] == 2 && ck[l][m][pre-1] == 2){ // - Float_t savel1 = clevel1->estrip[j4][m][l]; - if ( dexyc[l][m][j4] == 32767. ){ - savel1 += 5000.; - clevel2->nsatstrip += 1.; - }; + // y^ + // | + // | 6 7 8 + // | 3 4 5 + // | 0 1 2 + // | --------------------------------------> x // - tim = 100000.; - plo = m; - fbi = 0; - if ( savel1 > 0.99995 ){ - tim = 10000.; - plo = m; - fbi = 1; - }; - if ( savel1 > 9.9995 ){ - tim = 1000.; - plo = 22 + m; - fbi = 1; - }; - if ( savel1 > 99.995 ){ - tim = 100.; - plo = 22 + m; - fbi = 0; + Int_t si1 = 0; + Int_t si2 = 0; + Int_t si3 = 0; + if ( l == 0 ){ + if ( ladder == 0 ){ + si1 = 0; + si2 = 3; + si3 = 6; + }; + if ( ladder == 1 ){ + si1 = 1; + si2 = 4; + si3 = 7; + }; + if ( ladder == 2 ){ + si1 = 2; + si2 = 5; + si3 = 8; + }; + } else { + if ( ladder == 0 ){ + si1 = 0; + si2 = 1; + si3 = 2; + }; + if ( ladder == 1 ){ + si1 = 3; + si2 = 4; + si3 = 5; + }; + if ( ladder == 2 ){ + si1 = 6; + si2 = 7; + si3 = 8; + }; }; - if ( savel1 > 999.95 ){ - tim = 10.; - plo = 44 + m; - fbi = 0; + // + // Find the energy distribution along the considered plane looking at the two sandwiching plane of the other view. + // + Float_t sied[3] = {0.,0.,0.}; + Int_t othv = !l; + Int_t othpl1 = m - 1; + Int_t othpl2 = m + 1; + Float_t oprof[3] = {0.,0.,0.}; + for(Int_t s=0; s<3; s++){ + for(Int_t t=(s*32); t<32*(s + 1); t++){ + if ( othpl1 > -1 ) { + oprof[s] += clevel1->estrip[othv][othpl1][t]; + }; + if ( othpl2 < 22 ) { + oprof[s] += clevel1->estrip[othv][othpl2][t]; + }; + }; }; - if ( savel1 > 9999.5 ){ - tim = 1.; - plo = 66 + m; - fbi = 0; + Float_t otote = fabs(oprof[0]) + fabs(oprof[1]) + fabs(oprof[2]); + for(Int_t g=0; g<3; g++){ + if ( otote > 0. ){ + sied[g] = fabs(oprof[g])/otote; + } else { + sied[g] = 1./3.; + }; }; // - cle = (Int_t)lroundf(tim*savel1); // - if ( l == 0 ){ + // + Int_t it = 0; + Int_t jpre = 0; + Float_t nqsi = 0.; + Float_t snqsi = qsi[ladder]; + Float_t nqpre[2] = {0.,0.}; + Float_t deltaqsi = qsi[ladder]; + Float_t deltaqpre[2]; + deltaqpre[0] = qpre[l][m][pre-1]; + deltaqpre[1] = qpre[l][m][pre]; + // + if ( debug ) printf(" Limit for while: 0.07 it < 10 \n"); + // + // these values are empirically determined, usually the routine converge due to deltaqsi and the latest applied correction is based on less than 1 mip + // + while ( it < 10 && (deltaqsi > 0.07 || deltaqpre[0] > 0.07 || deltaqpre[1] > 0.07) ){ + nqsi = 0.; + nqpre[0] = 0.; + nqpre[1] = 0.; + for (Int_t j = ladder*32 ; j < (ladder+1)*32 ; j++){ + if ( debug ) printf(" CT STEP6 %i %i %i estrip %f ctprecor %f \n",j,m,l,clevel1->estrip[j][m][l],ctsicor[l][m][si2]); + ipre = 0; + if ( j > (ladder*32)+15 ) ipre = 1; + jpre = j/16 ; + // + // Silicon effect on the baseline when using the same pre previous baseline (ck = 2) + pre-amply effect + // + if ( !ctground ){ + if ( debug ) printf(" silicon correction: iteration %i deltaqsi[%i] %f ctsicor %f %f %f sied %f %f %f si %i %i %i TOTAL CORRECTION %f \n",it,ladder,deltaqsi,ctsicor[l][m][si1],ctsicor[l][m][si2],ctsicor[l][m][si3],sied[0],sied[1],sied[2],si1,si2,si3,deltaqsi * (ctsicor[l][m][si1] * sied[0] + ctsicor[l][m][si2] * sied[1] + ctsicor[l][m][si3] * sied[2])); + if ( debug ) printf(" pre correction: iteration %i deltaqpre[0] %f deltaqpre[1] %f ctprecor %f TOTAL CORRECTION %f \n",it,deltaqpre[0],deltaqpre[1],ctprecor[l][m][jpre],deltaqpre[ipre] * ctprecor[l][m][jpre]); + if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += (deltaqsi * (ctsicor[l][m][si1] * sied[0] + ctsicor[l][m][si2] * sied[1] + ctsicor[l][m][si3] * sied[2])/mip[l][m][j]) + deltaqpre[ipre] * ctprecor[l][m][jpre]; + } else { + if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += 0. + qpre[l][m][jpre] * 0.00478; // no correction + }; + if ( clevel1->estrip[j][m][l] > 0. ) nqsi += clevel1->estrip[j][m][l] ; + if ( clevel1->estrip[j][m][l] > 0. ) nqpre[ipre] += clevel1->estrip[j][m][l] ; + if ( debug ) printf(" CT STEP7 %i %i %i estrip %f ctprecor %f \n",j,m,l,clevel1->estrip[j][m][l],ctsicor[l][m][si2]); + }; + if ( ctground ) it = 100; + deltaqsi = nqsi-snqsi; + deltaqpre[0] = nqpre[0] - qpre[l][m][pre-1]; + deltaqpre[1] = nqpre[1] - qpre[l][m][pre]; // - // +-PPSSmmmm.mmmm + // Check for divergence and stop if it happens! [9RED bug noticed with plane 18X] // - svstrip[istrip] = fbi*1000000000 + plo*10000000 + j4*100000 + cle; - } else { - svstrip[istrip] = -(fbi*1000000000 + plo*10000000 + j4*100000 + cle); + if ( deltaqpre[0] > qpre[l][m][pre-1] || deltaqpre[1] > qpre[l][m][pre] || deltaqsi >snqsi ){ + if ( debug ) printf(" WARNING!! DIVERGING CORRECTION EXIT IMMEDIATLY FROM THE LOOP!! dqpre0 %f qpre0 %f // dqpre1 %f qpre1 %f // dqsi %f qsi %f \n",deltaqpre[0],qpre[l][m][pre-1],deltaqpre[1],qpre[l][m][pre],deltaqsi,snqsi); + it = 1000; + }; + // + snqsi = nqsi; + it++; + if ( debug ) printf(" BEFORE: qpre 0 %f qpre 1 %f \n",qpre[l][m][pre-1],qpre[l][m][pre]); + qpre[l][m][pre-1] = nqpre[0]; + qpre[l][m][pre] = nqpre[1]; + if ( debug ) printf(" AFTER: qpre 0 %f qpre 1 %f \n",qpre[l][m][pre-1],qpre[l][m][pre]); + }; + // + // + // +// for (Int_t j = ladder*32 ; j < (ladder+1)*32 ; j++){ +// ipre = j/16 ; +// // +// // pre-amplifier effect on baseline when using the same pre previous event baseline (ck=2) +// // +// if ( !ctground ){ +// if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += qpre[l][m][ipre] * ctprecor[l][m][ipre]; +// } else { +// if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += qpre[l][m][ipre] * 0.00478; +// }; +// }; + }; + }; + }; + }; + // + Int_t j4 = -4; + Int_t jjj = -3; + Int_t jj = -2; + Int_t jjpre = -1; + Int_t jjjpre = -1; + memset(ene, 0, 96*sizeof(Float_t)); + for (Int_t j = 0 ; j < 100 ; j++){ + jj++; + jjj++; + j4++; + if ( j < 96 ) ene[j] = clevel1->estrip[j][m][l]; + if ( crosst ){ + // + // "Real" crosstalk effect on the neighbour strips respect to the one which have seen the energy deposit + // + if ( jj >= 0 && jj < 96 ){ + if ( !ctground ){ + if ( jj%16 == 0 ) jjpre++; + if ( jj != 0 && jj != 32 && jj != 64 && ene[jj-1] != 0. ) ene[jj-1] += -clevel1->estrip[jj][m][l] * ctneigcor[l][m][jjpre]; + if ( jj != 31 && jj != 63 && jj != 95 && ene[jj+1] != 0. ) ene[jj+1] += -clevel1->estrip[jj][m][l] * ctneigcor[l][m][jjpre]; + } else { + if ( jj != 0 && jj != 32 && jj != 64 && ene[jj-1] != 0. ) ene[jj-1] += -clevel1->estrip[jj][m][l] * 0.01581; + if ( jj != 31 && jj != 63 && jj != 95 && ene[jj+1] != 0. ) ene[jj+1] += -clevel1->estrip[jj][m][l] * 0.01581; + }; + }; + if ( jjj >= 0 && jjj < 96 ){ + if ( !ctground ){ + if ( jjj%16 == 0 ) jjjpre++; + if ( jjj != 0 && jjj != 32 && jjj != 64 && clevel1->estrip[jjj-1][m][l] != 0. ) clevel1->estrip[jjj-1][m][l] += -ene[jjj] * ctneigcor[l][m][jjjpre]; + if ( jjj != 31 && jjj != 63 && jjj != 95 && clevel1->estrip[jjj+1][m][l] !=0. ) clevel1->estrip[jjj+1][m][l] += -ene[jjj] * ctneigcor[l][m][jjjpre]; + } else { + if ( jjj != 0 && jjj != 32 && jjj != 64 && clevel1->estrip[jjj-1][m][l] != 0. ) clevel1->estrip[jjj-1][m][l] += -ene[jjj] * 0.01581; + if ( jjj != 31 && jjj != 63 && jjj != 95 && clevel1->estrip[jjj+1][m][l] != 0. ) clevel1->estrip[jjj+1][m][l] += -ene[jjj] * 0.01581; + }; + }; + }; + if ( j4 >= 0 && j4 < 96 ){ + // + // CALOLEVEL1 CODING AND FILLING + // + // + // NOTICE: THE FOLLOWING LINE EXCLUDE ALL STRIPS FOR WHICH THE RMS*4 IS GREATER THAN 26 !!! <=============== IMPORTANT! =================> // not true anymore, now it trust parameter files + // + if ( obadmask[l][m][j4] == 1 || clevel1->estrip[j4][m][l] <= clevel1->emin || clevel1->estrip[j4][m][l] <= memin[l][m][j4] || calrms[l][m][j4] > maxrms[l][m] || (l==0 && m == 18 && mask18 ) ){ + clevel1->estrip[j4][m][l] = 0.; + }; + // + if ( debug ) printf(" STRIP: view %i plane %i strip %i energy: %f \n",l,m,j4,clevel1->estrip[j4][m][l]); + // + // code and save the energy for each strip in svstrip + // + if ( clevel1->estrip[j4][m][l] > clevel1->emin ){ + // + Float_t savel1 = clevel1->estrip[j4][m][l]; + // + if ( m == 18 && l == 0 ){ + if ( debug ) printf(" Resetting plane 18X for variable calculation: view %i plane %i strip %i \n",l,m,j4); + clevel1->estrip[j4][m][l] = 0.; // SAVE STRIPS VALUE FOR PLANE 18 X but DO NOT USE IT FOR VARIABLE CALCULATION + }; + if ( debug ) printf(" HIT STRIP: view %i plane %i strip %i energy: %f \n",l,m,j4,clevel1->estrip[j4][m][l]); + // if ( dexyc[l][m][j4] == 32767. ){ + if ( dexyc[l][m][j4] > 32000. || savel1 > 5000.){ // CaloLevel1 bug with plane 18X [9RED 14/04/2010] + if ( savel1 > 5000 ){ + if ( debug ) printf(" Absurd plane 18X energy... resetting value to 1100 MIP \n"); + savel1 = 1100.; // CaloLevel1 bug with plane 18x [9RED 14/04/2010] }; + savel1 += 5000.; + clevel2->nsatstrip += 1.; + }; + // + tim = 100000.; + plo = m; + fbi = 0; + if ( savel1 > 0.99995 ){ + tim = 10000.; + plo = m; + fbi = 1; + }; + if ( savel1 > 9.9995 ){ + tim = 1000.; + plo = 22 + m; + fbi = 1; + }; + if ( savel1 > 99.995 ){ + tim = 100.; + plo = 22 + m; + fbi = 0; + }; + if ( savel1 > 999.95 ){ + tim = 10.; + plo = 44 + m; + fbi = 0; + }; + if ( savel1 > 9999.5 ){ + tim = 1.; + plo = 66 + m; + fbi = 0; + }; + // + cle = (Int_t)lroundf(tim*savel1); + // + if ( l == 0 ){ // - // if ( ei >= -770 ) printf(" j %i l %i m %i estrip %f \n",j4,l,m,clevel1->estrip[j4][m][l]); - // if ( ei >= -770 ) printf(" num lim %i fbi %i tim %f plo %i cle %i \n",numeric_limits::max(),fbi,tim,plo,cle); - // if ( ei >= -770 ) printf(" svstrip %i \n",svstrip[istrip]); + // +-PPSSmmmm.mmmm // - istrip++; + svstrip[istrip] = fbi*1000000000 + plo*10000000 + j4*100000 + cle; + } else { + svstrip[istrip] = -(fbi*1000000000 + plo*10000000 + j4*100000 + cle); }; + if ( debug ) printf(" svstrip[%i] = %i fbi %i plo %i j4 %i cle %i \n",istrip,svstrip[istrip],fbi,plo,j4,cle); + // + istrip++; }; }; - // - } else { - for (Int_t nn = 0; nn < 96; nn++){ - clevel1->estrip[nn][m][l] = 0.; - }; }; + // }; }; + // + // store goodness flag + // if ( !pe ){ clevel2->good = 1; } else { clevel2->good = 0; }; + // + // done + // return(0); } @@ -1256,13 +1949,13 @@ void CaloLevel0::FillCommonVar(CaloLevel1 *c1, CaloLevel2 *ca){ // ca->good = clevel2->good; - if ( clevel2->trigty == 2. ){ - ca->selftrigger = 1; - } else { - ca->selftrigger = 0; - }; +// if ( clevel2->trigty == 2. ){ +// ca->selftrigger = 1; +// } else { +// ca->selftrigger = 0; +// }; // - ca->selftrigger += (Int_t)clevel2->wartrig; + ca->selftrigger = (Int_t)clevel2->trigty + (Int_t)clevel2->wartrig; // memcpy(ca->perr,clevel2->perr,sizeof(clevel2->perr)); memcpy(ca->swerr,clevel2->swerr,sizeof(clevel2->swerr)); @@ -1299,12 +1992,25 @@ }; // } +void CaloLevel0::FillCommonVar(CaloLevel1 *c1){ + if ( c1 ){ + c1->istrip = istrip; + c1->estrip = TArrayI(istrip,svstrip); + }; + // +} void CaloLevel0::ClearStructs(){ ClearTrkVar(); ClearCommonVar(); } +void CaloLevel0::Delete(Option_t *t){ + if ( de ) delete de; + delete this; +} + + void CaloLevel0::RunClose(){ l0tr->Delete(); ClearStructs(); @@ -1314,6 +2020,7 @@ memset(base, 0, 2*22*6*sizeof(Float_t)); memset(sbase, 0, 2*22*6*sizeof(Float_t)); memset(ctprecor, 0, 2*22*6*sizeof(Float_t)); + memset(ctsicor, 0, 2*22*9*sizeof(Float_t)); memset(ctneigcor, 0, 2*22*6*sizeof(Float_t)); // } @@ -1532,7 +2239,7 @@ // if ( !ncalibs ) return(-110); // - calo->GetEntry(calibno[s]); + if ( calo->GetEntry(calibno[s]) <= 0 ) throw -36; // if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) { for ( Int_t d=0 ; d<11 ;d++ ){