ANALYSIS/AliPhysicsSelection.cxx

   1 /* $Id: AliPhysicsSelection.cxx 35782 2009-10-22 11:54:31Z jgrosseo $ */
   2
   3 /**************************************************************************
   4  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   5  *                                                                        *
   6  * Author: The ALICE Off-line Project.                                    *
   7  * Contributors are mentioned in the code where appropriate.              *
   8  *                                                                        *
   9  * Permission to use, copy, modify and distribute this software and its   *
  10  * documentation strictly for non-commercial purposes is hereby granted   *
  11  * without fee, provided that the above copyright notice appears in all   *
  12  * copies and that both the copyright notice and this permission notice   *
  13  * appear in the supporting documentation. The authors make no claims     *
  14  * about the suitability of this software for any purpose. It is          *
  15  * provided "as is" without express or implied warranty.                  *
  16  **************************************************************************/
  17
  18 //-------------------------------------------------------------------------
  19 //                      Implementation of   Class AliPhysicsSelection
  20 // This class selects collision candidates from data runs, applying selection cuts on triggers
  21 // and background rejection based on the content of the ESD
  22 //
  23 // Usage:
  24 //
  25 // Create the object:
  26 //   fPhysicsSelection = new AliPhysicsSelection;
  27 //
  28 // For MC data, call
  29 //   fPhysicsSelection->SetAnalyzeMC()
  30 //
  31 // To check if an event is a collision candidate, use:
  32 //   fPhysicsSelection->IsCollisionCandidate(fESD)
  33 //
  34 // After processing save the resulting histograms to a file with (a folder physics_selection
  35 //   will be created that contains the histograms):
  36 //   fPhysicsSelection->SaveHistograms("physics_selection")
  37 //
  38 // To print statistics after processing use:
  39 //   fPhysicsSelection->Print();
  40 //
  41 // To seleect the BX ids corresponding to real bunches crossings p2 use:
  42 //   fPhysicsSelection->SetUseBXNumbers();
  43 // you cannot process runs with different filling schemes if you require this option.
  44 //
  45 // To compute the Background automatically using the control triggers
  46 // use:
  47 //   fPhysicsSelection->SetComputeBG();
  48 // this will show the value of the Beam Gas, accidentals and good
  49 // events as additional rows in the statistic tables, but it will NOT
  50 // subtract the background automatically.
  51 // This option will only work for runs taken with the CINT1
  52 // suite. This options enables automatically also the usage of BX
  53 // numbers. You can only process one run at a time if you require this
  54 // options, because it uses the bunch intensity estimated run by run.
  55 //
  56 // The BG will usually be more important in the so-called "bin 0": the
  57 // class can also compute the statistics table for events in this
  58 // bin. Since the definition of bin 0 may in general change from
  59 // analysis to analysis, the user needs to provide a callback
  60 // implementing the definition of bin zero. The callback should be
  61 // implemented as a method in the analysis task and should override
  62 // the IsEventInBinZero method of AliAnalysisTaskSE, and should thus
  63 // have the the following prototype:
  64 //   Bool_t IsEventInBinZero();
  65 // It should return true if the event is in the bin 0 and it is set by
  66 // passing to the physics selection the NAME of the task where the
  67 // callback is implemented:
  68 //   fPhysicsSelection->SetBin0Callback("MyTask").
  69 //
  70 //
  71 // Usually the class selects the trigger scheme by itself depending on the run number.
  72 // Nevertheless, you can do that manually by calling AddCollisionTriggerClass() and AddBGTriggerClass()
  73 // Example:
  74 // To define the class CINT1B-ABCE-NOPF-ALL as collision trigger (those will be accepted as
  75 // collision candidates when they pass the selection):
  76 //   AddCollisionTriggerClass("+CINT1B-ABCE-NOPF-ALL #769 #3119");
  77 // To select on bunch crossing IDs in addition, use:
  78 //   AddCollisionTriggerClass("+CINT1B-ABCE-NOPF-ALL #769 #3119");
  79 // To define the class CINT1A-ABCE-NOPF-ALL as a background trigger (those will only be counted
  80 // for the control histograms):
  81 //   AddBGTriggerClass("+CINT1A-ABCE-NOPF-ALL");
  82 // You can also specify more than one trigger class in a string or you can require that some are *not*
  83 // present. The following line would require CSMBA-ABCE-NOPF-ALL, but CSMBB-ABCE-NOPF-ALL is not allowed
  84 // to be present:
  85 //   AddBGTriggerClass("+CSMBA-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL");
  86 //
  87 //   Origin: Jan Fiete Grosse-Oetringhaus, CERN
  88 //           Michele Floris, CERN
  89 //-------------------------------------------------------------------------
  90
  91 #include <Riostream.h>
  92 #include <TH1F.h>
  93 #include <TH2F.h>
  94 #include <TList.h>
  95 #include <TIterator.h>
  96 #include <TDirectory.h>
  97 #include <TObjArray.h>
  98
  99 #include <AliPhysicsSelection.h>
 100
 101 #include <AliTriggerAnalysis.h>
 102 #include <AliLog.h>
 103
 104 #include <AliESDEvent.h>
 105 #include <AliAnalysisTaskSE.h>
 106 #include "AliAnalysisManager.h"
 107
 108 ClassImp(AliPhysicsSelection)
 109
 110 AliPhysicsSelection::AliPhysicsSelection() :
 111   AliAnalysisCuts("AliPhysicsSelection", "AliPhysicsSelection"),
 112   fCurrentRun(-1),
 113   fMC(kFALSE),
 114   fCollTrigClasses(),
 115   fBGTrigClasses(),
 116   fTriggerAnalysis(),
 117   fBackgroundIdentification(0),
 118   fHistBunchCrossing(0),
 119   fSkipTriggerClassSelection(0),
 120   fUsingCustomClasses(0),
 121   fSkipV0(0),
 122   fBIFactorA(1),
 123   fBIFactorC(1),
 124   fRatioBEEE(2),
 125   fComputeBG(0),
 126   fUseBXNumbers(0),
 127   fFillingScheme(""),
 128   fBin0CallBack(""),
 129   fBin0CallBackPointer(0)
 130 {
 131   // constructor
 132
 133   fCollTrigClasses.SetOwner(1);
 134   fBGTrigClasses.SetOwner(1);
 135   fTriggerAnalysis.SetOwner(1);
 136   fHistStatistics[0] = 0;
 137   fHistStatistics[1] = 0;
 138
 139   AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kWarning);
 140 }
 141
 142 AliPhysicsSelection::~AliPhysicsSelection()
 143 {
 144   // destructor
 145
 146   fCollTrigClasses.Delete();
 147   fBGTrigClasses.Delete();
 148   fTriggerAnalysis.Delete();
 149
 150   if (fHistStatistics[0])
 151   {
 152     delete fHistStatistics[0];
 153     fHistStatistics[0] = 0;
 154   }
 155   if (fHistStatistics[1])
 156   {
 157     delete fHistStatistics[1];
 158     fHistStatistics[1] = 0;
 159   }
 160
 161   if (fHistBunchCrossing)
 162   {
 163     delete fHistBunchCrossing;
 164     fHistBunchCrossing = 0;
 165   }
 166
 167 }
 168
 169 Bool_t AliPhysicsSelection::CheckTriggerClass(const AliESDEvent* aEsd, const char* trigger) const
 170 {
 171   // checks if the given trigger class(es) are found for the current event
 172   // format of trigger: +TRIGGER1 -TRIGGER2
 173   //   requires TRIGGER1 and rejects TRIGGER2
 174
 175   Bool_t foundBCRequirement = kFALSE;
 176   Bool_t foundCorrectBC = kFALSE;
 177
 178   TString str(trigger);
 179   TObjArray* tokens = str.Tokenize(" ");
 180
 181   for (Int_t i=0; i < tokens->GetEntries(); i++)
 182   {
 183     TString str2(((TObjString*) tokens->At(i))->String());
 184
 185     if (str2[0] == '+' || str2[0] == '-')
 186     {
 187       Bool_t flag = (str2[0] == '+');
 188
 189       str2.Remove(0, 1);
 190
 191       if (flag && !aEsd->IsTriggerClassFired(str2))
 192       {
 193         AliDebug(AliLog::kDebug, Form("Rejecting event because trigger class %s is not present", str2.Data()));
 194         delete tokens;
 195         return kFALSE;
 196       }
 197       if (!flag && aEsd->IsTriggerClassFired(str2))
 198       {
 199         AliDebug(AliLog::kDebug, Form("Rejecting event because trigger class %s is present", str2.Data()));
 200         delete tokens;
 201         return kFALSE;
 202       }
 203     }
 204     else if (str2[0] == '#')
 205     {
 206       foundBCRequirement = kTRUE;
 207
 208       str2.Remove(0, 1);
 209
 210       Int_t bcNumber = str2.Atoi();
 211       AliDebug(AliLog::kDebug, Form("Checking for bunch crossing number %d", bcNumber));
 212
 213       if (aEsd->GetBunchCrossNumber() == bcNumber)
 214       {
 215         foundCorrectBC = kTRUE;
 216         AliDebug(AliLog::kDebug, Form("Found correct bunch crossing %d", bcNumber));
 217       }
 218     }
 219     else
 220       AliFatal(Form("Invalid trigger syntax: %s", trigger));
 221   }
 222
 223   delete tokens;
 224
 225   if (foundBCRequirement && !foundCorrectBC)
 226     return kFALSE;
 227
 228   return kTRUE;
 229 }
 230
 231 Bool_t AliPhysicsSelection::IsCollisionCandidate(const AliESDEvent* aEsd)
 232 {
 233   // checks if the given event is a collision candidate
 234
 235   if (fCurrentRun != aEsd->GetRunNumber())
 236     if (!Initialize(aEsd->GetRunNumber()))
 237       AliFatal(Form("Could not initialize for run %d", aEsd->GetRunNumber()));
 238
 239   const AliESDHeader* esdHeader = aEsd->GetHeader();
 240   if (!esdHeader)
 241   {
 242     AliError("ESD Header could not be retrieved");
 243     return kFALSE;
 244   }
 245
 246   // check event type; should be PHYSICS = 7 for data and 0 for MC
 247   if (!fMC)
 248   {
 249     if (esdHeader->GetEventType() != 7)
 250       return kFALSE;
 251   }
 252   else
 253   {
 254     if (esdHeader->GetEventType() != 0)
 255       AliFatal(Form("Invalid event type for MC: %d", esdHeader->GetEventType()));
 256   }
 257
 258   Bool_t accept = kFALSE;
 259
 260   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
 261   for (Int_t i=0; i < count; i++)
 262   {
 263     const char* triggerClass = 0;
 264     if (i < fCollTrigClasses.GetEntries())
 265       triggerClass = ((TObjString*) fCollTrigClasses.At(i))->String();
 266     else
 267       triggerClass = ((TObjString*) fBGTrigClasses.At(i - fCollTrigClasses.GetEntries()))->String();
 268
 269     AliDebug(AliLog::kDebug, Form("Processing trigger class %s", triggerClass));
 270
 271     AliTriggerAnalysis* triggerAnalysis = static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i));
 272
 273     triggerAnalysis->FillTriggerClasses(aEsd);
 274
 275     if (CheckTriggerClass(aEsd, triggerClass))
 276     {
 277       triggerAnalysis->FillHistograms(aEsd);
 278
 279       Bool_t isBin0 = kFALSE;
 280       if (fBin0CallBack != "") {
 281         AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 282         if (!mgr) {
 283           AliError("Cannot get the analysis manager");
 284         }
 285         else {
 286           isBin0 = ((AliAnalysisTaskSE*)mgr->GetTask(fBin0CallBack.Data()))->IsEventInBinZero();
 287         }
 288       } else if (fBin0CallBackPointer) {
 289           isBin0 = (*fBin0CallBackPointer)(aEsd);
 290
 291       }
 292
 293
 294
 295       // hardware trigger (should only remove events for MC)
 296       // replay CINT1B hardware trigger
 297       // TODO this has to depend on the actual hardware trigger (and that depends on the run...)
 298       Int_t fastORHW = triggerAnalysis->SPDFiredChips(aEsd, 1); // SPD number of chips from trigger bits (!)
 299       Bool_t v0A       = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0A);
 300       Bool_t v0C       = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0C);
 301       Bool_t v0AHW     = fSkipV0 ? 0 :(triggerAnalysis->V0Trigger(aEsd, AliTriggerAnalysis::kASide, kTRUE) == AliTriggerAnalysis::kV0BB);// should replay hw trigger
 302       Bool_t v0CHW     = fSkipV0 ? 0 :(triggerAnalysis->V0Trigger(aEsd, AliTriggerAnalysis::kCSide, kTRUE) == AliTriggerAnalysis::kV0BB);// should replay hw trigger
 303       // offline trigger
 304       Int_t fastOROffline = triggerAnalysis->SPDFiredChips(aEsd, 0); // SPD number of chips from clusters (!)
 305       Bool_t v0ABG = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0ABG);
 306       Bool_t v0CBG = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0CBG);
 307       Bool_t v0BG = v0ABG || v0CBG;
 308
 309       // fmd
 310       Bool_t fmdA =  triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kFMDA);
 311       Bool_t fmdC =  triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kFMDC);
 312       Bool_t fmd  = fmdA || fmdC;
 313
 314       // SSD
 315       Int_t ssdClusters = triggerAnalysis->SSDClusters(aEsd);
 316
 317       // Some "macros"
 318       Bool_t mb1 = (fastOROffline > 0 || v0A || v0C) && (!v0BG);
 319       Bool_t mb1prime = (fastOROffline > 1 || (fastOROffline > 0 && (v0A || v0C)) || (v0A && v0C) ) && (!v0BG);
 320
 321       // Background rejection
 322       Bool_t bgID = ! fBackgroundIdentification->IsSelected(const_cast<AliESDEvent*> (aEsd));
 323
 324       Int_t ntrig = fastOROffline; // any 2 hits
 325       if(v0A)              ntrig += 1;
 326       if(v0C)              ntrig += 1; //v0C alone is enough
 327       if(fmd)              ntrig += 1;
 328       if(ssdClusters>1)    ntrig += 1;
 329
 330
 331       Bool_t hwTrig = fastORHW > 0 || v0AHW || v0CHW;
 332
 333       // fill statistics and return decision
 334       const Int_t nHistStat = 2;
 335       for(Int_t iHistStat = 0; iHistStat < nHistStat; iHistStat++){
 336         if (iHistStat == kStatIdxBin0 && !isBin0) continue; // skip the filling of bin0 stats if the event is not in the bin0
 337
 338         fHistStatistics[iHistStat]->Fill(kStatTriggerClass, i);
 339
 340
 341
 342         // We fill the rest only if hw trigger is ok
 343         if (!hwTrig)
 344           {
 345             AliDebug(AliLog::kDebug, "Rejecting event because hardware trigger is not fired");
 346             continue;
 347           } else {
 348           fHistStatistics[iHistStat]->Fill(kStatHWTrig, i);
 349         }
 350
 351
 352         // v0 BG stats
 353         if (v0ABG)
 354           fHistStatistics[iHistStat]->Fill(kStatV0ABG, i);
 355         if (v0CBG)
 356           fHistStatistics[iHistStat]->Fill(kStatV0CBG, i);
 357
 358         // We fill the rest only if mb1 && ! v0BG
 359         if (mb1)
 360           fHistStatistics[iHistStat]->Fill(kStatMB1, i);
 361         else continue;
 362
 363         if (mb1prime)
 364           fHistStatistics[iHistStat]->Fill(kStatMB1Prime, i);
 365
 366         if (fmd)
 367           fHistStatistics[iHistStat]->Fill(kStatFMD, i);
 368
 369         if(ssdClusters>1)
 370           fHistStatistics[iHistStat]->Fill(kStatSSD1, i);
 371
 372         if(ntrig >= 2 && !v0BG)
 373           fHistStatistics[iHistStat]->Fill(kStatAny2Hits, i);
 374
 375         if (fastOROffline > 0)
 376           fHistStatistics[iHistStat]->Fill(kStatFO1, i);
 377         if (fastOROffline > 1)
 378           fHistStatistics[iHistStat]->Fill(kStatFO2, i);
 379
 380         if (v0A)
 381           fHistStatistics[iHistStat]->Fill(kStatV0A, i);
 382         if (v0C)
 383           fHistStatistics[iHistStat]->Fill(kStatV0C, i);
 384
 385         //       if (fastOROffline > 1 && !v0BG)
 386         //         fHistStatistics[iHistStat]->Fill(kStatFO2NoBG, i);
 387
 388         if (fastOROffline > 0 && (v0A || v0C) && !v0BG)
 389           fHistStatistics[iHistStat]->Fill(kStatFO1AndV0, i);
 390
 391         if (v0A && v0C && !v0BG && !bgID)
 392           fHistStatistics[iHistStat]->Fill(kStatV0, i);
 393
 394
 395         if ( mb1 )
 396
 397           {
 398             if (!v0BG || fSkipV0)
 399               {
 400                 if (!v0BG) fHistStatistics[iHistStat]->Fill(kStatOffline, i);
 401
 402                 if (fBackgroundIdentification && bgID)
 403                   {
 404                     AliDebug(AliLog::kDebug, "Rejecting event because of background identification");
 405                     fHistStatistics[iHistStat]->Fill(kStatBG, i);
 406                   }
 407                 else
 408                   {
 409                     AliDebug(AliLog::kDebug, "Accepted event for histograms");
 410
 411                     fHistStatistics[iHistStat]->Fill(kStatAccepted, i);
 412                     if(iHistStat == kStatIdxAll) fHistBunchCrossing->Fill(aEsd->GetBunchCrossNumber(), i); // Fill only for all (avoid double counting)
 413                     if((i < fCollTrigClasses.GetEntries() || fSkipTriggerClassSelection) && (iHistStat==kStatIdxAll))
 414                       accept = kTRUE; // only set for "all" (should not really matter)
 415                   }
 416               }
 417             else
 418               AliDebug(AliLog::kDebug, "Rejecting event because of V0 BG flag");
 419           }
 420         else
 421           AliDebug(AliLog::kDebug, "Rejecting event because trigger condition is not fulfilled");
 422       }
 423     }
 424   }
 425
 426   if (accept)
 427     AliDebug(AliLog::kDebug, "Accepted event as collision candidate");
 428
 429   return accept;
 430 }
 431
 432 Int_t AliPhysicsSelection::GetTriggerScheme(UInt_t runNumber) const
 433 {
 434   // returns the current trigger scheme (classes that are accepted/rejected)
 435
 436   if (fMC)
 437     return 0;
 438
 439   // TODO dependent on run number
 440
 441   switch (runNumber)
 442   {
 443     // CSMBB triggers
 444     case 104044:
 445     case 105054:
 446     case 105057:
 447       return 2;
 448   }
 449
 450   // default: CINT1 suite
 451   return 1;
 452 }
 453
 454 const char * AliPhysicsSelection::GetFillingScheme(UInt_t runNumber)  {
 455
 456   if(fMC) return "MC";
 457
 458   if      (runNumber >= 104065 && runNumber <= 104160) {
 459     return "4x4a";
 460   }
 461   else if (runNumber >= 104315 && runNumber <= 104321) {
 462     return "4x4a*";
 463   }
 464   else if (runNumber >= 104792 && runNumber <= 104803) {
 465     return "4x4b";
 466   }
 467   else if (runNumber >= 104824 && runNumber <= 104892) {
 468     return "4x4c";
 469   }
 470   else if (runNumber == 105143 || runNumber == 105160) {
 471     return "16x16a";
 472   }
 473   else if (runNumber >= 105256 && runNumber <= 105268) {
 474     return "4x4c";
 475   }
 476   else {
 477     AliError(Form("Unknown filling scheme (run %d)", runNumber));
 478   }
 479
 480   return "Unknown";
 481 }
 482
 483 Int_t AliPhysicsSelection::GetRatioBBBE(Int_t runNumber) {
 484
 485
 486   if(fMC) return 1;
 487
 488   if (runNumber == 105143 || runNumber == 105160) {
 489     return 8;
 490   }
 491   else if (fComputeBG &&
 492            !(runNumber >= 105256 && runNumber <= 105268) &&
 493            !(runNumber >= 104065 && runNumber <= 104160) &&
 494            !(runNumber >= 104315 && runNumber <= 104321) &&
 495            !(runNumber >= 104792 && runNumber <= 104803) &&
 496            !(runNumber >= 104824 && runNumber <= 104892)
 497            ){
 498
 499     AliError(Form("Unknown run %d, assuming ratio BE/EE = 2",runNumber));
 500
 501   }
 502
 503   return 2;
 504 }
 505
 506
 507 const char * AliPhysicsSelection::GetBXIDs(UInt_t runNumber, const char * trigger)  {
 508
 509   if (!fUseBXNumbers || fMC) return "";
 510
 511   if      (runNumber >= 104065 && runNumber <= 104160) {
 512     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #2128 #3019";
 513     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #346 #3465";
 514     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1234 #1680";
 515     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
 516     else AliError(Form("Unknown trigger: %s", trigger));
 517   }
 518   else if (runNumber >= 104315 && runNumber <= 104321) {
 519     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #2000 #2891";
 520     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #218 #3337";
 521     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1106 #1552";
 522     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
 523     else AliError(Form("Unknown trigger: %s", trigger));
 524   }
 525   else if (runNumber >= 104792 && runNumber <= 104803) {
 526     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #2228 #3119";
 527     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2554 #446";
 528     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1334 #769";
 529     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
 530     else AliError(Form("Unknown trigger: %s", trigger));
 531   }
 532   else if (runNumber >= 104824 && runNumber <= 104892) {
 533     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #3119 #769";
 534     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2554 #446";
 535     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1334 #2228";
 536     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
 537     else AliError(Form("Unknown trigger: %s", trigger));
 538   }
 539   else if (runNumber == 105143 || runNumber == 105160) {
 540     fRatioBEEE = 8;
 541     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #1337 #1418 #2228 #2309 #3119 #3200 #446 #527";
 542     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #1580  #1742  #1904  #2066  #2630  #2792  #2954  #3362";
 543     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "  #845  #1007  #1169   #1577 #3359 #3521 #119  #281 ";
 544     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
 545     else AliError(Form("Unknown trigger: %s", trigger));
 546   }
 547   else if (runNumber >= 105256 && runNumber <= 105268) {
 548     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #3019 #669";
 549     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2454 #346";
 550     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1234 #2128";
 551     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
 552     else AliError(Form("Unknown trigger: %s", trigger));
 553   }
 554   else {
 555     AliError(Form("Unknown run %d, using all BXs!",runNumber));
 556   }
 557
 558   return "";
 559 }
 560
 561 Bool_t AliPhysicsSelection::Initialize(Int_t runNumber)
 562 {
 563   // initializes the object for the given run
 564   // TODO having the run number here and parameters hardcoded is clearly temporary, a way needs to be found to have a CDB-like configuration also for analysis
 565
 566   Bool_t oldStatus = TH1::AddDirectoryStatus();
 567   TH1::AddDirectory(kFALSE);
 568
 569   if(!fBin0CallBack)
 570     AliError("Bin0 Callback not set: will not fill the statistics for the bin 0");
 571
 572   if (fMC) {
 573     // ovverride BX and bg options in case of MC
 574     fComputeBG    = kFALSE;
 575     fUseBXNumbers = kFALSE;
 576   }
 577
 578   Int_t triggerScheme = GetTriggerScheme(runNumber);
 579   if (!fUsingCustomClasses && fCurrentRun != -1 && triggerScheme != GetTriggerScheme(fCurrentRun))
 580     AliFatal("Processing several runs with different trigger schemes is not supported");
 581
 582   if(fComputeBG && fCurrentRun != -1 && fCurrentRun != runNumber)
 583     AliFatal("Cannot process several runs because BG computation is requested");
 584
 585   if(fComputeBG && !fUseBXNumbers)
 586     AliFatal("Cannot compute BG id BX numbers are not used");
 587
 588   if(fUseBXNumbers && fFillingScheme != "" && fFillingScheme != GetFillingScheme(runNumber))
 589     AliFatal("Cannot process runs with different filling scheme if usage of BX numbers is requested");
 590
 591   fRatioBEEE          = GetRatioBBBE(runNumber);
 592   fFillingScheme      = GetFillingScheme(runNumber);
 593
 594   if(fComputeBG) SetBIFactors(runNumber);
 595
 596   AliInfo(Form("Initializing for run %d", runNumber));
 597
 598   // initialize first time?
 599   if (fCurrentRun == -1)
 600   {
 601     if (fUsingCustomClasses) {
 602       AliInfo("Using user-provided trigger classes");
 603     } else {
 604       switch (triggerScheme)
 605       {
 606       case 0:
 607         fCollTrigClasses.Add(new TObjString(""));
 608         break;
 609
 610       case 1:
 611         { // need a new scope to avoid cross-initialization errors
 612           TObjString * cint1b = new TObjString(Form("%s%s","+CINT1B-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1B-ABCE-NOPF-ALL")));
 613           TObjString * cint1a = new TObjString(Form("%s%s","+CINT1A-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1A-ABCE-NOPF-ALL")));
 614           TObjString * cint1c = new TObjString(Form("%s%s","+CINT1C-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1C-ABCE-NOPF-ALL")));
 615           TObjString * cint1e = new TObjString(Form("%s%s","+CINT1-E-NOPF-ALL",      GetBXIDs(runNumber,"CINT1-E-NOPF-ALL"))    );
 616           //
 617           fCollTrigClasses.Add(cint1b);
 618           fBGTrigClasses.Add(cint1a);
 619           fBGTrigClasses.Add(cint1c);
 620           fBGTrigClasses.Add(cint1e);
 621         }
 622         break;
 623
 624       case 2:
 625         fCollTrigClasses.Add(new TObjString("+CSMBB-ABCE-NOPF-ALL"));
 626         fBGTrigClasses.Add(new TObjString("+CSMBA-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL"));
 627         fBGTrigClasses.Add(new TObjString("+CSMBC-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL"));
 628         break;
 629
 630       default:
 631         AliFatal(Form("Unsupported trigger scheme %d", triggerScheme));
 632       }
 633     }
 634
 635     Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
 636
 637     for (Int_t i=0; i<count; i++)
 638     {
 639       AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis;
 640       triggerAnalysis->SetAnalyzeMC(fMC);
 641       triggerAnalysis->EnableHistograms();
 642       triggerAnalysis->SetSPDGFOThreshhold(1);
 643       fTriggerAnalysis.Add(triggerAnalysis);
 644     }
 645
 646     // TODO: shall I really delete this?
 647     if (fHistStatistics[0])
 648       delete fHistStatistics[0];
 649     if (fHistStatistics[1])
 650       delete fHistStatistics[1];
 651
 652     fHistStatistics[kStatIdxBin0] = BookHistStatistics("_Bin0");
 653     fHistStatistics[kStatIdxAll]  = BookHistStatistics("");
 654
 655     if (fHistBunchCrossing)
 656       delete fHistBunchCrossing;
 657
 658     fHistBunchCrossing = new TH2F("fHistBunchCrossing", "fHistBunchCrossing;bunch crossing number;", 4000, -0.5, 3999.5,  count, -0.5, -0.5 + count);
 659
 660     Int_t n = 1;
 661     for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
 662     {
 663       fHistBunchCrossing->GetYaxis()->SetBinLabel(n, ((TObjString*) fCollTrigClasses.At(i))->String());
 664       n++;
 665     }
 666     for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
 667     {
 668       fHistBunchCrossing->GetYaxis()->SetBinLabel(n, ((TObjString*) fBGTrigClasses.At(i))->String());
 669       n++;
 670     }
 671
 672
 673
 674   }
 675
 676   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
 677   for (Int_t i=0; i<count; i++)
 678   {
 679     AliTriggerAnalysis* triggerAnalysis = static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i));
 680
 681     switch (runNumber)
 682     {
 683       case 104315:
 684       case 104316:
 685       case 104320:
 686       case 104321:
 687       case 104439:
 688         triggerAnalysis->SetV0TimeOffset(7.5);
 689         break;
 690       default:
 691         triggerAnalysis->SetV0TimeOffset(0);
 692     }
 693   }
 694
 695   fCurrentRun = runNumber;
 696
 697   TH1::AddDirectory(oldStatus);
 698
 699   return kTRUE;
 700 }
 701
 702 TH2F * AliPhysicsSelection::BookHistStatistics(const char * tag) {
 703
 704     // add 6 rows to count for the estimate of good, accidentals and
 705     // BG and the ratio of BG and accidentals to total +ratio goot to
 706     // first col + 2 for error on good.
 707
 708   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
 709 #ifdef VERBOSE_STAT
 710   Int_t extrarows = fComputeBG ? 8 : 0;
 711 #else
 712   Int_t extrarows = fComputeBG ? 3 : 0;
 713 #endif
 714   TH2F * h = new TH2F(Form("fHistStatistics%s",tag), Form("fHistStatistics - %s ;;",tag), kStatAccepted, 0.5, kStatAccepted+0.5, count+extrarows, -0.5, -0.5 + count+extrarows);
 715
 716   h->GetXaxis()->SetBinLabel(kStatTriggerClass,  "Trigger class");
 717   h->GetXaxis()->SetBinLabel(kStatHWTrig,        "Hardware trigger");
 718   h->GetXaxis()->SetBinLabel(kStatFO1,           "FO >= 1");
 719   h->GetXaxis()->SetBinLabel(kStatFO2,           "FO >= 2");
 720   h->GetXaxis()->SetBinLabel(kStatV0A,           "V0A");
 721   h->GetXaxis()->SetBinLabel(kStatV0C,           "V0C");
 722   h->GetXaxis()->SetBinLabel(kStatFMD,           "FMD");
 723   h->GetXaxis()->SetBinLabel(kStatSSD1,          "SSD >= 2");
 724   h->GetXaxis()->SetBinLabel(kStatV0ABG,         "V0A BG");
 725   h->GetXaxis()->SetBinLabel(kStatV0CBG,         "V0C BG");
 726   h->GetXaxis()->SetBinLabel(kStatMB1,           "(FO >= 1 | V0A | V0C) & !V0 BG");
 727   h->GetXaxis()->SetBinLabel(kStatMB1Prime,      "(FO >= 2 | (FO >= 1 & (V0A | V0C)) | (V0A &v0C) ) & !V0 BG");
 728   h->GetXaxis()->SetBinLabel(kStatFO1AndV0,      "FO >= 1 & (V0A | V0C) & !V0 BG");
 729   h->GetXaxis()->SetBinLabel(kStatV0,            "V0A & V0C & !V0 BG & !BG ID");
 730   h->GetXaxis()->SetBinLabel(kStatOffline,       "Offline Trigger");
 731   h->GetXaxis()->SetBinLabel(kStatAny2Hits,      "2 Hits & !V0 BG");
 732   h->GetXaxis()->SetBinLabel(kStatBG,            "Background identification");
 733   h->GetXaxis()->SetBinLabel(kStatAccepted,      "Accepted");
 734
 735   Int_t n = 1;
 736   for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
 737     {
 738       h->GetYaxis()->SetBinLabel(n, ((TObjString*) fCollTrigClasses.At(i))->String());
 739       n++;
 740     }
 741   for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
 742     {
 743       h->GetYaxis()->SetBinLabel(n, ((TObjString*) fBGTrigClasses.At(i))->String());
 744       n++;
 745     }
 746
 747   if(fComputeBG) {
 748     h->GetYaxis()->SetBinLabel(n++, "BG (A+C)");
 749     h->GetYaxis()->SetBinLabel(n++, "ACC");
 750 #ifdef VERBOSE_STAT
 751     h->GetYaxis()->SetBinLabel(n++, "BG (A+C) %  (rel. to CINT1B)");
 752     h->GetYaxis()->SetBinLabel(n++, "ACC % (rel. to CINT1B)");
 753     h->GetYaxis()->SetBinLabel(n++, "ERR GOOD %");
 754     h->GetYaxis()->SetBinLabel(n++, "GOOD % (rel. to 1st col)");
 755     h->GetYaxis()->SetBinLabel(n++, "ERR GOOD");
 756 #endif
 757     h->GetYaxis()->SetBinLabel(n++, "GOOD");
 758   }
 759
 760   return h;
 761 }
 762
 763 void AliPhysicsSelection::Print(Option_t* /* option */) const
 764 {
 765   // print the configuration
 766
 767   Printf("Configuration initialized for run %d (MC: %d):", fCurrentRun, fMC);
 768
 769   Printf("Collision trigger classes:");
 770   for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
 771     Printf("%s", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
 772
 773   Printf("Background trigger classes:");
 774   for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
 775     Printf("%s", ((TObjString*) fBGTrigClasses.At(i))->String().Data());
 776
 777   AliTriggerAnalysis* triggerAnalysis = dynamic_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(0));
 778
 779   if (triggerAnalysis)
 780   {
 781     if (triggerAnalysis->GetV0TimeOffset() > 0)
 782       Printf("V0 time offset active: %.2f ns", triggerAnalysis->GetV0TimeOffset());
 783
 784     Printf("\nTotal available events:");
 785
 786     triggerAnalysis->PrintTriggerClasses();
 787   }
 788
 789   if (fHistStatistics[kStatIdxAll] && fCollTrigClasses.GetEntries() > 0)
 790   {
 791     Printf("\nSelection statistics for first collision trigger (%s):", ((TObjString*) fCollTrigClasses.First())->String().Data());
 792
 793     Printf("Total events with correct trigger class: %d", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(1, 1));
 794     Printf("Selected collision candidates: %d", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(fHistStatistics[kStatIdxAll]->GetXaxis()->FindBin("Accepted"), 1));
 795   }
 796
 797   if (fHistBunchCrossing)
 798   {
 799     Printf("\nBunch crossing statistics:");
 800
 801     for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
 802     {
 803       TString str;
 804       str.Form("Trigger %s has accepted events in the bunch crossings: ", fHistBunchCrossing->GetYaxis()->GetBinLabel(i));
 805
 806       for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
 807         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
 808           str += Form("%d, ", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
 809
 810       Printf("%s", str.Data());
 811     }
 812
 813     for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
 814     {
 815       Int_t count = 0;
 816       for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
 817       {
 818         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
 819           count++;
 820       }
 821       if (count > 1)
 822         Printf("WARNING: Bunch crossing %d has more than one trigger class active. Check BPTX functioning for this run!", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
 823     }
 824   }
 825
 826   if (fUsingCustomClasses)
 827     Printf("WARNING: Using custom trigger classes!");
 828   if (fSkipTriggerClassSelection)
 829     Printf("WARNING: Skipping trigger class selection!");
 830   if (fSkipV0)
 831     Printf("WARNING: Ignoring V0 information in selection");
 832   if(!fBin0CallBack)
 833     Printf("WARNING: Callback not set: will not fill the statistics for the bin 0");
 834
 835 }
 836
 837 Long64_t AliPhysicsSelection::Merge(TCollection* list)
 838 {
 839   // Merge a list of AliMultiplicityCorrection objects with this (needed for
 840   // PROOF).
 841   // Returns the number of merged objects (including this).
 842
 843   if (!list)
 844     return 0;
 845
 846   if (list->IsEmpty())
 847     return 1;
 848
 849   TIterator* iter = list->MakeIterator();
 850   TObject* obj;
 851
 852   // collections of all histograms
 853   const Int_t nHists = 9;
 854   TList collections[nHists];
 855
 856   Int_t count = 0;
 857   while ((obj = iter->Next())) {
 858
 859     AliPhysicsSelection* entry = dynamic_cast<AliPhysicsSelection*> (obj);
 860     if (entry == 0)
 861       continue;
 862
 863     collections[0].Add(&(entry->fTriggerAnalysis));
 864     if (entry->fHistStatistics[0])
 865       collections[1].Add(entry->fHistStatistics[0]);
 866     if (entry->fHistStatistics[1])
 867       collections[2].Add(entry->fHistStatistics[1]);
 868     if (entry->fHistBunchCrossing)
 869       collections[3].Add(entry->fHistBunchCrossing);
 870     if (entry->fBackgroundIdentification)
 871       collections[4].Add(entry->fBackgroundIdentification);
 872
 873     count++;
 874   }
 875
 876   fTriggerAnalysis.Merge(&collections[0]);
 877   if (fHistStatistics[0])
 878     fHistStatistics[0]->Merge(&collections[1]);
 879   if (fHistStatistics[1])
 880     fHistStatistics[1]->Merge(&collections[2]);
 881   if (fHistBunchCrossing)
 882     fHistBunchCrossing->Merge(&collections[3]);
 883   if (fBackgroundIdentification)
 884     fBackgroundIdentification->Merge(&collections[4]);
 885
 886   delete iter;
 887
 888   return count+1;
 889 }
 890
 891 void AliPhysicsSelection::SaveHistograms(const char* folder) const
 892 {
 893   // write histograms to current directory
 894
 895   if (!fHistStatistics[0] || !fHistStatistics[1])
 896     return;
 897
 898   if (folder)
 899   {
 900     gDirectory->mkdir(folder);
 901     gDirectory->cd(folder);
 902   }
 903
 904
 905   // Fill the last rows of fHistStatistics before saving
 906   if (fComputeBG) {
 907     Int_t triggerScheme = GetTriggerScheme(UInt_t(fCurrentRun));
 908     if(triggerScheme != 1){
 909       AliWarning("BG estimate only supported for trigger scheme \"1\" (CINT1 suite)");
 910     } else {
 911       Int_t nHistStat = 2;
 912       // TODO: get number of rows in a more flexible way
 913       // 1. loop over all cols
 914
 915       for(Int_t iHistStat = 0; iHistStat < nHistStat; iHistStat++){
 916
 917         Int_t ncol = fHistStatistics[iHistStat]->GetNbinsX();
 918         Float_t good1 = 0;
 919         for(Int_t icol = 1; icol <= ncol; icol++) {
 920           Int_t cint1B = (Int_t) fHistStatistics[iHistStat]->GetBinContent(icol,1);
 921           Int_t cint1A = (Int_t) fHistStatistics[iHistStat]->GetBinContent(icol,2);
 922           Int_t cint1C = (Int_t) fHistStatistics[iHistStat]->GetBinContent(icol,3);
 923           Int_t cint1E = (Int_t) fHistStatistics[iHistStat]->GetBinContent(icol,4);
 924
 925           if (cint1B>0) {
 926             Int_t acc  = fRatioBEEE*cint1E;
 927             //      Int_t bg   = cint1A + cint1C - 2*acc;
 928             Float_t bg   = fBIFactorA*(cint1A-acc) + fBIFactorC*(cint1C-acc) ;
 929             Float_t good = Float_t(cint1B) - bg - acc;
 930             if (icol ==1) good1 = good;
 931             //      Float_t errGood     = TMath::Sqrt(2*(cint1A+cint1C+cint1E));// Error on the number of goods assuming only bg fluctuates
 932             //      DeltaG^2 = B + FA^2 A + FC^2 C + Ratio^2 (FA+FC-1)^2 E.
 933             Float_t errGood     = TMath::Sqrt( cint1B +
 934                                                fBIFactorA*fBIFactorA*cint1A +
 935                                                fBIFactorC*fBIFactorC*cint1C +
 936                                                fRatioBEEE * fRatioBEEE *
 937                                                (fBIFactorA + fBIFactorC - 1)*(fBIFactorA + fBIFactorC - 1)*cint1E);
 938             Float_t errBG = TMath::Sqrt(fBIFactorA*fBIFactorA*cint1A+
 939                                         fBIFactorC*fBIFactorC*cint1C+
 940                                         fRatioBEEE*fRatioBEEE*(fBIFactorA+fBIFactorC)*(fBIFactorA+fBIFactorC)*cint1E);
 941
 942
 943             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowBG,bg);
 944             fHistStatistics[iHistStat]->SetBinError  (icol,kStatRowBG,errBG);
 945             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowAcc,acc);
 946             fHistStatistics[iHistStat]->SetBinError  (icol,kStatRowAcc,TMath::Sqrt(fRatioBEEE*fRatioBEEE*cint1E));
 947             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowGood,good);
 948             fHistStatistics[iHistStat]->SetBinError  (icol,kStatRowGood,errGood);
 949
 950 #ifdef VERBOSE_STAT
 951             Float_t accFrac   = Float_t(acc) / cint1B  *100;
 952             Float_t bgFrac    = Float_t(bg)  / cint1B  *100;
 953             Float_t goodFrac  = Float_t(good)  / good1 *100;
 954             Float_t errGoodFrac = errGood/good1 * 100;
 955             Float_t errFracBG = bg > 0 ? TMath::Sqrt(errBG/bg + 1/TMath::Sqrt(cint1B))*bgFrac : 0;
 956             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowBGFrac,bgFrac);
 957             fHistStatistics[iHistStat]->SetBinError  (icol,kStatRowBGFrac,errFracBG);
 958             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowAccFrac,accFrac);
 959             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowAccFrac,errGoodFrac);
 960             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowGoodFrac,goodFrac);
 961             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowErrGood,errGood);
 962 #endif
 963           }
 964         }
 965       }
 966     }
 967   }
 968
 969   fHistStatistics[0]->Write();
 970   fHistStatistics[1]->Write();
 971   fHistBunchCrossing->Write();
 972
 973   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
 974   for (Int_t i=0; i < count; i++)
 975   {
 976     TString triggerClass = "trigger_histograms_";
 977     if (i < fCollTrigClasses.GetEntries())
 978       triggerClass += ((TObjString*) fCollTrigClasses.At(i))->String();
 979     else
 980       triggerClass += ((TObjString*) fBGTrigClasses.At(i - fCollTrigClasses.GetEntries()))->String();
 981
 982     gDirectory->mkdir(triggerClass);
 983     gDirectory->cd(triggerClass);
 984
 985     static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i))->SaveHistograms();
 986
 987     gDirectory->cd("..");
 988   }
 989
 990   if (fBackgroundIdentification)
 991   {
 992     gDirectory->mkdir("background_identification");
 993     gDirectory->cd("background_identification");
 994
 995     fBackgroundIdentification->GetOutput()->Write();
 996
 997     gDirectory->cd("..");
 998   }
 999
