loon: LILinResp.cxx Source File

00001 #define LILinResp_cxx
00002 
00003 #include <LILinResp.h>
00004 #include <TF1.h>
00005 
00006 using namespace std;
00007 
00008 std::vector<Float_t> LILinResp::fXrefPINLG1;
00009 std::vector<Float_t> LILinResp::fXrefPINLG2;
00010 std::vector<Float_t> LILinResp::fXrefPINLG3;
00011 std::vector<Float_t> LILinResp::fXrefPINLG4;
00012 std::vector<Float_t> LILinResp::fXrefPINLG5;
00013 std::vector<Float_t> LILinResp::fXrefPINLG6;
00014 
00015 std::vector<Float_t> LILinResp::fXrefPINHG1;
00016 std::vector<Float_t> LILinResp::fXrefPINHG2;
00017 std::vector<Float_t> LILinResp::fXrefPINHG3;
00018 std::vector<Float_t> LILinResp::fXrefPINHG4;
00019 std::vector<Float_t> LILinResp::fXrefPINHG5;
00020 std::vector<Float_t> LILinResp::fXrefPINHG6;
00021 
00022 CalVaLinearity* LILinResp::myCV = 0;
00023 
00024 
00025 CVSID("$Id: LILinResp.cxx,v 1.4 2005/04/18 14:24:54 anatael Exp $");
00026 
00027 //-------------------------------------------------------------
00028 
00029 LILinResp::LILinResp() { 
00030 
00031   MSG("LILinResp",Msg::kDebug) << "Making a LILinResp..." << endl;
00032 
00033   entryHG = 0;
00034   entryLG = 0;
00035 
00036   addressKey  = -1;
00037 
00038   slopeHG       = -1;
00039   interceptHG   = -1;
00040   slopeHG_e     = -1;
00041   interceptHG_e = -1;
00042 
00043   slopeLG       = -1;
00044   interceptLG   = -1;
00045   slopeLG_e     = -1;
00046   interceptLG_e = -1;
00047 
00048   gain  = -1;
00049   drift = -1;
00050   led   = -1;
00051 
00052   //Zeroing all pointers:
00053   resPlotHG_g = 0;
00054   resPlotLG_g = 0;
00055 }
00056 
00057 //-------------------------------------------------------------
00058 
00059 LILinResp::~LILinResp() {
00060 
00061   if(resPlotHG_g) delete resPlotHG_g; 
00062   if(resPlotLG_g) delete resPlotLG_g; 
00063 
00064   delete myCV;
00065   myCV = 0;
00066 
00067   MSG("LILinResp",Msg::kInfo) << "Killing a LILinResp" << endl;
00068 }
00069 
00070 //-------------------------------------------------------------
00071 
00072 Bool_t LILinResp::GetLinearity() {
00073 
00074   Bool_t success = true;
00075 
00076   //
00077   // Create VA-nonLinearity object -> linearise PINs, which are FD electronis:
00078   //
00079   if(!myCV) myCV = new CalVaLinearity( RawChannelId(),
00080                                        24,
00081                                        16061.3,
00082                                          4.63,
00083                                        13683.9,
00084                                        6650 );
00085 
00086   //
00087   // Fill TGraph with ALL the information
00088   // for the linearity calculation:
00089   //
00090   vector<Float_t>::iterator xIter  = fX.begin();
00091   vector<Float_t>::iterator xeIter = fX_error.begin();
00092 
00093   vector<Float_t>::iterator xRefIterLG;// = NULL;
00094   vector<Float_t>::iterator xRefIterHG;// = NULL;
00095   
00096   //Use the LED number to ID the corresponding PIN (reference X):
00097   switch (led) {
00098   case 1:
00099     xRefIterLG = fXrefPINLG1.begin();
00100     xRefIterHG = fXrefPINHG1.begin();
00101     break;
00102   case 2:
00103     xRefIterLG = fXrefPINLG2.begin();
00104     xRefIterHG = fXrefPINHG2.begin();
00105     break;
00106   case 3:
00107     xRefIterLG = fXrefPINLG3.begin();
00108     xRefIterHG = fXrefPINHG3.begin();
00109     break;
00110   case 4:
00111     xRefIterLG = fXrefPINLG4.begin();
00112     xRefIterHG = fXrefPINHG4.begin();
00113     break;
00114   case 5:
00115     xRefIterLG = fXrefPINLG5.begin();
00116     xRefIterHG = fXrefPINHG5.begin();
00117     break;
00118   case 6:
00119     xRefIterLG = fXrefPINLG6.begin();
00120     xRefIterHG = fXrefPINHG6.begin();
00121     break;
00122   }
00123 
00124   //
00125   //LOOP over vectors:
00126   //
00127   while( xIter != fX.end() ) {
00128 
00129     //Avoid using entries with -1's
00130     if((*xRefIterHG) > 0  && (*xIter) > 2) {
00131       
00132       success = this->FillGraphHG( myCV->Linearize( (*xRefIterHG) ),
00133                                    0.,
00134                                    (*xIter),
00135                                    (*xeIter) );  
00136     }
00137     
00138     if((*xRefIterLG) > 0 && (*xIter) > 2) {
00139       
00140       success = this->FillGraphLG( myCV->Linearize( (*xRefIterLG) ),
00141                                    0., 
00142                                    (*xIter),
00143                                    (*xeIter) );  
00144     }
00145 
00146     xRefIterHG++;
00147     xRefIterLG++;
00148     xIter++;
00149     xeIter++;
00150   }    
00151 
00152   //
00153   //Get linearity:
00154   //
00155 
00156   TF1* lin = new TF1("lin","[0]+[1]*x"); //"pol1");
00157 
00158 
00159   lin->FixParameter(0,0);
00160   //lin->SetParLimits(0,0.001,-0.001);
00161   lin->SetParameter(1,0.3);
00162 
00163 
00164   Int_t out = -1;
00165 
00166   //This complication is needed to tune the fit to the lowest possible region, which is LED dependent
00167   if(typeFEE == ElecType::kVA) {
00168     switch( led ) {
00169     case 1:
00170       //if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,4000);
00171       break;
00172     case 2:
00173       //if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,2200);
00174       break;
00175     case 3:
00176       //if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,5500);
00177       break;
00178     case 4:
00179       if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,3800);
00180       break;
00181     case 5:
00182       if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,5000);
00183       break;
00184     case 6:
00185       if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,3000);
00186       break;
00187     }
00188   }
00189   else {
00190     switch( led ) {
00191     case 1:
00192       if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,4000);
00193       break;
00194     case 2:
00195       if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,2200);
00196       break;
00197     case 3:
00198       if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,5500);
00199       break;
00200     case 4:
00201       //if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,3500);
00202       break;
00203     case 5:
00204       //if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,5000);
00205       break;
00206     case 6:
00207       //if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,4000);
00208       break;
00209     }
00210   }
00211 
00212   slopeHG       = lin->GetParameter(1);
00213   slopeHG_e     = lin->GetParError(1);
00214   interceptHG   = lin->GetParameter(0);
00215   interceptHG_e = lin->GetParError(0);
00216 
00217 
00218   lin->FixParameter(0,0);
00219   lin->SetParameter(1,0.3);
00220   
00221 
00222   //This complication is needed to tune the fit to the lowest possible region, which is LED dependent
00223   if(typeFEE == ElecType::kVA) {
00224     switch( led ) {
00225     case 1:
00226       //if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,500);
00227       break;
00228     case 2:
00229       //if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,500);
00230       break;
00231     case 3:
00232       //if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,500);
00233       break;
00234     case 4:
00235       if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,800);
00236       break;
00237     case 5:
00238       if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,1000);
00239       break;
00240     case 6:
00241       if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,600);
00242       break;
00243     }
00244   }
00245   else {
00246     switch( led ) {
00247     case 1:
00248       if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,800);//1000);
00249       break;
00250     case 2:
00251       if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,800);//1100);
00252       break;
00253     case 3:
00254       if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,800);//1000,2300); //Structure at low charges
00255       break;
00256     case 4:
00257       if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,800);//900);
00258       break;
00259     case 5:
00260       if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,800);//1200);
00261       break;
00262     case 6:
00263       if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,800);//1200);
00264       break;
00265     }
00266   }
00267 
00268 
00269   //Check quality of LG, HG info is not used
00270   if(slopeLG==0) return false;
00271   if(slopeLG>0.29) return false;
00272 
00273 
00274   slopeLG       = lin->GetParameter(1);
00275   slopeLG_e     = lin->GetParError(1);
00276   interceptLG   = lin->GetParameter(0);
00277   interceptLG_e = lin->GetParError(0);
00278 
00279   delete lin;
00280 
00281 
00282 
00283   xIter  = fX.begin();
00284   //xeIter = fX_error.begin();
00285 
00286   //vector<Float_t>::iterator xRefIterLG;
00287   //vector<Float_t>::iterator xRefIterHG;
00288   
00289   //Use the LED number to ID the corresponding PIN (reference X):
00290   switch (led) {
00291   case 1:
00292     xRefIterLG = fXrefPINLG1.begin();
00293     xRefIterHG = fXrefPINHG1.begin();
00294     break;
00295   case 2:
00296     xRefIterLG = fXrefPINLG2.begin();
00297     xRefIterHG = fXrefPINHG2.begin();
00298     break;
00299   case 3:
00300     xRefIterLG = fXrefPINLG3.begin();
00301     xRefIterHG = fXrefPINHG3.begin();
00302     break;
00303   case 4:
00304     xRefIterLG = fXrefPINLG4.begin();
00305     xRefIterHG = fXrefPINHG4.begin();
00306     break;
00307   case 5:
00308     xRefIterLG = fXrefPINLG5.begin();
00309     xRefIterHG = fXrefPINHG5.begin();
00310     break;
00311   case 6:
00312     xRefIterLG = fXrefPINLG6.begin();
00313     xRefIterHG = fXrefPINHG6.begin();
00314     break;
00315   }
00316 
00317   //
00318   //LOOP over the responses of each LILinResp object:
00319   //
00320   while( xIter != fX.end() ) {
00321     
00322     //If there was bad data (ID by "-1") -> label as bad data.
00323     if( (*xRefIterLG)==-1 || (*xRefIterHG)==-1 || (*xIter)==-1 ) {
00324       fResidual.push_back( -1 );
00325       fResidualPMT.push_back( -1 );
00326     }
00327     else {
00328 
00329       //Furthermore, HG suffers from non-linearity of VA: BEWARE 
00330       //Float_t func = slopeHG * (*xRefIterHG) + interceptHG; 
00331       Float_t func = slopeLG * (*xRefIterLG) + interceptLG;
00332       
00333       //Calculating residuals:
00334       Float_t res    = (*xIter) - func;
00335       //% residual
00336       //Normalise by LIN:
00337       if( func != 0 ) res /= func;
00338       //Normalise by NON-LIN:
00339       //if( (*xIter)!= 0 ) res /= (*xIter);
00340       else res = -1;
00341 
00342       //Get residuals into vector of residuals:
00343       fResidual.push_back( res );
00344 
00345 
00346       //Remove linearity from PMT response on FD side:
00347       if(typeFEE == ElecType::kVA) {
00348 
00349         const Float_t linResp = myCV->Linearize( (*xIter)*gain )/gain;
00350         //Calculating residuals:
00351         Float_t resPMT = linResp - func;
00352 
00353         //Normalise by LIN:
00354         if( func != 0 ) resPMT /= func;
00355         //Normalise by Non-LIN:
00356         //if( linResp != 0 ) resPMT /= linResp;
00357         else resPMT = -1;
00358         fResidualPMT.push_back( resPMT );
00359 
00360       }      
00361       else 
00362         fResidualPMT.push_back( -1 ); //For ND.
00363 
00364     }
00365     
00366     xIter++;
00367     xRefIterHG++;
00368     xRefIterLG++;
00369   }
00370 
00371   return success;
00372 }
00373 
00374 //--------------------------------------------------------
00375 
00376 Bool_t LILinResp::FillGraphHG( const Float_t x,
00377                              const Float_t xe,
00378                              const Float_t y,
00379                              const Float_t ye ) {
00380 
00381   //Check whether created or not:
00382   if(!resPlotHG_g) {
00383 
00384     MSG("LILinResp",Msg::kVerbose) << "Getting TGraph..." << endl;
00385 
00386     resPlotHG_g = new TGraphErrors();
00387 
00388     //Setting Style of Graph:
00389     resPlotHG_g->SetLineWidth(4);
00390     resPlotHG_g->SetLineColor(1);
00391     resPlotHG_g->SetMarkerSize(1);
00392     resPlotHG_g->SetMarkerStyle(20);
00393     resPlotHG_g->SetMarkerColor(2);
00394   }
00395 
00396   if(resPlotHG_g) {
00397 
00398     resPlotHG_g->SetPoint(entryHG,x,y);
00399     resPlotHG_g->SetPointError(entryHG,xe,ye);
00400 
00401     entryHG++;
00402     return true;
00403   }
00404 
00405   else return false;
00406 }
00407 
00408 //--------------------------------------------------------
00409 
00410 Bool_t LILinResp::FillGraphLG( const Float_t x,
00411                              const Float_t xe,
00412                              const Float_t y,
00413                              const Float_t ye ) {
00414 
00415   //Check whether created or not:
00416   if(!resPlotLG_g) {
00417 
00418     MSG("LILinResp",Msg::kVerbose) << "Getting TGraph..." << endl;
00419 
00420     resPlotLG_g = new TGraphErrors();
00421 
00422     //Setting Style of Graph:
00423     resPlotLG_g->SetLineWidth(4);
00424     resPlotLG_g->SetLineColor(1);
00425     resPlotLG_g->SetMarkerSize(1);
00426     resPlotLG_g->SetMarkerStyle(20);
00427     resPlotLG_g->SetMarkerColor(2);
00428   }
00429 
00430   if(resPlotLG_g) {
00431 
00432     resPlotLG_g->SetPoint(entryLG,x,y);
00433     resPlotLG_g->SetPointError(entryLG,xe,ye);
00434 
00435     entryLG++;
00436     //resPlotLG_g->Print();
00437     return true;
00438   }
00439 
00440   else return false;
00441 }
00442 
00443 //--------------------------------------------------------
00444 
00445 void LILinResp::PrintMe() {
00446 
00447   //cout << gain << "\t" << drift << "\t" << led << "\t" << addressKey << "\t" << entryHG << "\t" << entryLG << "\t" << typeFEE << endl;
00448 
00449   //resPlotHG_g->Print();
00450   //resPlotLG_g->Print();
00451 
00452   //cout << "X:    " << fX.size() << " (" << (&fX) << ")" << endl;
00453   //cout << "X_e:  " << fX_error.size()  << " (" << (&fX_error)  << ")" << endl;
00454   //cout << "Xres: " << fResidual.size() << " (" << (&fResidual) << ")" << endl;
00455 
00456   //cout << LILinResp::fXrefPINHG1.size() << " (" << (&fXrefPINHG1) << ")" << endl;
00457   //cout << LILinResp::fXrefPINHG2.size() << " (" << (&fXrefPINHG2) << ")" << endl;
00458   //cout << LILinResp::fXrefPINHG3.size() << " (" << (&fXrefPINHG3) << ")" << endl;
00459   //cout << LILinResp::fXrefPINHG4.size() << " (" << (&fXrefPINHG4) << ")" << endl;
00460   //cout << LILinResp::fXrefPINHG5.size() << " (" << (&fXrefPINHG5) << ")" << endl;
00461   //cout << LILinResp::fXrefPINHG6.size() << " (" << (&fXrefPINHG6) << ")" << endl;
00462 
00463   //cout << LILinResp::fXrefPINLG1.size() << endl;
00464   //cout << LILinResp::fXrefPINLG2.size() << endl;
00465   //cout << LILinResp::fXrefPINLG3.size() << endl;
00466   //cout << LILinResp::fXrefPINLG4.size() << endl;
00467   //cout << LILinResp::fXrefPINLG5.size() << endl;
00468   //cout << LILinResp::fXrefPINLG6.size() << endl;
00469 
00470   //cout << "SlopeHG: " << slopeHG << "\t" << interceptHG << endl;
00471   //cout << "SlopeHG_e: " << slopeHG_e << "\t" << interceptHG_e << endl;
00472 
00473   //cout << "SlopeLG: " << slopeLG << "\t" << interceptLG << endl;
00474   //cout << "SlopeLG_e: " << slopeLG_e << "\t" << interceptLG_e << endl << endl;
00475 
00476 
00477   return;
00478 
00479 }
00480 
00481 //EOF