SFHMarana: TMarJetCal.C Source File

00001 
00002 // Class TMarJetCal
00003 //
00004 // Author     : Emmanuel Sauvan
00005 // Created    : 12/2002
00006 // Last update: $Date: 2005/03/11 08:52:43 $
00007 //          by: $Author: sauvan $
00008 // Comment: Contains definition of MarJetCalibration class
00010 
00011 #ifndef __TMARJETCAL
00012 #include "Marana/TMarJetCal.h"
00013 #endif
00014 
00015 #define MARJETCALDEBUG 0
00016 
00017 
00018 
00019 //ClassImp(TMarJetCal)
00020 
00021 
00022 TMarJetCal::TMarJetCal()
00023 {
00024  for(Int_t i=0;i<NBTHETABINS;i++) {
00025         fStep1Coef[i].val=1.;
00026         fStep2CoefData[i].val[0]=1.;
00027         fStep2CoefData[i].val[1]=0.;
00028         fStep2CoefData[i].val[2]=0.0000001;
00029         fStep2CoefMC[i].val[0]=1.;
00030         fStep2CoefMC[i].val[1]=0.;
00031         fStep2CoefMC[i].val[2]=0.0000001;
00032  }
00033  
00034 }
00035 
00036 
00037 TMarJetCal::TMarJetCal(Char_t *FileStep1,Char_t *FileStep2Data,Char_t *FileStep2MC)
00038 {
00040 //     constructor, coefficients read from text files
00042 
00043 //---- Init all
00044  for(Int_t i=0;i<NBTHETABINS;i++) {
00045         fStep1Coef[i].val=1.;
00046         fStep2CoefData[i].val[0]=1.;
00047         fStep2CoefData[i].val[1]=0.;
00048         fStep2CoefData[i].val[2]=0.0000001;
00049         fStep2CoefMC[i].val[0]=1.;
00050         fStep2CoefMC[i].val[1]=0.;
00051         fStep2CoefMC[i].val[2]=0.0000001;
00052  }
00053 
00054 printf(">>> JetCalibration from %s ... \n",FileStep1);
00055 //---- Read step 1 coefficients
00056  ifstream in(FileStep1,ios::in);        // open file
00057  if(!in) {
00058         cerr<<">> MarJetCal : File "<<FileStep1<<" for step 1 not found"<<endl;
00059         return;
00060  }
00061  Int_t ientries=0;
00062  while(1) {
00063         in >> fStep1Coef[ientries].min >> fStep1Coef[ientries].max >> fStep1Coef[ientries].val;
00064         printf("%f < Theta < %f : %f \n",fStep1Coef[ientries].min,fStep1Coef[ientries].max,fStep1Coef[ientries].val);
00065         if (!in.good()||ientries>=NBTHETABINS) break;
00066         ientries++;
00067  }
00068  in.close();
00069 
00070 //---- Read step 2 Data 
00071 printf(">>> JetCalibration from %s ... \n",FileStep2Data);
00072  ifstream in2d(FileStep2Data,ios::in);  // open file
00073  if(!in2d) {
00074         cerr<<">> MarJetCal : File "<<FileStep2Data<<" for step 2 not found"<<endl;
00075         return;
00076  }
00077  ientries=0;
00078  while(1) {
00079         in2d >> fStep2CoefData[ientries].min >> fStep2CoefData[ientries].max 
00080              >> fStep2CoefData[ientries].val[0] >> fStep2CoefData[ientries].val[1] 
00081              >> fStep2CoefData[ientries].val[2];
00082         printf("%f < Theta < %f : A=%f B=%f C=%f D=%f E=%f\n",fStep2CoefData[ientries].min,fStep2CoefData[ientries].max,
00083         fStep2CoefData[ientries].val[0],fStep2CoefData[ientries].val[1],fStep2CoefData[ientries].val[2]);
00084         if (!in2d.good()||ientries>=NBTHETABINS) break;
00085         ientries++;
00086  }
00087  in2d.close();
00088 
00089 //---- Read step 2 MC 
00090 printf(">>> JetCalibration from %s ... \n",FileStep2MC);
00091  ifstream in2mc(FileStep2MC,ios::in);   // open file
00092  if(!in2mc) {
00093         cerr<<">> MarJetCal : File "<<FileStep2Data<<" for step 2 not found"<<endl;
00094         return;
00095  }
00096  ientries=0;
00097  while(1) {
00098         in2mc >> fStep2CoefMC[ientries].min >> fStep2CoefMC[ientries].max 
00099              >> fStep2CoefMC[ientries].val[0] >> fStep2CoefMC[ientries].val[1] 
00100              >> fStep2CoefMC[ientries].val[2] ;
00101         printf("%f < Theta < %f : A=%f B=%f C=%f D=%f E=%f\n",fStep2CoefMC[ientries].min,fStep2CoefMC[ientries].max,
00102         fStep2CoefMC[ientries].val[0],fStep2CoefMC[ientries].val[1],fStep2CoefMC[ientries].val[2]
00103         ,fStep2CoefMC[ientries].val[3],fStep2CoefMC[ientries].val[4],fStep2CoefMC[ientries].val[5]);
00104         if (!in2mc.good()||ientries>=NBTHETABINS) break;
00105         ientries++;
00106  }
00107  in2mc.close();
00108 
00109 }
00110 
00111 
00112 TMarJetCal::TMarJetCal(Char_t *FileStep2Data,Char_t *FileStep2MC)
00113 {
00115 //     constructor, coefficients read from text files
00117 
00118 //---- Init all
00119  for(Int_t i=0;i<NBTHETABINS;i++) {
00120         fStep2CoefData[i].val[0]=1.;
00121         fStep2CoefData[i].val[1]=1.;
00122         fStep2CoefData[i].val[2]=10.;
00123         fStep2CoefMC[i].val[0]=1.;
00124         fStep2CoefMC[i].val[1]=1.;
00125         fStep2CoefMC[i].val[2]=10.;
00126  }
00127 
00128 printf(">>> JetCalibration from %s ... \n",FileStep2Data);
00129 
00130  ifstream in2d(FileStep2Data,ios::in);  // open file
00131  if(!in2d) {
00132         cerr<<">> MarJetCal : File "<<FileStep2Data<<" for step 2 not found"<<endl;
00133         return;
00134  }
00135 Int_t ientries=0;
00136  while(1) {
00137         in2d >> fStep2CoefData[ientries].min >> fStep2CoefData[ientries].max 
00138              >> fStep2CoefData[ientries].val[0] >> fStep2CoefData[ientries].val[1] 
00139              >> fStep2CoefData[ientries].val[2];
00140         printf("%f < Theta < %f : A=%f B=%f C=%f\n",fStep2CoefData[ientries].min,fStep2CoefData[ientries].max,
00141         fStep2CoefData[ientries].val[0],fStep2CoefData[ientries].val[1],fStep2CoefData[ientries].val[2]);
00142         if (!in2d.good()||ientries>=NBTHETABINS) break;
00143         ientries++;
00144  }
00145  in2d.close();
00146 
00147 //---- Read step 2 MC 
00148 printf(">>> JetCalibration from %s ... \n",FileStep2MC);
00149  ifstream in2mc(FileStep2MC,ios::in);   // open file
00150  if(!in2mc) {
00151         cerr<<">> MarJetCal : File "<<FileStep2Data<<" for step 2 not found"<<endl;
00152         return;
00153  }
00154  ientries=0;
00155  while(1) {
00156         in2mc >> fStep2CoefMC[ientries].min >> fStep2CoefMC[ientries].max 
00157              >> fStep2CoefMC[ientries].val[0] >> fStep2CoefMC[ientries].val[1] 
00158              >> fStep2CoefMC[ientries].val[2] ;
00159         printf("%f < Theta < %f : A=%f B=%f C=%f\n",fStep2CoefMC[ientries].min,fStep2CoefMC[ientries].max,
00160         fStep2CoefMC[ientries].val[0],fStep2CoefMC[ientries].val[1],fStep2CoefMC[ientries].val[2]);
00161         if (!in2mc.good()||ientries>=NBTHETABINS) break;
00162         ientries++;
00163  }
00164  in2mc.close();
00165 
00166 }
00167 
00168 
00169 
00170 TMarJetCal::TMarJetCal(Char_t *FileStep2Data,Char_t *FileStep2MC,Char_t *FileStep2DataRatio, Char_t *FileStep2MCRatio)
00171 {
00173 //     constructor, coefficients read from text files
00175 
00176 //---- Init all
00177  for(Int_t i=0;i<NBTHETABINS;i++) {
00178         fStep2CoefData[i].val[0]=1.;
00179         fStep2CoefData[i].val[1]=1.;
00180         fStep2CoefData[i].val[2]=10.;
00181         fStep2CoefMC[i].val[0]=1.;
00182         fStep2CoefMC[i].val[1]=1.;
00183         fStep2CoefMC[i].val[2]=10.;
00184         fStep2CoefDataRatio[i].val[0]=1.;
00185         fStep2CoefDataRatio[i].val[1]=1.;
00186         fStep2CoefDataRatio[i].val[2]=10.;
00187         fStep2CoefMCRatio[i].val[0]=1.;
00188         fStep2CoefMCRatio[i].val[1]=1.;
00189         fStep2CoefMCRatio[i].val[2]=10.;
00190  }
00191 
00192 printf(">>> JetCalibration from %s ... \n",FileStep2Data);
00193 
00194  ifstream in2d(FileStep2Data,ios::in);  // open file
00195  if(!in2d) {
00196         cerr<<">> MarJetCal : File "<<FileStep2Data<<" for step 2 not found"<<endl;
00197         return;
00198  }
00199 Int_t ientries=0;
00200 // cout << "NBTHETABINS " << NBTHETABINS <<endl;
00201  while(1) {
00202    //   cout << " ientries " << ientries << endl;
00203         in2d >> fStep2CoefData[ientries].min >> fStep2CoefData[ientries].max 
00204              >> fStep2CoefData[ientries].val[0] >> fStep2CoefData[ientries].val[1] 
00205              >> fStep2CoefData[ientries].val[2];
00206         printf("%f < Theta < %f : A=%f B=%f C=%f\n",fStep2CoefData[ientries].min,fStep2CoefData[ientries].max,
00207         fStep2CoefData[ientries].val[0],fStep2CoefData[ientries].val[1],fStep2CoefData[ientries].val[2]);
00208         
00209         ientries++;
00210         if (!in2d.good()||ientries>=NBTHETABINS) break;
00211 
00212  }
00213  in2d.close();
00214 
00215 printf(">>> JetCalibration from %s ... \n",FileStep2DataRatio);
00216 
00217  ifstream in2dratio(FileStep2DataRatio,ios::in);        // open file
00218  if(!in2dratio) {
00219         cerr<<">> MarJetCal : File "<<FileStep2DataRatio<<" for step 2 not found"<<endl;
00220         return;
00221  }
00222 Int_t ientries2=0;
00223  while(1) {
00224         in2dratio >> fStep2CoefDataRatio[ientries2].min >> fStep2CoefDataRatio[ientries2].max 
00225              >> fStep2CoefDataRatio[ientries2].val[0] >> fStep2CoefDataRatio[ientries2].val[1] 
00226              >> fStep2CoefDataRatio[ientries2].val[2];
00227         printf("%f < Theta < %f : A=%f B=%f C=%f\n",fStep2CoefDataRatio[ientries2].min,fStep2CoefDataRatio[ientries2].max,
00228         fStep2CoefDataRatio[ientries2].val[0],fStep2CoefDataRatio[ientries2].val[1],fStep2CoefDataRatio[ientries2].val[2]);
00229         ientries2++;
00230         if (!in2dratio.good()||ientries2>=NBTHETABINS) break;
00231 
00232  }
00233  in2dratio.close();
00234 
00235 
00236 
00237 //---- Read step 2 MC 
00238 printf(">>> JetCalibration from %s ... \n",FileStep2MC);
00239  ifstream in2mc(FileStep2MC,ios::in);   // open file
00240  if(!in2mc) {
00241         cerr<<">> MarJetCal : File "<<FileStep2Data<<" for step 2 not found"<<endl;
00242         return;
00243  }
00244  ientries=0;
00245  while(1) {
00246         in2mc >> fStep2CoefMC[ientries].min >> fStep2CoefMC[ientries].max 
00247              >> fStep2CoefMC[ientries].val[0] >> fStep2CoefMC[ientries].val[1] 
00248              >> fStep2CoefMC[ientries].val[2] ;
00249         printf("%f < Theta < %f : A=%f B=%f C=%f\n",fStep2CoefMC[ientries].min,fStep2CoefMC[ientries].max,
00250         fStep2CoefMC[ientries].val[0],fStep2CoefMC[ientries].val[1],fStep2CoefMC[ientries].val[2]);
00251         ientries++;
00252         if (!in2mc.good()||ientries>=NBTHETABINS) break;
00253         
00254  }
00255  in2mc.close();
00256 
00257 printf(">>> JetCalibration from %s ... \n",FileStep2MCRatio);
00258  ifstream in2mcratio(FileStep2MCRatio,ios::in); // open file
00259  if(!in2mcratio) {
00260         cerr<<">> MarJetCal : File "<<FileStep2Data<<" for step 2 not found"<<endl;
00261         return;
00262  }
00263  ientries=0;
00264  while(1) {
00265         in2mcratio >> fStep2CoefMCRatio[ientries].min >> fStep2CoefMCRatio[ientries].max 
00266              >> fStep2CoefMCRatio[ientries].val[0] >> fStep2CoefMCRatio[ientries].val[1] 
00267              >> fStep2CoefMCRatio[ientries].val[2] ;
00268         printf("%f < Theta < %f : A=%f B=%f C=%f\n",fStep2CoefMCRatio[ientries].min,fStep2CoefMCRatio[ientries].max,
00269         fStep2CoefMCRatio[ientries].val[0],fStep2CoefMCRatio[ientries].val[1],fStep2CoefMCRatio[ientries].val[2]);
00270         ientries++;
00271         if (!in2mcratio.good()||ientries>=NBTHETABINS) break;
00272      
00273  }
00274  in2mcratio.close();
00275 
00276 }
00277 
00278 
00279 // Standard destructor
00280 TMarJetCal::~TMarJetCal()
00281 {
00282 
00283 };
00284 
00285 
00286 
00287 
00288 Float_t TMarJetCal::ApplyCalib(Float_t Ptjetnocal,Float_t Thetajet,Int_t MCFlag,Int_t option=2)
00289 {
00291 //     Apply jet calibrations  (option = 1(step one only), 2(2 steps), 0 nocal
00292 //              Theta in radians !!!!!
00294 
00295  Float_t Ptjet=Ptjetnocal;
00296 
00297  if(option==0) return Ptjet;    //---- no calibration applied
00298  //---- Apply step 1 calibration (only for data)
00299 //  if(MCFlag==0&&option!=2) {
00300 //      for(Int_t i=0;i<NBTHETABINS;i++) {
00301 //              if(Thetajet*MYR2D>fStep1Coef[i].min&&Thetajet*MYR2D<=fStep1Coef[i].max) { 
00302 //                      if(fStep1Coef[i].val) Ptjet/=fStep1Coef[i].val;
00303 //                      break;
00304 //              }
00305 //      }
00306 //  }
00307 // printf("Correction step 1 : %f to %f \n",Ptjetnocal,Ptjet);
00308  //--- Apply step 2 calibration
00309 
00310   if(option!=2) return Ptjet;   //---- no step2 calibration applied
00311   if(Ptjet<5||Thetajet*MYR2D>155) return Ptjet; //---- no calibration if Ptjet < 10 GeV
00312   Float_t A=0,B=0,C=0;
00313   Float_t factor=0;
00314 
00315   for(Int_t j=0;j<1;j++) {      //---- make 2 iterations for correction
00316     if(MCFlag==0) {
00317       for(Int_t i=0;i<NBTHETABINS;i++) {
00318             if(Thetajet*MYR2D>fStep2CoefData[i].min&&Thetajet*MYR2D<=fStep2CoefData[i].max) {
00319                     A=fStep2CoefData[i].val[0];
00320                     B=fStep2CoefData[i].val[1];
00321                     C=fStep2CoefData[i].val[2];
00322         
00323                     factor=1;
00324                     break;
00325             }
00326       }
00327      }else{
00328       for(Int_t i=0;i<NBTHETABINS;i++) {
00329             if(Thetajet*MYR2D>fStep2CoefMC[i].min&&Thetajet*MYR2D<=fStep2CoefMC[i].max) {
00330                     A=fStep2CoefMC[i].val[0];
00331                     B=fStep2CoefMC[i].val[1];
00332                     C=fStep2CoefMC[i].val[2]; 
00333                     factor=1;
00334                     break;
00335             }
00336       }
00337     }
00338     if(factor==1) {
00339         Float_t Ptjetcor=Ptjet;
00340 //      if(C) factor=A*(1-exp(-B-Ptjet/C))+D*exp(-E*Ptjet);
00341 //      if(factor) Ptjetcor=Ptjet/factor;
00342         if (A !=0. && C !=0.) factor=A*(1-exp(-B-Ptjetcor/C));
00343         if(factor) Ptjet/=factor;
00344         
00345     }
00346     //        printf("Calib factor : %f \n",factor);
00347 
00348   }
00349 
00350 return Ptjet;
00351 
00352 }
00353 
00354 
00355 
00356 
00357 
00358 
00359 // TMarJetCal::TMarJetCal()
00360 // {
00361 //  for(Int_t i=0;i<NBTHETABINS;i++) {
00362 //      fStep1Coef[i].val=1.;
00363 //      fStep2CoefData[i].val[0]=1.;
00364 //      fStep2CoefData[i].val[1]=0.;
00365 //      fStep2CoefData[i].val[2]=0.0000001;
00366 //      fStep2CoefData[i].val[3]=1.;
00367 //      fStep2CoefData[i].val[4]=0.0;
00368 //      fStep2CoefData[i].val[5]=0.000001;
00369 //      fStep2CoefMC[i].val[0]=1.;
00370 //      fStep2CoefMC[i].val[1]=0.;
00371 //      fStep2CoefMC[i].val[2]=0.0000001;
00372 //      fStep2CoefMC[i].val[3]=1.0;
00373 //      fStep2CoefMC[i].val[4]=0.0;
00374 //      fStep2CoefMC[i].val[5]=0.0000001;
00375 //  }
00376  
00377 // }
00378 
00379 
00380 // TMarJetCal::TMarJetCal(Char_t *FileStep1,Char_t *FileStep2Data,Char_t *FileStep2MC)
00381 // {
00382 // ///////////////////////////////////////////////////////////// 
00383 // //     constructor, coefficients read from text files
00384 // ///////////////////////////////////////////////////////////// 
00385 
00386 // //---- Init all
00387 //  for(Int_t i=0;i<NBTHETABINS;i++) {
00388 //      fStep1Coef[i].val=1.;
00389 //      fStep2CoefData[i].val[0]=1.;
00390 //      fStep2CoefData[i].val[1]=0.;
00391 //      fStep2CoefData[i].val[2]=0.0000001;
00392 //      fStep2CoefData[i].val[3]=1.0;
00393 //      fStep2CoefData[i].val[4]=0.0;
00394 //      fStep2CoefData[i].val[5]=0.00001;
00395 //      fStep2CoefMC[i].val[0]=1.;
00396 //      fStep2CoefMC[i].val[1]=0.;
00397 //      fStep2CoefMC[i].val[2]=0.0000001;
00398 //      fStep2CoefMC[i].val[3]=1.;
00399 //      fStep2CoefMC[i].val[4]=0.0;
00400 //      fStep2CoefMC[i].val[5]=0.000001;
00401 //  }
00402 
00403 // printf(">>> JetCalibration from %s ... \n",FileStep1);
00404 // //---- Read step 1 coefficients
00405 //  ifstream in(FileStep1,ios::in);     // open file
00406 //  if(!in) {
00407 //      cerr<<">> MarJetCal : File "<<FileStep1<<" for step 1 not found"<<endl;
00408 //      return;
00409 //  }
00410 //  Int_t ientries=0;
00411 //  while(1) {
00412 //      in >> fStep1Coef[ientries].min >> fStep1Coef[ientries].max >> fStep1Coef[ientries].val;
00413 //      printf("%f < Theta < %f : %f \n",fStep1Coef[ientries].min,fStep1Coef[ientries].max,fStep1Coef[ientries].val);
00414 //      if (!in.good()||ientries>=NBTHETABINS) break;
00415 //      ientries++;
00416 //  }
00417 //  in.close();
00418 
00419 // //---- Read step 2 Data 
00420 // printf(">>> JetCalibration from %s ... \n",FileStep2Data);
00421 //  ifstream in2d(FileStep2Data,ios::in);       // open file
00422 //  if(!in2d) {
00423 //      cerr<<">> MarJetCal : File "<<FileStep2Data<<" for step 2 not found"<<endl;
00424 //      return;
00425 //  }
00426 //  ientries=0;
00427 //  while(1) {
00428 //      in2d >> fStep2CoefData[ientries].min >> fStep2CoefData[ientries].max 
00429 //           >> fStep2CoefData[ientries].val[0] >> fStep2CoefData[ientries].val[1] 
00430 //           >> fStep2CoefData[ientries].val[2] >> fStep2CoefData[ientries].val[3]
00431 //           >> fStep2CoefData[ientries].val[4]  >> fStep2CoefData[ientries].val[5];
00432 //      printf("%f < Theta < %f : A=%f B=%f C=%f D=%f E=%f\n",fStep2CoefData[ientries].min,fStep2CoefData[ientries].max,
00433 //      fStep2CoefData[ientries].val[0],fStep2CoefData[ientries].val[1],fStep2CoefData[ientries].val[2]
00434 //      ,fStep2CoefData[ientries].val[3],fStep2CoefData[ientries].val[4],fStep2CoefData[ientries].val[5]);
00435 //      if (!in2d.good()||ientries>=NBTHETABINS) break;
00436 //      ientries++;
00437 //  }
00438 //  in2d.close();
00439 
00440 // //---- Read step 2 MC 
00441 // printf(">>> JetCalibration from %s ... \n",FileStep2MC);
00442 //  ifstream in2mc(FileStep2MC,ios::in);        // open file
00443 //  if(!in2mc) {
00444 //      cerr<<">> MarJetCal : File "<<FileStep2Data<<" for step 2 not found"<<endl;
00445 //      return;
00446 //  }
00447 //  ientries=0;
00448 //  while(1) {
00449 //      in2mc >> fStep2CoefMC[ientries].min >> fStep2CoefMC[ientries].max 
00450 //           >> fStep2CoefMC[ientries].val[0] >> fStep2CoefMC[ientries].val[1] 
00451 //           >> fStep2CoefMC[ientries].val[2] >> fStep2CoefMC[ientries].val[3] 
00452 //           >> fStep2CoefMC[ientries].val[4] >> fStep2CoefMC[ientries].val[5];
00453 //      printf("%f < Theta < %f : A=%f B=%f C=%f D=%f E=%f\n",fStep2CoefMC[ientries].min,fStep2CoefMC[ientries].max,
00454 //      fStep2CoefMC[ientries].val[0],fStep2CoefMC[ientries].val[1],fStep2CoefMC[ientries].val[2]
00455 //      ,fStep2CoefMC[ientries].val[3],fStep2CoefMC[ientries].val[4],fStep2CoefMC[ientries].val[5]);
00456 //      if (!in2mc.good()||ientries>=NBTHETABINS) break;
00457 //      ientries++;
00458 //  }
00459 //  in2mc.close();
00460 
00461 // }
00462 
00463 
00464 
00465 // // Standard destructor
00466 // TMarJetCal::~TMarJetCal()
00467 // {
00468 
00469 // };
00470 
00471 
00472 
00473 
00474 // Float_t TMarJetCal::ApplyCalib(Float_t Ptjetnocal,Float_t Thetajet,Int_t MCFlag,Int_t option=2)
00475 // {
00476 // ///////////////////////////////////////////////////////////// 
00477 // //     Apply jet calibrations  (option = 1(step one only), 2(2 steps), 0 nocal
00478 // //           Theta in radians !!!!!
00479 // ///////////////////////////////////////////////////////////// 
00480 
00481 //  Float_t Ptjet=Ptjetnocal;
00482 
00483 //  if(option==0) return Ptjet; //---- no calibration applied
00484 //  //---- Apply step 1 calibration (only for data)
00485 // //  if(MCFlag==0&&option!=2) {
00486 // //   for(Int_t i=0;i<NBTHETABINS;i++) {
00487 // //           if(Thetajet*MYR2D>fStep1Coef[i].min&&Thetajet*MYR2D<=fStep1Coef[i].max) { 
00488 // //                   if(fStep1Coef[i].val) Ptjet/=fStep1Coef[i].val;
00489 // //                   break;
00490 // //           }
00491 // //   }
00492 // //  }
00493 // // printf("Correction step 1 : %f to %f \n",Ptjetnocal,Ptjet);
00494 //  //--- Apply step 2 calibration
00495 
00496 //   if(option!=2) return Ptjet;        //---- no step2 calibration applied
00497 //   if(Ptjet<5||Thetajet*MYR2D>155) return Ptjet;      //---- no calibration if Ptjet < 10 GeV
00498 //   Float_t A=0,B=0,C=0,D=0,E=0,F=0;
00499 //   Float_t factor=0;
00500 
00501 //   for(Int_t j=0;j<1;j++) {   //---- make 2 iterations for correction
00502 //     if(MCFlag==0) {
00503 //       for(Int_t i=0;i<NBTHETABINS;i++) {
00504 //          if(Thetajet*MYR2D>fStep2CoefData[i].min&&Thetajet*MYR2D<=fStep2CoefData[i].max) {
00505 //                  A=fStep2CoefData[i].val[0];
00506 //                  B=fStep2CoefData[i].val[1];
00507 //                  C=fStep2CoefData[i].val[2];
00508 //                  D=fStep2CoefData[i].val[3];
00509 //                  E=fStep2CoefData[i].val[4];
00510 //                  F=fStep2CoefData[i].val[5];
00511 //                  factor=1;
00512 //                  break;
00513 //          }
00514 //       }
00515 //      }else{
00516 //       for(Int_t i=0;i<NBTHETABINS;i++) {
00517 //          if(Thetajet*MYR2D>fStep2CoefMC[i].min&&Thetajet*MYR2D<=fStep2CoefMC[i].max) {
00518 //                  A=fStep2CoefMC[i].val[0];
00519 //                  B=fStep2CoefMC[i].val[1];
00520 //                  C=fStep2CoefMC[i].val[2];
00521 //                  D=fStep2CoefMC[i].val[3];
00522 //                  E=fStep2CoefMC[i].val[4];
00523 //                  F=fStep2CoefMC[i].val[5];
00524 //                  factor=1;
00525 //                  break;
00526 //          }
00527 //       }
00528 //     }
00529 //     if(factor==1) {
00530 //         Float_t Ptjetcor=Ptjet;
00531 // //           if(C) factor=A*(1-exp(-B-Ptjet/C))+D*exp(-E*Ptjet);
00532 // //           if(factor) Ptjetcor=Ptjet/factor;
00533 //      if (A !=0 && C !=0 && D !=0 && F !=0)
00534 //       factor= (1+A*(1-exp(-B-Ptjetcor/C))*(1-D*(1-exp(-E-Ptjetcor/F))))/(A*(1-exp(-B-Ptjetcor/C))*D*(1-exp(-E-Ptjetcor/F));
00535 //      if(factor) Ptjet*=factor;
00536         
00537 //     }
00538 //     printf("Calib factor : %f \n",factor);
00539 
00540 //   }
00541 
00542 // return Ptjet;
00543 
00544 // }
00545 
00546 
00547 TLorentzVector TMarJetCal::ApplyHadroo2Calib(TLorentzVector tracks,TLorentzVector clusters,Float_t Thetajet,Short_t MCFlag)
00548 {
00550 //     Apply jet calibrations  (option = 1(step one only), 2(2 steps), 0 nocal
00551 //              Theta in radians !!!!!
00553   TLorentzVector hadrons=tracks+clusters;
00554   TLorentzVector calib(0.,0.,0.,0.);
00555 
00556  Float_t Ptjet=hadrons.Pt();
00557  // Float_t Thetajet=hadrons.Theta();
00558 
00559  // cout << " start hadroo2 calib  " << endl;
00560  //  cout << " theta " << Thetajet*MYR2D <<  endl;
00561  if(Ptjet<5||Thetajet*MYR2D>155||Thetajet*MYR2D<7) {
00562    //   cout << " low pt " << Ptjet << " theta " << Thetajet*MYR2D << endl;
00563   return hadrons;       //---- no calibration if Ptjet < 5
00564 
00565 }
00566   Float_t A=0,B=0,C=0,D=0,E=0,F=0;
00567   Float_t factor=0;
00568 
00569   // let us try without
00570   // for(Int_t j=0;j<1;j++) {   //---- make 2 iterations for correction
00571 
00572   // Ptjet=hadrons.Pt();
00573   // Thetajet=hadrons.Theta();
00574   // cout << "  pt " << Ptjet << " theta " << Thetajet*MYR2D << "  MCFlag " << MCFlag <<  endl;
00575     if(MCFlag==0) {
00576       for(Int_t i=0;i<NBTHETABINS;i++) {
00577             if(Thetajet*MYR2D>fStep2CoefData[i].min&&Thetajet*MYR2D<=fStep2CoefData[i].max) {
00578                     A=fStep2CoefData[i].val[0];
00579                     //    cout << "ici A= " << A << " bin i " << i << endl;
00580                     B=fStep2CoefData[i].val[1];
00581                     C=fStep2CoefData[i].val[2];
00582                     D=fStep2CoefDataRatio[i].val[0];
00583                     E=fStep2CoefDataRatio[i].val[1];
00584                     F=fStep2CoefDataRatio[i].val[2];
00585                     factor=1;
00586                     break;
00587             }
00588       }
00589      }else{
00590       for(Int_t i=0;i<NBTHETABINS;i++) {
00591             if(Thetajet*MYR2D>fStep2CoefMC[i].min&&Thetajet*MYR2D<=fStep2CoefMC[i].max) {
00592                     A=fStep2CoefMC[i].val[0];
00593                     //                  cout << "et non ici A= " << A << endl;
00594                     B=fStep2CoefMC[i].val[1];
00595                     C=fStep2CoefMC[i].val[2];
00596                     D=fStep2CoefMCRatio[i].val[0];
00597                     E=fStep2CoefMCRatio[i].val[1];
00598                     F=fStep2CoefMCRatio[i].val[2];
00599                     
00600                     factor=1;
00601                     break;
00602             }
00603       }
00604     }
00605     //  if(factor==1) {
00606    
00607     //   cout << "here " <<" A " <<A <<" B " <<B <<" C " <<C  <<" D" <<D <<" E " <<E <<" F " <<F <<endl;
00608     Float_t Ptjetcor=Ptjet;
00609 //      if(C) factor=A*(1-exp(-B-Ptjet/C))+D*exp(-E*Ptjet);
00610 //      if(factor) Ptjetcor=Ptjet/factor;
00611         Float_t Fptbal,Fc;
00612 
00613         Fptbal=A*(1-exp(-B-Ptjetcor/C));
00614         if (Thetajet*MYR2D<30) {
00615                 Fc=D*(1-exp(-E-Ptjetcor/F));
00616         }else{
00617                 Fc=D*(1+exp(-E+Ptjetcor/F));
00618         }
00619 
00620         //factor=(1+Fptbal*(Fc-1))/(Fptbal*Fc);
00621         //Ptjetcor=(tracks+factor*clusters).Pt();
00622         
00623         factor=1/Fptbal;
00624         Ptjetcor=(factor*(tracks+clusters)).Pt();
00625         
00626         Fptbal=A*(1-exp(-B-Ptjetcor/C));
00627         if (Thetajet*MYR2D<30) {
00628                 Fc=D*(1-exp(-E-Ptjetcor/F));
00629         }else{
00630                 Fc=D*(1+exp(-E+Ptjetcor/F));
00631         }
00632         //factor=(1+Fptbal*(Fc-1))/(Fptbal*Fc);
00633 
00634         
00635         //      cout << factor << endl;
00636         //      clusters*=factor;
00637 
00638         factor=1/Fptbal;
00639         //hadrons=tracks+factor*clusters;
00640         hadrons=(tracks+clusters)*factor;
00641         
00642         
00643         //      if (A !=0. && C !=0.) factor=A*(1-exp(-B-Ptjetcor/C));
00644         //if(factor) Ptjet/=factor;
00645         
00646         //  }
00647         //      printf("here Calib factor : %f \n",factor);
00648         calib=hadrons;
00649 
00650     
00651     //  }
00652   return calib;
00653   //return Ptjet;
00654 
00655 }
00656 
00657 
00658 
00659 // --------------------------------------------------------------------------------------------------------------------------
00660 // ------------------------------------ My Hadroo2 Kt Jet Calibration ----------------------------------
00661 //---------------------------------------------------------------------------------------------------------------------------
00662 
00663 TLorentzVector TMarJetCal::MyHadroo2KtJetCalibration(H1PartJet* UncalJet,Short_t MCFlag)
00664 {
00665   // Form the uncalibrated Jet, access fourvectors of tracks and clusters, 
00666   // array of fourvectors ...
00667   // calibrates the tracks component, 
00668 
00669   // do not calibrate jets outside LAr or Low Pt jets.
00670   // say, below 5 GeV
00671  if(UncalJet->GetPt()<4.
00672     ||UncalJet->GetTheta()*MYR2D>155
00673     ||UncalJet->GetTheta()*MYR2D<7) {
00674    return UncalJet->GetFourVector();    
00675 }
00676 
00677  Double_t Thetajet=UncalJet->GetTheta();
00678 
00679   Float_t A=0,B=0,C=0,D=0,E=0,F=0;
00680   Float_t factor=0;
00681     if(MCFlag==0) {
00682       for(Int_t i=0;i<NBTHETABINS;i++) {
00683             if(Thetajet*MYR2D>fStep2CoefData[i].min&&Thetajet*MYR2D<=fStep2CoefData[i].max) {
00684                     A=fStep2CoefData[i].val[0];
00685                     //    cout << "ici A= " << A << " bin i " << i << endl;
00686                     B=fStep2CoefData[i].val[1];
00687                     C=fStep2CoefData[i].val[2];
00688                     D=fStep2CoefDataRatio[i].val[0];
00689                     E=fStep2CoefDataRatio[i].val[1];
00690                     F=fStep2CoefDataRatio[i].val[2];
00691                     factor=1;
00692                     break;
00693             }
00694       }
00695      }else{
00696       for(Int_t i=0;i<NBTHETABINS;i++) {
00697             if(Thetajet*MYR2D>fStep2CoefMC[i].min&&Thetajet*MYR2D<=fStep2CoefMC[i].max) {
00698                     A=fStep2CoefMC[i].val[0];
00699                     //                  cout << "et non ici A= " << A << endl;
00700                     B=fStep2CoefMC[i].val[1];
00701                     C=fStep2CoefMC[i].val[2];
00702                     D=fStep2CoefMCRatio[i].val[0];
00703                     E=fStep2CoefMCRatio[i].val[1];
00704                     F=fStep2CoefMCRatio[i].val[2];
00705                     
00706                     factor=1;
00707                     break;
00708             }
00709       }
00710     }
00711 
00712     
00713 
00714     Float_t Ptjetcor=UncalJet->GetPt();
00715     Float_t Fptbal;
00716     Float_t Fc;
00717    
00718     Fptbal=A*(1-exp(-B-Ptjetcor/C));
00719    const TObjArray *JetParticles = UncalJet->GetParticles();
00720     Int_t NumParticles = JetParticles->GetEntriesFast();
00721     TLorentzVector tracks(0,0,0,0),clusters(0,0,0,0);
00722 
00723   
00724     // moi : et si on calculais Fc ??????????????? 
00725 //    if (Thetajet*MYR2D<30) {
00726 //       Fc=D*(1-exp(-E-Ptjetcor/F));
00727 //     }else{
00728 //       Fc=D*(1+exp(-E+Ptjetcor/F));
00729 //     }
00730     
00731     
00732 
00733       
00734 
00735        // First loop over jets daughters to get the Pt of tracks and clusters
00736 
00737  
00738     // loop over jet particles
00739     for (Int_t ipart=0 ; ipart<NumParticles ; ipart++){
00740          const H1PartCand *jetpartcand = UncalJet->GetPartCand(ipart);
00741          if(!jetpartcand->IsSelTrack()){
00742            clusters+=jetpartcand->GetFourVector();
00743          } else {
00744            tracks+=jetpartcand->GetFourVector();
00745          }
00746          //      jetpart4vect[ipart]->SetPtEtaPhiE(pt,eta,phi,E);
00747     } // end loop over jet particles
00748     Fc=clusters.Pt()/(tracks.Pt()+clusters.Pt());
00749     factor=(1+Fptbal*(Fc-1))/(Fptbal*Fc);
00750     Ptjetcor=tracks.Pt()+factor*(clusters.Pt());
00751 
00752       //        Ptjetcor=(factor*(tracks+clusters)).Pt();       
00753 
00754         Fptbal=A*(1-exp(-B-Ptjetcor/C));
00755         //      if (Thetajet*MYR2D<30) {
00756         //      Fc=D*(1-exp(-E-Ptjetcor/F));
00757         //      }else{
00758         //      Fc=D*(1+exp(-E+Ptjetcor/F));
00759         //      }
00760         factor=(1+Fptbal*(Fc-1))/(Fptbal*Fc);
00761         //      factor=1;
00762         //  const TObjArray *JetParticles = UncalJet->GetParticles();
00763     //  Int_t NumParticles = JetParticles->GetEntriesFast();
00764 
00765 
00766 
00767 
00768 
00769 
00770     TLorentzVector *jetpart4vect= new TLorentzVector[NumParticles];
00771 
00772 
00773     // loop over jet particles
00774     for (Int_t ipart=0 ; ipart<NumParticles ; ipart++){
00775          const H1PartCand *jetpartcand = UncalJet->GetPartCand(ipart);
00776          Double_t pt=jetpartcand->GetPt();
00777          Double_t E=jetpartcand->GetE();
00778          Double_t eta=jetpartcand->GetEta();
00779          Double_t phi=jetpartcand->GetPhi();
00780          
00781          if(!jetpartcand->IsSelTrack()){
00782            pt*=factor;
00783            E*=factor;
00784          }
00785          jetpart4vect[ipart].SetPtEtaPhiE(pt,eta,phi,E);
00786 
00787     } // end loop over jet particles
00788   
00789     // now compute jet total 4 vector.
00790 
00791     TLorentzVector CalJet(0,0,0,0);
00792     //    Double_t Ecal=0;
00793     Double_t Ptcal=0;
00794     Double_t Etacal=0;
00795     Double_t Phical=0;
00796 
00797     //    Double_t E1=0;
00798     Double_t Pt1=0;
00799     Double_t Eta1=0;
00800     Double_t Phi1=0;
00801 
00802     //  Double_t E2=0;
00803     Double_t Pt2=0;
00804     Double_t Eta2=0;
00805     Double_t Phi2=0;
00806 
00807 
00808 
00809     for (Int_t i=0 ; i<NumParticles ; i++ ) {
00810       Pt1=Ptcal;
00811       Eta1=Etacal;
00812       Phi1=Phical;
00813 
00814       Pt2 =jetpart4vect[i].Pt();
00815       Eta2=jetpart4vect[i].Eta();
00816       Phi2=jetpart4vect[i].Phi();
00817       
00818       Ptcal=Pt1+Pt2;
00819       Etacal=(Eta1*Pt1+Eta2*Pt2)/(Pt1+Pt2);
00820       Phical=(Phi1*Pt1+Phi2*Pt2)/(Pt1+Pt2);
00821     }
00822     
00823     CalJet.SetPtEtaPhiM(Ptcal,Etacal,Phical,0.);
00824  
00825   //   cout << "UncalJet->GetE()  " << UncalJet->GetE() << endl;
00826 //     cout << "CalJet.E()  " << CalJet.E() << endl;
00827 //     cout << "UncalJet->GetPt()  " << UncalJet->GetPt() << endl;
00828 //     cout << "CalJet.Pt()  " << CalJet.Pt() << endl;
00829 //     cout << "UncalJet->GetTheta()  " << UncalJet->GetTheta() << endl;
00830 //     cout << "CalJet.Theta()  " << CalJet.Theta() << endl;
00831 //     cout << "UncalJet->GetPhi()  " << UncalJet->GetPhi() << endl;
00832 //     cout << "CalJet.Phi()  " << CalJet.Phi() << endl;
00833     return CalJet;
00834 
00835 }
00836 
00837 
00838 
00839 TLorentzVector TMarJetCal::MyHadroo2TestJetCalibration(H1PartJet* UncalJet,Short_t MCFlag)
00840 {
00841   // Form the uncalibrated Jet, access fourvectors of tracks and clusters, 
00842   // array of fourvectors ...
00843   // calibrates the tracks component, 
00844 
00845   // do not calibrate jets outside LAr or Low Pt jets.
00846  if(UncalJet->GetPt()<5
00847     ||UncalJet->GetTheta()*MYR2D>157
00848     ||UncalJet->GetTheta()*MYR2D<6) {
00849    return UncalJet->GetFourVector();    
00850 }
00851 
00852  Double_t Thetajet=UncalJet->GetTheta();
00853 
00854   Float_t A=0,B=0,C=0,D=0,E=0,F=0;
00855   Float_t factor=0;
00856     if(MCFlag==0) {
00857       for(Int_t i=0;i<NBTHETABINS;i++) {
00858             if(Thetajet*MYR2D>fStep2CoefData[i].min&&Thetajet*MYR2D<=fStep2CoefData[i].max) {
00859                     A=fStep2CoefData[i].val[0];
00860                     //    cout << "ici A= " << A << " bin i " << i << endl;
00861                     B=fStep2CoefData[i].val[1];
00862                     C=fStep2CoefData[i].val[2];
00863                     D=fStep2CoefDataRatio[i].val[0];
00864                     E=fStep2CoefDataRatio[i].val[1];
00865                     F=fStep2CoefDataRatio[i].val[2];
00866                     factor=1;
00867                     break;
00868             }
00869       }
00870      }else{
00871       for(Int_t i=0;i<NBTHETABINS;i++) {
00872             if(Thetajet*MYR2D>fStep2CoefMC[i].min&&Thetajet*MYR2D<=fStep2CoefMC[i].max) {
00873                     A=fStep2CoefMC[i].val[0];
00874                     //                  cout << "et non ici A= " << A << endl;
00875                     B=fStep2CoefMC[i].val[1];
00876                     C=fStep2CoefMC[i].val[2];
00877                     D=fStep2CoefMCRatio[i].val[0];
00878                     E=fStep2CoefMCRatio[i].val[1];
00879                     F=fStep2CoefMCRatio[i].val[2];
00880                     
00881                     factor=1;
00882                     break;
00883             }
00884       }
00885     }
00886 
00887     
00888 
00889     Float_t Ptjetcor=UncalJet->GetPt();
00890     Float_t Fptbal;
00891 
00892     Fptbal=A*(1-exp(-B-Ptjetcor/C));
00893     
00894     Ptjetcor/=Fptbal;
00895     
00896     Fptbal=A*(1-exp(-B-Ptjetcor/C));
00897   
00898     // now compute jet total 4 vector.
00899 
00900     TLorentzVector CalJet(0,0,0,0);
00901     Double_t Ecal=0;
00902     Double_t Pxcal=0;
00903     Double_t Pycal=0;
00904     Double_t Pzcal=0;
00905 
00906     Ecal=UncalJet->GetE()/Fptbal;
00907     Pxcal=UncalJet->GetPx()/Fptbal;
00908     Pycal=UncalJet->GetPy()/Fptbal;
00909     Pzcal=UncalJet->GetPz()/Fptbal;
00910     
00911     CalJet.SetPxPyPzE(Pxcal,Pycal,Pzcal,Ecal);
00912     
00913     return CalJet;
00914 
00915 }