1000   if (folder)
1001     gDirectory->cd("..");
1002 }
1003
1004 void AliPhysicsSelection::SetBIFactors(Int_t run) {
1005
1006   switch(run) {
1007   case 104155:
1008     fBIFactorA = 0.961912722908;
1009     fBIFactorC = 1.04992336081;
1010     break;
1011   case 104157:
1012     fBIFactorA = 0.947312854998;
1013     fBIFactorC = 1.01599706417;
1014     break;
1015   case 104159:
1016     fBIFactorA = 0.93659320151;
1017     fBIFactorC = 0.98580804207;
1018     break;
1019   case 104160:
1020     fBIFactorA = 0.929664189926;
1021     fBIFactorC = 0.963467679851;
1022     break;
1023   case 104315:
1024     fBIFactorA = 1.08939104979;
1025     fBIFactorC = 0.931113921925;
1026     break;
1027   case 104316:
1028     fBIFactorA = 1.08351880974;
1029     fBIFactorC = 0.916068345845;
1030     break;
1031   case 104320:
1032     fBIFactorA = 1.07669281245;
1033     fBIFactorC = 0.876818744763;
1034     break;
1035   case 104321:
1036     fBIFactorA = 1.00971079602;
1037     fBIFactorC = 0.773781299076;
1038     break;
1039   case 104792:
1040     fBIFactorA = 0.787215863962;
1041     fBIFactorC = 0.778253173071;
1042     break;
1043   case 104793:
1044     fBIFactorA = 0.692211363661;
1045     fBIFactorC = 0.733152456667;
1046     break;
1047   case 104799:
1048     fBIFactorA = 1.04027825161;
1049     fBIFactorC = 1.00530825942;
1050     break;
1051   case 104800:
1052     fBIFactorA = 1.05309910671;
1053     fBIFactorC = 1.00376801855;
1054     break;
1055   case 104801:
1056     fBIFactorA = 1.0531231922;
1057     fBIFactorC = 0.992439666758;
1058     break;
1059   case 104802:
1060     fBIFactorA = 1.04191478134;
1061     fBIFactorC = 0.979368585208;
1062     break;
1063   case 104803:
1064     fBIFactorA = 1.03121314094;
1065     fBIFactorC = 0.973379962609;
1066     break;
1067   case 104824:
1068     fBIFactorA = 0.969945926722;
1069     fBIFactorC = 0.39549745806;
1070     break;
1071   case 104825:
1072     fBIFactorA = 0.968627213937;
1073     fBIFactorC = 0.310100412205;
1074     break;
1075   case 104841:
1076     fBIFactorA = 0.991601393212;
1077     fBIFactorC = 0.83762204722;
1078     break;
1079   case 104845:
1080     fBIFactorA = 0.98040863886;
1081     fBIFactorC = 0.694824205793;
1082     break;
1083   case 104867:
1084     fBIFactorA = 1.10646173412;
1085     fBIFactorC = 0.841407246916;
1086     break;
1087   case 104876:
1088     fBIFactorA = 1.12063452421;
1089     fBIFactorC = 0.78726542895;
1090     break;
1091   case 104890:
1092     fBIFactorA = 1.02346137453;
1093     fBIFactorC = 1.03355663595;
1094     break;
1095   case 104892:
1096     fBIFactorA = 1.05406025913;
1097     fBIFactorC = 1.00029166135;
1098     break;
1099   case 105143:
1100     fBIFactorA = 0.947343384349;
1101     fBIFactorC = 0.972637444408;
1102     break;
1103   case 105160:
1104     fBIFactorA = 0.908854622177;
1105     fBIFactorC = 0.958851103977;
1106     break;
1107   case 105256:
1108     fBIFactorA = 0.810076150206;
1109     fBIFactorC = 0.884663561883;
1110     break;
1111   case 105257:
1112     fBIFactorA = 0.80974912303;
1113     fBIFactorC = 0.878859123479;
1114     break;
1115   case 105268:
1116     fBIFactorA = 0.809052110679;
1117     fBIFactorC = 0.87233890989;
1118     break;
1119   default:
1120     fBIFactorA = 1;
1121     fBIFactorC = 1;
1122   }
1123
1124
1125 }