Modifications to the computation of residul background in the physics selection,...
[u/mrichter/AliRoot.git] / 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 // The BX ids corresponding to real bunches crossings p2 are
42 // automatically selected. You cannot process runs with different
43 // filling schemes if this option is set. If you want to disable this,
44 // use: 
45 //   fPhysicsSelection->SetUseBXNumbers(0);
46 //
47 //
48 // If you are analizing muons and you want to keep the muon triggers
49 // besides the CINT1B you can set:
50 //   fPhysicsSelection->SetUseMuonTriggers();
51 //
52 //
53 // To compute the Background automatically using the control triggers
54 // use: 
55 //   fPhysicsSelection->SetComputeBG();
56 // this will show the value of the Beam Gas, accidentals and good
57 // events as additional rows in the statistic tables, but it will NOT
58 // subtract the background automatically.
59 // This option will only work for runs taken with the CINT1
60 // suite. This options enables automatically also the usage of BX
61 // numbers. You can only process one run at a time if you require this
62 // options, because it uses the bunch intensity estimated run by run.
63 // 
64 // The BG will usually be more important in the so-called "bin 0": the
65 // class can also compute the statistics table for events in this
66 // bin. Since the definition of bin 0 may in general change from
67 // analysis to analysis, the user needs to provide a callback
68 // implementing the definition of bin zero. The callback should be
69 // implemented as a method in the analysis task and should override
70 // the IsEventInBinZero method of AliAnalysisTaskSE, and should thus
71 // have the the following prototype:
72 //   Bool_t IsEventInBinZero(); 
73 // It should return true if the event is in the bin 0 and it is set by
74 // passing to the physics selection the NAME of the task where the
75 // callback is implemented: 
76 //   fPhysicsSelection->SetBin0Callback("MyTask").
77 //
78 //
79 // Usually the class selects the trigger scheme by itself depending on the run number.
80 // Nevertheless, you can do that manually by calling AddCollisionTriggerClass() and AddBGTriggerClass()
81 // Example:
82 // To define the class CINT1B-ABCE-NOPF-ALL as collision trigger (those will be accepted as  
83 // collision candidates when they pass the selection):
84 //   AddCollisionTriggerClass("+CINT1B-ABCE-NOPF-ALL #769 #3119");
85 // To select on bunch crossing IDs in addition, use:
86 //   AddCollisionTriggerClass("+CINT1B-ABCE-NOPF-ALL #769 #3119");
87 // To define the class CINT1A-ABCE-NOPF-ALL as a background trigger (those will only be counted
88 // for the control histograms):
89 //   AddBGTriggerClass("+CINT1A-ABCE-NOPF-ALL");
90 // You can also specify more than one trigger class in a string or you can require that some are *not*
91 // present. The following line would require CSMBA-ABCE-NOPF-ALL, but CSMBB-ABCE-NOPF-ALL is not allowed
92 // to be present:
93 //   AddBGTriggerClass("+CSMBA-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL");
94 //
95 //   Origin: Jan Fiete Grosse-Oetringhaus, CERN 
96 //           Michele Floris, CERN
97 //-------------------------------------------------------------------------
98
99 #include <Riostream.h>
100 #include <TH1F.h>
101 #include <TH2F.h>
102 #include <TList.h>
103 #include <TIterator.h>
104 #include <TDirectory.h>
105 #include <TObjArray.h>
106
107 #include <AliPhysicsSelection.h>
108
109 #include <AliTriggerAnalysis.h>
110 #include <AliLog.h>
111
112 #include <AliESDEvent.h>
113 #include <AliAnalysisTaskSE.h>
114 #include "AliAnalysisManager.h"
115 #include "TPRegexp.h"
116
117 ClassImp(AliPhysicsSelection)
118
119 AliPhysicsSelection::AliPhysicsSelection() :
120   AliAnalysisCuts("AliPhysicsSelection", "AliPhysicsSelection"),
121   fCurrentRun(-1),
122   fMC(kFALSE),
123   fCollTrigClasses(),
124   fBGTrigClasses(),
125   fTriggerAnalysis(),
126   fBackgroundIdentification(0),
127   fHistBunchCrossing(0),
128   fHistTriggerPattern(0),
129   fSkipTriggerClassSelection(0),
130   fUsingCustomClasses(0),
131   fSkipV0(0),
132   fBIFactorA(1),
133   fBIFactorC(1),
134   fBIFactorAC(1), 
135   fComputeBG(0),
136   fBGStatOffset(0),
137   fUseBXNumbers(1),
138   fUseMuonTriggers(0),
139   fFillingScheme(""),
140   fBin0CallBack(""),
141   fBin0CallBackPointer(0)
142 {
143   // constructor
144   
145   fCollTrigClasses.SetOwner(1);
146   fBGTrigClasses.SetOwner(1);
147   fTriggerAnalysis.SetOwner(1);
148   fHistStatistics[0] = 0;
149   fHistStatistics[1] = 0;
150   
151   AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kWarning);
152 }
153     
154 AliPhysicsSelection::~AliPhysicsSelection()
155 {
156   // destructor
157
158   fCollTrigClasses.Delete();
159   fBGTrigClasses.Delete();
160   fTriggerAnalysis.Delete();
161
162   if (fHistStatistics[0])
163   {
164     delete fHistStatistics[0];
165     fHistStatistics[0] = 0;
166   }
167   if (fHistStatistics[1])
168   {
169     delete fHistStatistics[1];
170     fHistStatistics[1] = 0;
171   }
172   
173   if (fHistBunchCrossing)
174   {
175     delete fHistBunchCrossing;
176     fHistBunchCrossing = 0;
177   }
178   if (fHistTriggerPattern)
179   {
180     delete fHistTriggerPattern;
181     fHistTriggerPattern = 0;
182   }
183
184 }
185
186 UInt_t AliPhysicsSelection::CheckTriggerClass(const AliESDEvent* aEsd, const char* trigger) const
187 {
188   // checks if the given trigger class(es) are found for the current event
189   // format of trigger: +TRIGGER1 -TRIGGER2 [#XXX] [&YY]
190   //   requires TRIGGER1 and rejects TRIGGER2
191   //   in bunch crossing XXX
192   //   if successful, YY is returned (for association between entry in fCollTrigClasses and AliVEvent::EOfflineTriggerTypes)
193   
194   Bool_t foundBCRequirement = kFALSE;
195   Bool_t foundCorrectBC = kFALSE;
196   
197   UInt_t returnCode = AliVEvent::kUserDefined;
198   
199   TString str(trigger);
200   TObjArray* tokens = str.Tokenize(" ");
201   
202   for (Int_t i=0; i < tokens->GetEntries(); i++)
203   {
204     TString str2(((TObjString*) tokens->At(i))->String());
205     
206     if (str2[0] == '+' || str2[0] == '-')
207     {
208       Bool_t flag = (str2[0] == '+');
209       
210       str2.Remove(0, 1);
211       
212       if (flag && !aEsd->IsTriggerClassFired(str2))
213       {
214         AliDebug(AliLog::kDebug, Form("Rejecting event because trigger class %s is not present", str2.Data()));
215         delete tokens;
216         return kFALSE;
217       }
218       if (!flag && aEsd->IsTriggerClassFired(str2))
219       {
220         AliDebug(AliLog::kDebug, Form("Rejecting event because trigger class %s is present", str2.Data()));
221         delete tokens;
222         return kFALSE;
223       }
224     }
225     else if (str2[0] == '#')
226     {
227       foundBCRequirement = kTRUE;
228     
229       str2.Remove(0, 1);
230       
231       Int_t bcNumber = str2.Atoi();
232       AliDebug(AliLog::kDebug, Form("Checking for bunch crossing number %d", bcNumber));
233       
234       if (aEsd->GetBunchCrossNumber() == bcNumber)
235       {
236         foundCorrectBC = kTRUE;
237         AliDebug(AliLog::kDebug, Form("Found correct bunch crossing %d", bcNumber));
238       }
239     }
240     else if (str2[0] == '&' && !fUsingCustomClasses)
241     {
242       str2.Remove(0, 1);
243       
244       returnCode = str2.Atoll();
245     }
246     else
247       AliFatal(Form("Invalid trigger syntax: %s", trigger));
248   }
249   
250   delete tokens;
251   
252   if (foundBCRequirement && !foundCorrectBC)
253     return kFALSE;
254   
255   return returnCode;
256 }
257     
258 UInt_t AliPhysicsSelection::IsCollisionCandidate(const AliESDEvent* aEsd)
259 {
260   // checks if the given event is a collision candidate
261   //
262   // returns a bit word describing the fired offline triggers (see AliVEvent::EOfflineTriggerTypes)
263   
264   if (fCurrentRun != aEsd->GetRunNumber()) {
265     if (!Initialize(aEsd->GetRunNumber()))
266       AliFatal(Form("Could not initialize for run %d", aEsd->GetRunNumber()));
267     if(fComputeBG) SetBIFactors(aEsd); // is this safe here?
268   }
269   const AliESDHeader* esdHeader = aEsd->GetHeader();
270   if (!esdHeader)
271   {
272     AliError("ESD Header could not be retrieved");
273     return kFALSE;
274   }
275   
276   // check event type; should be PHYSICS = 7 for data and 0 for MC
277   if (!fMC)
278   {
279     if (esdHeader->GetEventType() != 7)
280       return kFALSE;
281   }
282   else
283   {
284     if (esdHeader->GetEventType() != 0)
285       AliFatal(Form("Invalid event type for MC: %d", esdHeader->GetEventType()));
286   }
287   
288   UInt_t accept = 0;
289     
290   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
291   for (Int_t i=0; i < count; i++)
292   {
293     const char* triggerClass = 0;
294     if (i < fCollTrigClasses.GetEntries())
295       triggerClass = ((TObjString*) fCollTrigClasses.At(i))->String();
296     else
297       triggerClass = ((TObjString*) fBGTrigClasses.At(i - fCollTrigClasses.GetEntries()))->String();
298   
299     AliDebug(AliLog::kDebug, Form("Processing trigger class %s", triggerClass));
300   
301     AliTriggerAnalysis* triggerAnalysis = static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i));
302   
303     triggerAnalysis->FillTriggerClasses(aEsd);
304     
305     UInt_t singleTriggerResult = CheckTriggerClass(aEsd, triggerClass);
306     if (singleTriggerResult)
307     {
308       triggerAnalysis->FillHistograms(aEsd);
309   
310       Bool_t isBin0 = kFALSE;
311       if (fBin0CallBack != "") {
312         AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
313         if (!mgr) {
314           AliError("Cannot get the analysis manager");
315         }  
316         else {
317           isBin0 = ((AliAnalysisTaskSE*)mgr->GetTask(fBin0CallBack.Data()))->IsEventInBinZero();
318         }
319       } else if (fBin0CallBackPointer) {
320           isBin0 = (*fBin0CallBackPointer)(aEsd);
321         
322       }
323       
324
325       
326       // hardware trigger (should only remove events for MC)
327       // replay CINT1B hardware trigger
328       // TODO this has to depend on the actual hardware trigger (and that depends on the run...)
329       Int_t fastORHW = triggerAnalysis->SPDFiredChips(aEsd, 1); // SPD number of chips from trigger bits (!)
330       Bool_t v0A       = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0A);
331       Bool_t v0C       = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0C);
332       Bool_t v0AHW     = fSkipV0 ? 0 :(triggerAnalysis->V0Trigger(aEsd, AliTriggerAnalysis::kASide, kTRUE) == AliTriggerAnalysis::kV0BB);// should replay hw trigger
333       Bool_t v0CHW     = fSkipV0 ? 0 :(triggerAnalysis->V0Trigger(aEsd, AliTriggerAnalysis::kCSide, kTRUE) == AliTriggerAnalysis::kV0BB);// should replay hw trigger
334       // offline trigger
335       Int_t fastOROffline = triggerAnalysis->SPDFiredChips(aEsd, 0); // SPD number of chips from clusters (!)
336       Bool_t v0ABG = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0ABG);
337       Bool_t v0CBG = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0CBG);
338       Bool_t v0BG = v0ABG || v0CBG;
339
340       // fmd
341       Bool_t fmdA =  triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kFMDA);
342       Bool_t fmdC =  triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kFMDC);
343       Bool_t fmd  = fmdA || fmdC;
344     
345       // SSD
346       Int_t ssdClusters = triggerAnalysis->SSDClusters(aEsd);
347
348       // Some "macros"
349       Bool_t mb1 = (fastOROffline > 0 || v0A || v0C) && (!v0BG);
350       Bool_t mb1prime = (fastOROffline > 1 || (fastOROffline > 0 && (v0A || v0C)) || (v0A && v0C) ) && (!v0BG);
351
352       // Background rejection
353       Bool_t bgID = kFALSE;
354       if (fBackgroundIdentification)
355         bgID = ! fBackgroundIdentification->IsSelected(const_cast<AliESDEvent*> (aEsd));
356       
357       Int_t ntrig = fastOROffline; // any 2 hits
358       if(v0A)              ntrig += 1;
359       if(v0C)              ntrig += 1; //v0C alone is enough
360       if(fmd)              ntrig += 1;
361       if(ssdClusters>1)    ntrig += 1;
362
363
364       Bool_t hwTrig = fastORHW > 0 || v0AHW || v0CHW;
365
366       // Fill trigger pattern histo
367       Int_t tpatt = 0;
368       if (fastORHW>0) tpatt+=1;
369       if (v0AHW)      tpatt+=2;
370       if (v0CHW)      tpatt+=4;
371       fHistTriggerPattern->Fill( tpatt );
372
373       // fill statistics and return decision
374       const Int_t nHistStat = 2;
375       for(Int_t iHistStat = 0; iHistStat < nHistStat; iHistStat++){
376         if (iHistStat == kStatIdxBin0 && !isBin0) continue; // skip the filling of bin0 stats if the event is not in the bin0
377       
378         fHistStatistics[iHistStat]->Fill(kStatTriggerClass, i);
379
380
381
382         // We fill the rest only if hw trigger is ok
383         if (!hwTrig)
384           {
385             AliDebug(AliLog::kDebug, "Rejecting event because hardware trigger is not fired");
386             continue;
387           } else {       
388           fHistStatistics[iHistStat]->Fill(kStatHWTrig, i);
389         }
390       
391
392         // v0 BG stats
393         if (v0ABG)
394           fHistStatistics[iHistStat]->Fill(kStatV0ABG, i);
395         if (v0CBG)
396           fHistStatistics[iHistStat]->Fill(kStatV0CBG, i);
397
398         // We fill the rest only if mb1 && ! v0BG
399         if (mb1)
400           fHistStatistics[iHistStat]->Fill(kStatMB1, i);
401         else continue;
402
403         if (mb1prime)
404           fHistStatistics[iHistStat]->Fill(kStatMB1Prime, i);
405
406         if (fmd)
407           fHistStatistics[iHistStat]->Fill(kStatFMD, i);
408
409         if(ssdClusters>1)
410           fHistStatistics[iHistStat]->Fill(kStatSSD1, i);
411
412         if(ntrig >= 2 && !v0BG) 
413           fHistStatistics[iHistStat]->Fill(kStatAny2Hits, i);
414
415         if (fastOROffline > 0)
416           fHistStatistics[iHistStat]->Fill(kStatFO1, i);
417         if (fastOROffline > 1)
418           fHistStatistics[iHistStat]->Fill(kStatFO2, i);
419         
420         if (v0A)
421           fHistStatistics[iHistStat]->Fill(kStatV0A, i);
422         if (v0C)
423           fHistStatistics[iHistStat]->Fill(kStatV0C, i);
424         
425         //       if (fastOROffline > 1 && !v0BG)
426         //         fHistStatistics[iHistStat]->Fill(kStatFO2NoBG, i);
427             
428         if (fastOROffline > 0 && (v0A || v0C) && !v0BG)
429           fHistStatistics[iHistStat]->Fill(kStatFO1AndV0, i);
430   
431         if (v0A && v0C && !v0BG && !bgID)
432           fHistStatistics[iHistStat]->Fill(kStatV0, i);
433
434
435         if ( mb1 )
436         
437           {
438             if (!v0BG || fSkipV0)
439               {
440                 if (!v0BG) fHistStatistics[iHistStat]->Fill(kStatOffline, i);
441       
442                 if (fBackgroundIdentification && bgID)
443                   {
444                     AliDebug(AliLog::kDebug, "Rejecting event because of background identification");
445                     fHistStatistics[iHistStat]->Fill(kStatBG, i);
446                   }
447                 else
448                   {
449                     AliDebug(AliLog::kDebug, "Accepted event for histograms");
450             
451                     fHistStatistics[iHistStat]->Fill(kStatAccepted, i);
452                     if(iHistStat == kStatIdxAll) fHistBunchCrossing->Fill(aEsd->GetBunchCrossNumber(), i); // Fill only for all (avoid double counting)
453                     if((i < fCollTrigClasses.GetEntries() || fSkipTriggerClassSelection) && (iHistStat==kStatIdxAll))
454                       accept |= singleTriggerResult; // only set for "all" (should not really matter)
455                   }
456               }
457             else
458               AliDebug(AliLog::kDebug, "Rejecting event because of V0 BG flag");
459           }
460         else
461           AliDebug(AliLog::kDebug, "Rejecting event because trigger condition is not fulfilled");
462       }
463     }
464   }
465  
466   if (accept)
467     AliDebug(AliLog::kDebug, Form("Accepted event as collision candidate with bit mask %d", accept));
468   
469   return accept;
470 }
471
472 Int_t AliPhysicsSelection::GetTriggerScheme(UInt_t runNumber) const
473 {
474   // returns the current trigger scheme (classes that are accepted/rejected)
475   
476   if (fMC)
477     return 0;
478     
479   // TODO dependent on run number
480   
481   switch (runNumber)
482   {
483     // CSMBB triggers
484     case 104044:
485     case 105054:
486     case 105057:
487       return 2;
488   }
489   
490   // default: CINT1 suite
491   return 1;
492 }  
493
494 const char * AliPhysicsSelection::GetFillingScheme(UInt_t runNumber)  {
495
496   if(fMC) return "MC";
497
498   if      (runNumber >= 104065 && runNumber <= 104160) {
499     return "4x4a";
500   } 
501   else if (runNumber >= 104315 && runNumber <= 104321) {
502     return "4x4a*";
503   }
504   else if (runNumber >= 104792 && runNumber <= 104803) {
505     return "4x4b";
506   }
507   else if (runNumber >= 104824 && runNumber <= 104892) {
508     return "4x4c";
509   }
510   else if (runNumber == 105143 || runNumber == 105160) {
511     return "16x16a";
512   }
513   else if (runNumber >= 105256 && runNumber <= 105268) {
514     return "4x4c";
515   } 
516   else if (runNumber >= 114786 && runNumber <= 116684) {
517     return "Single_2b_1_1_1";
518   }
519   else if (runNumber >= 117048 && runNumber <= 117120) {
520     return "Single_3b_2_2_2";
521   }
522   else if (runNumber >= 117220 && runNumber <= 119163) {
523     return "Single_2b_1_1_1";
524   }
525   else if (runNumber >= 119837 && runNumber <= 119862) {
526     return "Single_4b_2_2_2";
527   }
528   else if (runNumber >= 119902 && runNumber <= 120691) {
529     return "Single_6b_3_3_3";
530   }
531   else if (runNumber >= 120741 && runNumber <= 122375) {
532     return "Single_13b_8_8_8";
533   }
534   else if (runNumber >= 130148 && runNumber <= 130375) {
535     return "125n_48b_36_16_36";
536   } 
537   else {
538     AliError(Form("Unknown filling scheme (run %d)", runNumber));
539   }
540
541   return "Unknown";
542 }
543
544 const char * AliPhysicsSelection::GetBXIDs(UInt_t runNumber, const char * trigger)  {
545
546   if (!fUseBXNumbers || fMC) return "";
547
548   if      (runNumber >= 104065 && runNumber <= 104160) {
549     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #2128 #3019";
550     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #346 #3465";
551     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1234 #1680";
552     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
553     else AliError(Form("Unknown trigger: %s", trigger));
554   } 
555   else if (runNumber >= 104315 && runNumber <= 104321) {
556     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #2000 #2891";
557     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #218 #3337";
558     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1106 #1552";
559     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
560     else AliError(Form("Unknown trigger: %s", trigger));
561   }
562   else if (runNumber >= 104792 && runNumber <= 104803) {
563     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #2228 #3119";
564     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2554 #446";
565     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1334 #769";
566     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
567     else AliError(Form("Unknown trigger: %s", trigger));
568   }
569   else if (runNumber >= 104824 && runNumber <= 104892) {
570     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #3119 #769";
571     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2554 #446";
572     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1334 #2228";
573     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
574     else AliError(Form("Unknown trigger: %s", trigger));
575   }
576   else if (runNumber == 105143 || runNumber == 105160) {
577     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #1337 #1418 #2228 #2309 #3119 #3200 #446 #527";
578     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #1580  #1742  #1904  #2066  #2630  #2792  #2954  #3362";
579     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "  #845  #1007  #1169   #1577 #3359 #3521 #119  #281 ";
580     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
581     else AliError(Form("Unknown trigger: %s", trigger));
582   }
583   else if (runNumber >= 105256 && runNumber <= 105268) {
584     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #3019 #669";
585     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2454 #346";
586     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1234 #2128";
587     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
588     else AliError(Form("Unknown trigger: %s", trigger));
589   } else if (runNumber >= 114786 && runNumber <= 116684) { // 7 TeV 2010, assume always the same filling scheme
590     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #346";
591     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2131";
592     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #3019";
593     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #1238";
594     else AliError(Form("Unknown trigger: %s", trigger));
595   }
596   else if (runNumber >= 117048 && runNumber <= 117120) {
597     //    return "Single_3b_2_2_2";
598    if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346  #1240 ";
599    else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #2131 ";
600    else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019 ";
601    else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1238";
602    else AliError(Form("Unknown trigger: %s", trigger));
603
604   }
605   else if (runNumber >= 117220 && runNumber <= 119163) {
606     //    return "Single_2b_1_1_1";
607     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346 ";
608     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #2131 ";
609     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019 ";
610     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1238 ";
611     else AliError(Form("Unknown trigger: %s", trigger));                                                    
612   }
613   else if (runNumber >= 119837 && runNumber <= 119862) {
614     //    return "Single_4b_2_2_2";
615     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #669  #3019 ";
616     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #346  #2454 ";
617     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #1234  #2128 ";
618     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1681 #3463";
619     else AliError(Form("Unknown trigger: %s", trigger));
620
621   }
622   else if (runNumber >= 119902 && runNumber <= 120691) {
623     //    return "Single_6b_3_3_3";
624     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346  #546  #746 ";
625     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #2131  #2331  #2531 ";
626     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019  #3219  #3419 ";
627     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1296 #1670";
628     else AliError(Form("Unknown trigger: %s", trigger));
629   }
630   else if (runNumber >= 120741 && runNumber <= 122375) {
631     //    return "Single_13b_8_8_8";
632     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346  #446  #546  #646  #1240  #1340  #1440  #1540 ";
633     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #946  #2131  #2231  #2331  #2431 ";
634     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019  #3119  #3219  #3319  #3519 ";
635     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1835 #2726";
636     else AliError(Form("Unknown trigger: %s", trigger));
637     
638   } 
639   else if (runNumber >= 130148 && runNumber <= 130375) {
640     TString triggerString = trigger;
641     static TString returnString = " ";
642     returnString = "";
643     if (triggerString.Contains("B")) returnString += "   #346  #396  #446  #496  #546  #596  #646  #696  #1240  #1290  #1340  #1390  #1440  #1490  #1540  #1590 ";
644     if (triggerString.Contains("A")) returnString += "   #755  #805  #855  #905  #955  #1005  #1799  #1849  #1899  #2131  #2181  #2231  #2281  #2331  #2381  #2431  #2481  #2531  #2581  #2631  #2846  #3016  #3066  #3116  #3166  #3216  #3266  #3316  #3366  #3425  #3475  #3525 ";
645     if (triggerString.Contains("C")) returnString += "   #3019  #3069  #3119  #3169  #3219  #3269  #3319  #3369  #14  #64  #114  #746  #796  #846  #908  #958  #1008  #1640  #1690  #1740  #2055  #2125  #2175  #2225  #2275  #2325  #2375  #2425  #2475  #2534  #2584  #2634 ";
646     // Printf("0x%x",returnString.Data());
647     // Printf("%s",returnString.Data());
648     return returnString.Data();
649   }
650
651   else {
652     AliError(Form("Unknown run %d, using all BXs!",runNumber));
653   }
654
655   return "";
656 }
657     
658 Bool_t AliPhysicsSelection::Initialize(Int_t runNumber)
659 {
660   // initializes the object for the given run
661   // 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
662   
663   Bool_t oldStatus = TH1::AddDirectoryStatus();
664   TH1::AddDirectory(kFALSE);
665   
666   if(!fBin0CallBack) 
667     AliError("Bin0 Callback not set: will not fill the statistics for the bin 0");
668
669   if (fMC) {
670     // override BX and bg options in case of MC
671     fComputeBG    = kFALSE;
672     fUseBXNumbers = kFALSE;
673   }
674
675   Int_t triggerScheme = GetTriggerScheme(runNumber);
676   if (!fUsingCustomClasses && fCurrentRun != -1 && triggerScheme != GetTriggerScheme(fCurrentRun))
677     AliFatal("Processing several runs with different trigger schemes is not supported");
678   
679   if(fComputeBG && fCurrentRun != -1 && fCurrentRun != runNumber) 
680     AliFatal("Cannot process several runs because BG computation is requested");
681
682   if(fComputeBG && !fUseBXNumbers) 
683     AliFatal("Cannot compute BG id BX numbers are not used");
684   
685   if(fUseBXNumbers && fFillingScheme != "" && fFillingScheme != GetFillingScheme(runNumber))
686     AliFatal("Cannot process runs with different filling scheme if usage of BX numbers is requested");
687
688   fFillingScheme      = GetFillingScheme(runNumber);
689
690   AliInfo(Form("Initializing for run %d", runNumber));
691   
692   // initialize first time?
693   if (fCurrentRun == -1)
694   {
695     if (fUsingCustomClasses) {
696       AliInfo("Using user-provided trigger classes");
697     } else {
698       switch (triggerScheme)
699       {
700       case 0:
701         fCollTrigClasses.Add(new TObjString(Form("&%u", (UInt_t) AliVEvent::kMB)));
702         break;
703         
704       case 1:
705         // trigger classes used before August 2010
706         fCollTrigClasses.Add(new TObjString(Form("%s%s &%u","+CINT1B-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1B-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
707         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CINT1A-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1A-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
708         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CINT1C-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1C-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
709         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CINT1-E-NOPF-ALL",      GetBXIDs(runNumber,"CINT1-E-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
710
711         // Muon trigger have the same BXIDs of the corresponding CINT triggers
712         fCollTrigClasses.Add(new TObjString(Form("%s%s &%u","+CMUS1B-ABCE-NOPF-MUON",  GetBXIDs(runNumber,"CINT1B-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
713         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CMUS1A-ABCE-NOPF-MUON",  GetBXIDs(runNumber,"CINT1A-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
714         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CMUS1C-ABCE-NOPF-MUON",  GetBXIDs(runNumber,"CINT1C-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
715         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CMUS1-E-NOPF-MUON"    ,  GetBXIDs(runNumber,"CINT1-E-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
716         
717         // triggers classes used from August 2010
718         // MB
719         fCollTrigClasses.Add(new TObjString(Form("%s%s &%u", "+CINT1-B-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "B"),  (UInt_t) AliVEvent::kMB)));
720         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CINT1-AC-NOPF-ALLNOTRD", GetBXIDs(runNumber, "AC"), (UInt_t) AliVEvent::kMB)));
721         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CINT1-E-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "E"),  (UInt_t) AliVEvent::kMB)));
722                                                                                       
723         // MUON                                                                       
724         fCollTrigClasses.Add(new TObjString(Form("%s%s &%u", "+CMUS1-B-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "B"),  (UInt_t) AliVEvent::kMUON)));
725         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CMUS1-AC-NOPF-ALLNOTRD", GetBXIDs(runNumber, "AC"), (UInt_t) AliVEvent::kMUON)));
726         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CMUS1-E-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "E"),  (UInt_t) AliVEvent::kMUON)));
727                                                                                      
728         // High Multiplicity                                                         
729         fCollTrigClasses.Add(new TObjString(Form("%s%s &%u", "+CSH1-B-NOPF-ALLNOTRD"  , GetBXIDs(runNumber, "B"),  (UInt_t) AliVEvent::kHighMult)));
730         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CSH1-AC-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "AC"), (UInt_t) AliVEvent::kHighMult)));
731         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CSH1-E-NOPF-ALLNOTRD"  , GetBXIDs(runNumber, "E"),  (UInt_t) AliVEvent::kHighMult)));
732
733         // WARNING: IF YOU ADD MORE TRIGGER CLASSES, PLEASE CHECK THAT THE REGULAR EXPRESSION IN GetStatRow IS STILL VALID
734
735         break;
736         
737       case 2:
738         fCollTrigClasses.Add(new TObjString(Form("+CSMBB-ABCE-NOPF-ALL &%u", (UInt_t) AliVEvent::kMB)));
739         fBGTrigClasses.Add(new TObjString(Form("+CSMBA-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL &%u", (UInt_t) AliVEvent::kMB)));
740         fBGTrigClasses.Add(new TObjString(Form("+CSMBC-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL &%u", (UInt_t) AliVEvent::kMB)));
741         break;
742         
743       case 3:
744         // 
745         break;
746         
747       default:
748         AliFatal(Form("Unsupported trigger scheme %d", triggerScheme));
749       }
750     }
751     
752     Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
753     
754     for (Int_t i=0; i<count; i++)
755     {
756       AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis;
757       triggerAnalysis->SetAnalyzeMC(fMC);
758       triggerAnalysis->EnableHistograms();
759       triggerAnalysis->SetSPDGFOThreshhold(1);
760       fTriggerAnalysis.Add(triggerAnalysis);
761     }
762       
763     // TODO: shall I really delete this?
764     if (fHistStatistics[0])
765       delete fHistStatistics[0];
766     if (fHistStatistics[1])
767       delete fHistStatistics[1];
768   
769     fHistStatistics[kStatIdxBin0] = BookHistStatistics("_Bin0");
770     fHistStatistics[kStatIdxAll]  = BookHistStatistics("");
771     
772     if (fHistBunchCrossing)
773       delete fHistBunchCrossing;
774   
775     fHistBunchCrossing = new TH2F("fHistBunchCrossing", "fHistBunchCrossing;bunch crossing number;", 4000, -0.5, 3999.5,  count, -0.5, -0.5 + count);
776
777     if (fHistTriggerPattern)
778       delete fHistTriggerPattern;
779     
780     const int ntrig=3;
781     Int_t n = 1;
782     const Int_t nbinTrig = TMath::Nint(TMath::Power(2,ntrig));
783
784     fHistTriggerPattern = new TH1F("fHistTriggerPattern", "Trigger pattern: FO + 2*v0A + 4*v0C", 
785                                    nbinTrig, -0.5, nbinTrig-0.5);    
786     fHistTriggerPattern->GetXaxis()->SetBinLabel(1,"NO TRIG");
787     fHistTriggerPattern->GetXaxis()->SetBinLabel(2,"FO");
788     fHistTriggerPattern->GetXaxis()->SetBinLabel(3,"v0A");
789     fHistTriggerPattern->GetXaxis()->SetBinLabel(4,"FO & v0A");
790     fHistTriggerPattern->GetXaxis()->SetBinLabel(5,"v0C");
791     fHistTriggerPattern->GetXaxis()->SetBinLabel(6,"FO & v0C");
792     fHistTriggerPattern->GetXaxis()->SetBinLabel(7,"v0A & v0C");
793     fHistTriggerPattern->GetXaxis()->SetBinLabel(8,"FO & v0A & v0C");
794
795   
796     n = 1;
797     for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
798     {
799       fHistBunchCrossing->GetYaxis()->SetBinLabel(n, ((TObjString*) fCollTrigClasses.At(i))->String());
800       n++;
801     }
802     for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
803     {
804       fHistBunchCrossing->GetYaxis()->SetBinLabel(n, ((TObjString*) fBGTrigClasses.At(i))->String());
805       n++;
806     }
807
808     
809
810   }
811     
812   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
813   for (Int_t i=0; i<count; i++)
814   {
815     AliTriggerAnalysis* triggerAnalysis = static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i));
816   
817     switch (runNumber)
818     {
819       case 104315:
820       case 104316:
821       case 104320:
822       case 104321:
823       case 104439:
824         triggerAnalysis->SetV0TimeOffset(7.5);
825         break;
826       default:
827         triggerAnalysis->SetV0TimeOffset(0);
828     }
829   }
830     
831   fCurrentRun = runNumber;
832   
833   TH1::AddDirectory(oldStatus);
834   
835   return kTRUE;
836 }
837
838 TH2F * AliPhysicsSelection::BookHistStatistics(const char * tag) {
839   // add 6 rows to count for the estimate of good, accidentals and
840   // BG and the ratio of BG and accidentals to total +ratio goot to
841   // first col + 2 for error on good.
842   // TODO: Remember the the indexes of rows for the BG selection. Add new member fBGRows[] and use kStat as indexes
843
844   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
845 #ifdef VERBOSE_STAT
846   Int_t extrarows = fComputeBG != 0 ? 8 : 0;
847 #else
848   Int_t extrarows = fComputeBG != 0 ? 3 : 0;
849 #endif
850   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);
851
852   h->GetXaxis()->SetBinLabel(kStatTriggerClass,  "Trigger class");
853   h->GetXaxis()->SetBinLabel(kStatHWTrig,        "Hardware trigger");
854   h->GetXaxis()->SetBinLabel(kStatFO1,           "FO >= 1");
855   h->GetXaxis()->SetBinLabel(kStatFO2,           "FO >= 2");
856   h->GetXaxis()->SetBinLabel(kStatV0A,           "V0A");
857   h->GetXaxis()->SetBinLabel(kStatV0C,           "V0C");
858   h->GetXaxis()->SetBinLabel(kStatFMD,           "FMD");
859   h->GetXaxis()->SetBinLabel(kStatSSD1,          "SSD >= 2");
860   h->GetXaxis()->SetBinLabel(kStatV0ABG,         "V0A BG");
861   h->GetXaxis()->SetBinLabel(kStatV0CBG,         "V0C BG");
862   h->GetXaxis()->SetBinLabel(kStatMB1,           "(FO >= 1 | V0A | V0C) & !V0 BG");
863   h->GetXaxis()->SetBinLabel(kStatMB1Prime,      "(FO >= 2 | (FO >= 1 & (V0A | V0C)) | (V0A &v0C) ) & !V0 BG");
864   h->GetXaxis()->SetBinLabel(kStatFO1AndV0,      "FO >= 1 & (V0A | V0C) & !V0 BG");
865   h->GetXaxis()->SetBinLabel(kStatV0,            "V0A & V0C & !V0 BG & !BG ID");
866   h->GetXaxis()->SetBinLabel(kStatOffline,       "Offline Trigger");
867   h->GetXaxis()->SetBinLabel(kStatAny2Hits,      "2 Hits & !V0 BG");
868   h->GetXaxis()->SetBinLabel(kStatBG,            "Background identification");
869   h->GetXaxis()->SetBinLabel(kStatAccepted,      "Accepted");
870
871   Int_t n = 1;
872   for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
873     {
874       h->GetYaxis()->SetBinLabel(n, ((TObjString*) fCollTrigClasses.At(i))->String());
875       n++;
876     }
877   for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
878     {
879       h->GetYaxis()->SetBinLabel(n, ((TObjString*) fBGTrigClasses.At(i))->String());
880       n++;
881     }
882
883   if(fComputeBG) {
884     fBGStatOffset = n;
885     h->GetYaxis()->SetBinLabel(n++, Form("BG (A+C) (Mask [0x%x])", fComputeBG));
886     h->GetYaxis()->SetBinLabel(n++, "ACC");
887 #ifdef VERBOSE_STAT
888     h->GetYaxis()->SetBinLabel(n++, "BG (A+C) %  (rel. to CINT1B)");
889     h->GetYaxis()->SetBinLabel(n++, "ACC % (rel. to CINT1B)");
890     h->GetYaxis()->SetBinLabel(n++, "ERR GOOD %");
891     h->GetYaxis()->SetBinLabel(n++, "GOOD % (rel. to 1st col)");
892     h->GetYaxis()->SetBinLabel(n++, "ERR GOOD");
893 #endif
894     h->GetYaxis()->SetBinLabel(n++, "GOOD");
895   }
896
897   return h;
898 }
899
900 void AliPhysicsSelection::Print(Option_t* /* option */) const
901 {
902   // print the configuration
903   
904   Printf("Configuration initialized for run %d (MC: %d):", fCurrentRun, fMC);
905   
906   Printf("Collision trigger classes:");
907   for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
908     Printf("%s", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
909   
910   Printf("Background trigger classes:");
911   for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
912     Printf("%s", ((TObjString*) fBGTrigClasses.At(i))->String().Data());
913
914   AliTriggerAnalysis* triggerAnalysis = dynamic_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(0));
915   
916   if (triggerAnalysis)
917   {
918     if (triggerAnalysis->GetV0TimeOffset() > 0)
919       Printf("V0 time offset active: %.2f ns", triggerAnalysis->GetV0TimeOffset());
920     
921     Printf("\nTotal available events:");
922     
923     triggerAnalysis->PrintTriggerClasses();
924   }
925   
926   if (fHistStatistics[kStatIdxAll])
927   {
928     for (Int_t i=0; i<fCollTrigClasses.GetEntries(); i++)
929     {
930       Printf("\nSelection statistics for collision trigger %s:", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
931       
932       Printf("Total events with correct trigger class: %d", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(1, i+1));
933       Printf("Selected collision candidates: %d", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(fHistStatistics[kStatIdxAll]->GetXaxis()->FindBin("Accepted"), i+1));
934     }
935   }
936   
937   if (fHistBunchCrossing)
938   {
939     Printf("\nBunch crossing statistics:");
940     
941     for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
942     {
943       TString str;
944       str.Form("Trigger %s has accepted events in the bunch crossings: ", fHistBunchCrossing->GetYaxis()->GetBinLabel(i));
945       
946       for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
947         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
948           str += Form("%d, ", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
949        
950       Printf("%s", str.Data());
951     }
952     
953     for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
954     {
955       Int_t countColl = 0;
956       Int_t countBG = 0;
957       for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
958       {
959         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
960         {
961           if (fCollTrigClasses.FindObject(fHistBunchCrossing->GetYaxis()->GetBinLabel(i)))
962             countColl++;
963           if (fBGTrigClasses.FindObject(fHistBunchCrossing->GetYaxis()->GetBinLabel(i)))
964             countBG++;
965         }
966       }
967       if (countColl > 0 && countBG > 0)
968         Printf("WARNING: Bunch crossing %d has collision and BG trigger classes active. Check BPTX functioning for this run!", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
969     }
970   }
971
972   if (fUsingCustomClasses)        
973     Printf("WARNING: Using custom trigger classes!");
974   if (fSkipTriggerClassSelection) 
975     Printf("WARNING: Skipping trigger class selection!");
976   if (fSkipV0) 
977     Printf("WARNING: Ignoring V0 information in selection");
978   if(!fBin0CallBack) 
979     Printf("WARNING: Callback not set: will not fill the statistics for the bin 0");
980
981 }
982
983 Long64_t AliPhysicsSelection::Merge(TCollection* list)
984 {
985   // Merge a list of AliMultiplicityCorrection objects with this (needed for
986   // PROOF).
987   // Returns the number of merged objects (including this).
988
989   if (!list)
990     return 0;
991
992   if (list->IsEmpty())
993     return 1;
994
995   TIterator* iter = list->MakeIterator();
996   TObject* obj;
997   
998   // collections of all histograms
999   const Int_t nHists = 9;
1000   TList collections[nHists];
1001
1002   Int_t count = 0;
1003   while ((obj = iter->Next())) {
1004
1005     AliPhysicsSelection* entry = dynamic_cast<AliPhysicsSelection*> (obj);
1006     if (entry == 0) 
1007       continue;
1008     // Update run number. If this one is not initialized (-1) take the one from 
1009     // the next physics selection to be merged with. In case of 2 different run
1010     // numbers issue a warning (should physics selections from different runs be 
1011     // merged together) A.G.
1012     Int_t currentRun = entry->GetCurrentRun();
1013     // Nothing to merge with since run number was not initialized.
1014     if (currentRun < 0) continue;
1015     if (fCurrentRun < 0) fCurrentRun = currentRun;
1016     if (fCurrentRun != currentRun)
1017        AliWarning(Form("Current run %d not matching the one to be merged with %d", fCurrentRun, currentRun));
1018     
1019     collections[0].Add(&(entry->fTriggerAnalysis));
1020     if (entry->fHistStatistics[0])
1021       collections[1].Add(entry->fHistStatistics[0]);
1022     if (entry->fHistStatistics[1])
1023       collections[2].Add(entry->fHistStatistics[1]);
1024     if (entry->fHistBunchCrossing)
1025       collections[3].Add(entry->fHistBunchCrossing);
1026     if (entry->fHistTriggerPattern)
1027       collections[4].Add(entry->fHistTriggerPattern);
1028     if (entry->fBackgroundIdentification)
1029       collections[5].Add(entry->fBackgroundIdentification);
1030
1031     count++;
1032   }
1033
1034   fTriggerAnalysis.Merge(&collections[0]);
1035   if (fHistStatistics[0])
1036     fHistStatistics[0]->Merge(&collections[1]);
1037   if (fHistStatistics[1])
1038     fHistStatistics[1]->Merge(&collections[2]);
1039   if (fHistBunchCrossing)
1040     fHistBunchCrossing->Merge(&collections[3]);
1041   if (fHistTriggerPattern)
1042     fHistTriggerPattern->Merge(&collections[4]);
1043   if (fBackgroundIdentification)
1044     fBackgroundIdentification->Merge(&collections[5]);
1045   
1046   delete iter;
1047
1048   return count+1;
1049 }
1050
1051 void AliPhysicsSelection::SaveHistograms(const char* folder) const
1052 {
1053   // write histograms to current directory
1054   
1055   if (!fHistStatistics[0] || !fHistStatistics[1])
1056     return;
1057     
1058   if (folder)
1059   {
1060     gDirectory->mkdir(folder);
1061     gDirectory->cd(folder);
1062   }
1063   
1064
1065   // Fill the last rows of fHistStatistics before saving
1066   if (fComputeBG) {
1067     Int_t triggerScheme = GetTriggerScheme(UInt_t(fCurrentRun));
1068     if(triggerScheme != 1){
1069       AliWarning("BG estimate only supported for trigger scheme \"1\" (CINT1 suite)");
1070     } else if (fUseMuonTriggers) {
1071       AliWarning("BG estimate with muon triggers to be implemented");
1072     } else {      
1073       AliInfo("BG estimate assumes that for a given run you only have A and C triggers separately or"
1074               " toghether as a AC class! Make sure this assumption holds in your case"); 
1075       
1076       // use an anum for the different trigger classes, to make loops easier to read
1077       enum {kClassB =0 , kClassA, kClassC, kClassAC, kClassE, kNClasses};
1078       const char * classFlags[] = {"B", "A", "C", "AC", "E"}; // labels
1079
1080       UInt_t * rows[kNClasses] = {0}; // Array of matching rows
1081       Int_t nrows[kNClasses] = {0};
1082       // Get rows matching the requested trigger bits for all trigger classes
1083       for(Int_t iTrigClass = 0; iTrigClass < kNClasses; iTrigClass++){
1084         nrows[iTrigClass] = GetStatRow(classFlags[iTrigClass],fComputeBG,&rows[iTrigClass]);    
1085       }
1086
1087       // 0. Determine the ratios of triggers E/B, A/B, C/B from the stat histogram
1088       // Those are used to rescale the different classes to the same number of bx ids
1089       // TODO: pass names of the rows for B, CA and E and look names of the rows. How do I handle the case in which both AC are in the same row?         
1090       Int_t nBXIds[kNClasses] = {0};
1091       cout <<"Computing BG:" << endl;
1092       
1093       for(Int_t iTrigClass = 0; iTrigClass < kNClasses; iTrigClass++){
1094         for(Int_t irow = 0; irow < nrows[iTrigClass]; irow++) {
1095           if(irow==0) cout << "- Class " << classFlags[iTrigClass] << endl;   
1096           for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++) {
1097             if (fHistBunchCrossing->GetBinContent(j, rows[iTrigClass][irow]) > 0) {
1098               nBXIds[iTrigClass]++;      
1099             }
1100           }
1101           if(nBXIds[iTrigClass]>0) cout << "   Using row " << rows[iTrigClass][irow] <<  ": " 
1102                                         << fHistBunchCrossing->GetYaxis()->GetBinLabel(rows[iTrigClass][irow]) 
1103                                         << " (nBXID "<< nBXIds[iTrigClass] << ")"<< endl;
1104
1105         }
1106
1107       }
1108
1109       Float_t ratioToB[kNClasses];
1110       ratioToB[kClassE]  = nBXIds[kClassE]  >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassE]   : 0;
1111       ratioToB[kClassA]  = nBXIds[kClassA]  >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassA]   : 0;
1112       ratioToB[kClassC]  = nBXIds[kClassC]  >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassC]   : 0;
1113       ratioToB[kClassAC] = nBXIds[kClassAC] >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassAC]  : 0;
1114       Printf("Ratio between the BX ids in the different trigger classes:");
1115       Printf("  B/E  = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassE], ratioToB[kClassE] );
1116       Printf("  B/A  = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassA], ratioToB[kClassA] );
1117       Printf("  B/C  = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassC], ratioToB[kClassC] );
1118       Printf("  B/AC = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassAC],ratioToB[kClassAC]);
1119       Int_t nHistStat = 2;
1120
1121       // 1. loop over all cols
1122       for(Int_t iHistStat = 0; iHistStat < nHistStat; iHistStat++){
1123         Int_t ncol = fHistStatistics[iHistStat]->GetNbinsX();
1124         Float_t good1 = 0;      
1125         for(Int_t icol = 1; icol <= ncol; icol++) {
1126           Int_t nEvents[kNClasses] = {0}; // number of events should be reset at every column
1127           // For all trigger classes, add up over row matching trigger mask (as selected before)
1128           for(Int_t iTrigClass = 0; iTrigClass < kNClasses; iTrigClass++){
1129             for(Int_t irow = 0; irow < nrows[iTrigClass]; irow++) {       
1130               nEvents[iTrigClass] += (Int_t) fHistStatistics[iHistStat]->GetBinContent(icol,rows[iTrigClass][irow]);    
1131             }
1132             //      cout << "Events " << classFlags[iTrigClass] << " ("<<icol<<") " << nEvents[iTrigClass] << endl;         
1133           }
1134           if (nEvents[kClassB]>0) {
1135             Int_t acc  = ratioToB[kClassE]*nEvents[kClassE]; 
1136             Double_t acc_err = TMath::Sqrt(ratioToB[kClassE]*ratioToB[kClassE]*nEvents[kClassE]);
1137             //      Int_t bg   = cint1A + cint1C - 2*acc;
1138             // cout << "-----------------------" << endl;
1139             // cout << "Factors: " << fBIFactorA << " " << fBIFactorC << " " << fBIFactorAC << endl;
1140             // cout << "Ratios: "  << ratioToB[kClassA] << " " << ratioToB[kClassC] << " " << ratioToB[kClassAC] << endl;
1141             // cout << "Evts:   "  << nEvents[kClassA] << " " << nEvents[kClassC] << " " << nEvents[kClassAC] << " " <<  nEvents[kClassB] << endl;
1142             // cout << "Acc: " << acc << endl;
1143             
1144              // Assuming that for a given class the triggers are either recorded as A+C or AC
1145             Float_t bg   = nEvents[kClassAC] > 0 ?
1146               fBIFactorAC*(ratioToB[kClassAC]*nEvents[kClassAC] - 2*acc):
1147               fBIFactorA* (ratioToB[kClassA]*nEvents[kClassA]-acc) + 
1148               fBIFactorC* (ratioToB[kClassC]*nEvents[kClassC]-acc) ;
1149
1150             // cout << "BG: " << bg << endl;
1151             // cout  << "  " <<   fBIFactorA* (ratioToB[kClassA]*nEvents[kClassA]-acc) <<endl;
1152             // cout  << "  " <<   fBIFactorC* (ratioToB[kClassC]*nEvents[kClassC]-acc) <<endl;
1153             // cout  << "  " <<   fBIFactorAC*(ratioToB[kClassAC]*nEvents[kClassAC] - 2*acc) << endl;
1154             // cout << "-----------------------" << endl;
1155
1156             Float_t good = Float_t(nEvents[kClassB]) - bg - acc;
1157             if (icol ==1) good1 = good;
1158             //      Float_t errGood     = TMath::Sqrt(2*(nEvents[kClassA]+nEvents[kClassC]+nEvents[kClassE]));// Error on the number of goods assuming only bg fluctuates
1159             //      DeltaG^2 = B + FA^2 A + FC^2 C + Ratio^2 (FA+FC-1)^2 E.
1160             Float_t errGood     = nEvents[kClassAC] > 0 ? 
1161               TMath::Sqrt( nEvents[kClassB] +
1162                            fBIFactorAC*fBIFactorAC*ratioToB[kClassAC]*ratioToB[kClassAC]*nEvents[kClassAC] +
1163                            ratioToB[kClassE] * ratioToB[kClassE] * 
1164                            (fBIFactorAC - 1)*(fBIFactorAC - 1)*nEvents[kClassE]) :  
1165               TMath::Sqrt( nEvents[kClassB] + 
1166                            fBIFactorA*fBIFactorA*ratioToB[kClassA]*ratioToB[kClassA]*nEvents[kClassA] +
1167                            fBIFactorC*fBIFactorC*ratioToB[kClassC]*ratioToB[kClassC]*nEvents[kClassC] +
1168                            ratioToB[kClassE] * ratioToB[kClassE] * 
1169                            (fBIFactorA + fBIFactorC - 1)*(fBIFactorA + fBIFactorC - 1)*nEvents[kClassE]);
1170
1171             Float_t errBG = nEvents[kClassAC] > 0 ? 
1172               TMath::Sqrt(fBIFactorAC*fBIFactorAC*ratioToB[kClassAC]*ratioToB[kClassAC]*nEvents[kClassAC]+
1173                           4*ratioToB[kClassE]*ratioToB[kClassE]*(fBIFactorAC*fBIFactorAC)*nEvents[kClassE]) :
1174               TMath::Sqrt(fBIFactorA*fBIFactorA*ratioToB[kClassA]*ratioToB[kClassA]*nEvents[kClassA]+
1175                           fBIFactorC*fBIFactorC*ratioToB[kClassC]*ratioToB[kClassC]*nEvents[kClassC]+
1176                           ratioToB[kClassE]*ratioToB[kClassE]*(fBIFactorA+fBIFactorC)*(fBIFactorA+fBIFactorC)*nEvents[kClassE]);
1177         
1178         
1179             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowBG,bg);        
1180             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowBG,errBG);     
1181             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowAcc,acc);      
1182             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowAcc,acc_err);  
1183             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowGood,good);    
1184             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowGood,errGood);    
1185
1186 #ifdef VERBOSE_STAT
1187             //kStatRowBG=0,kStatRowAcc,kStatRowBGFrac,kStatRowAccFrac,kStatRowErrGoodFrac,kStatRowGoodFrac,kStatRowGood,kStatRowErrGood
1188             Float_t accFrac   = Float_t(acc) / nEvents[kClassB]  *100;
1189             Float_t errAccFrac= Float_t(acc_err) / nEvents[kClassB]  *100;
1190             Float_t bgFrac    = Float_t(bg)  / nEvents[kClassB]  *100;
1191             Float_t goodFrac  = Float_t(good)  / good1 *100;
1192             Float_t errGoodFrac = errGood/good1 * 100;
1193             Float_t errFracBG = bg > 0 ? TMath::Sqrt((errBG/bg)*(errBG/bg) + 1/nEvents[kClassB])*bgFrac : 0;
1194             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowBGFrac,bgFrac);        
1195             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowBGFrac,errFracBG);     
1196             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowAccFrac,accFrac);    
1197             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowAccFrac,errAccFrac);    
1198             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowGoodFrac,goodFrac);    
1199             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowErrGoodFrac,errGoodFrac);    
1200             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowErrGood,errGood);    
1201 #endif
1202           }
1203         }
1204       }
1205     }  
1206   }
1207
1208   fHistStatistics[0]->Write();
1209   fHistStatistics[1]->Write();
1210   if(fHistBunchCrossing ) fHistBunchCrossing ->Write();
1211   if(fHistTriggerPattern) fHistTriggerPattern->Write();
1212   
1213   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
1214   for (Int_t i=0; i < count; i++)
1215   {
1216     TString triggerClass = "trigger_histograms_";
1217     if (i < fCollTrigClasses.GetEntries())
1218       triggerClass += ((TObjString*) fCollTrigClasses.At(i))->String();
1219     else
1220       triggerClass += ((TObjString*) fBGTrigClasses.At(i - fCollTrigClasses.GetEntries()))->String();
1221   
1222     gDirectory->mkdir(triggerClass);
1223     gDirectory->cd(triggerClass);
1224   
1225     static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i))->SaveHistograms();
1226     
1227     gDirectory->cd("..");
1228   }
1229  
1230   if (fBackgroundIdentification)
1231   {
1232     gDirectory->mkdir("background_identification");
1233     gDirectory->cd("background_identification");
1234       
1235     fBackgroundIdentification->GetOutput()->Write();
1236       
1237     gDirectory->cd("..");
1238   }
1239   
1240   if (folder)
1241     gDirectory->cd("..");
1242 }
1243
1244 Int_t AliPhysicsSelection::GetStatRow(const char * triggerBXClass, UInt_t offlineTriggerType, UInt_t ** rowIDs) const {
1245   // Puts inside the array rowIDs the row number for a given offline
1246   // trigger in a given bx class. Returns the total number of lines
1247   // matching the selection
1248   // triggerBXClass can be either "A", "AC", "B" or "E"
1249   // offlineTriggerType is one of the types defined in AliVEvent
1250   // User should delete rowIDs if no longer needed
1251
1252   if(!fHistStatistics[0]) {
1253     AliWarning("Not initialized, returning 0");
1254     return 0;
1255   }
1256   const Int_t nrows = fHistStatistics[0]->GetNbinsY();
1257
1258   // allocate memory for at maximum nrows
1259   Int_t nMatches = 0;
1260   (*rowIDs) = new UInt_t[nrows];
1261
1262   // Build regular expression. look for a +, followed by the beginning
1263   // of a word. Within the word, look for the class id after any
1264   // number of any char, but at most one dash ("[^-]*-?"), followed by
1265   // a - and then any char (".*") and at the class id ("\\&(\\d)")
1266   // The class id is stored.
1267   // WARNING: please check this if the trigger classes change
1268   TPRegexp re(Form("\\+\\b[^-]*-?%s-.*\\&(\\d)",triggerBXClass));  
1269   // Loop over rows and find matching ones:
1270   for(Int_t irow = 1; irow <= nrows; irow++){
1271     TObjArray * matches = re.MatchS(fHistStatistics[0]->GetYaxis()->GetBinLabel(irow));
1272     if (matches->GetEntries()) {
1273       TString s = ((TObjString*)matches->At(1))->GetString();      
1274       if(UInt_t(s.Atoi()) & offlineTriggerType) { // bitwise comparison with the requested mask
1275         //      cout << "Marching " << s.Data() << " " << offlineTriggerType << " " << fHistStatistics[0]->GetYaxis()->GetBinLabel(irow) << endl;       
1276         (*rowIDs)[nMatches] = irow;
1277         nMatches++;     
1278       }
1279     }
1280     delete matches;
1281   }
1282   
1283   return nMatches;
1284 }
1285
1286 void AliPhysicsSelection::SetBIFactors(const AliESDEvent * aESD) {
1287   // Set factors for realtive bunch intesities
1288   if(!aESD) { 
1289     AliFatal("ESD not given");
1290   }
1291   Int_t run = aESD->GetRunNumber();
1292   if (run > 105268) {
1293     // intensities stored in the ESDs
1294     const AliESDRun* esdRun = aESD->GetESDRun();
1295     Double_t intAB = esdRun->GetMeanIntensityIntecting(0);
1296     Double_t intCB = esdRun->GetMeanIntensityIntecting(1);
1297     Double_t intAA = esdRun->GetMeanIntensityNonIntecting(0);
1298     Double_t intCC = esdRun->GetMeanIntensityNonIntecting(1);
1299
1300     // cout << "INT " <<intAB <<endl;
1301     // cout << "INT " <<intCB <<endl;
1302     // cout << "INT " <<intAA <<endl;
1303     // cout << "INT " <<intCC <<endl;
1304
1305     if (intAB > -1 && intAA > -1) {
1306       fBIFactorA = intAB/intAA;
1307     } else {
1308       AliWarning("Cannot set fBIFactorA, assuming 1");
1309     }
1310     
1311     if (intCB > -1 && intCC > -1) {
1312       fBIFactorC = intCB/intCC;
1313     } else {
1314       AliWarning("Cannot set fBIFactorC, assuming 1");
1315     }
1316       
1317     if (intAB > -1 && intAA > -1 &&
1318         intCB > -1 && intCC > -1) {
1319       fBIFactorAC = (intAB+intCB)/(intAA+intCC);
1320     } else {
1321       AliWarning("Cannot set fBIFactorAC, assuming 1");
1322     }
1323         
1324   }  
1325   else {
1326     // First runs. Intensities hardcoded
1327     switch(run) {
1328     case 104155:
1329       fBIFactorA = 0.961912722908;
1330       fBIFactorC = 1.04992336081;
1331       break;
1332     case 104157:
1333       fBIFactorA = 0.947312854998;
1334       fBIFactorC = 1.01599706417;
1335       break;
1336     case 104159:
1337       fBIFactorA = 0.93659320151;
1338       fBIFactorC = 0.98580804207;
1339       break;
1340     case 104160:
1341       fBIFactorA = 0.929664189926;
1342       fBIFactorC = 0.963467679851;
1343       break;
1344     case 104315:
1345       fBIFactorA = 1.08939104979;
1346       fBIFactorC = 0.931113921925;
1347       break;
1348     case 104316:
1349       fBIFactorA = 1.08351880974;
1350       fBIFactorC = 0.916068345845;
1351       break;
1352     case 104320:
1353       fBIFactorA = 1.07669281245;
1354       fBIFactorC = 0.876818744763;
1355       break;
1356     case 104321:
1357       fBIFactorA = 1.00971079602;
1358       fBIFactorC = 0.773781299076;
1359       break;
1360     case 104792:
1361       fBIFactorA = 0.787215863962;
1362       fBIFactorC = 0.778253173071;
1363       break;
1364     case 104793:
1365       fBIFactorA = 0.692211363661;
1366       fBIFactorC = 0.733152456667;
1367       break;
1368     case 104799:
1369       fBIFactorA = 1.04027825161;
1370       fBIFactorC = 1.00530825942;
1371       break;
1372     case 104800:
1373       fBIFactorA = 1.05309910671;
1374       fBIFactorC = 1.00376801855;
1375       break;
1376     case 104801:
1377       fBIFactorA = 1.0531231922;
1378       fBIFactorC = 0.992439666758;
1379       break;
1380     case 104802:
1381       fBIFactorA = 1.04191478134;
1382       fBIFactorC = 0.979368585208;
1383       break;
1384     case 104803:
1385       fBIFactorA = 1.03121314094;
1386       fBIFactorC = 0.973379962609;
1387       break;
1388     case 104824:
1389       fBIFactorA = 0.969945926722;
1390       fBIFactorC = 0.39549745806;
1391       break;
1392     case 104825:
1393       fBIFactorA = 0.968627213937;
1394       fBIFactorC = 0.310100412205;
1395       break;
1396     case 104841:
1397       fBIFactorA = 0.991601393212;
1398       fBIFactorC = 0.83762204722;
1399       break;
1400     case 104845:
1401       fBIFactorA = 0.98040863886;
1402       fBIFactorC = 0.694824205793;
1403       break;
1404     case 104867:
1405       fBIFactorA = 1.10646173412;
1406       fBIFactorC = 0.841407246916;
1407       break;
1408     case 104876:
1409       fBIFactorA = 1.12063452421;
1410       fBIFactorC = 0.78726542895;
1411       break;
1412     case 104890:
1413       fBIFactorA = 1.02346137453;
1414       fBIFactorC = 1.03355663595;
1415       break;
1416     case 104892:
1417       fBIFactorA = 1.05406025913;
1418       fBIFactorC = 1.00029166135;
1419       break;
1420     case 105143:
1421       fBIFactorA = 0.947343384349;
1422       fBIFactorC = 0.972637444408;
1423       break;
1424     case 105160:
1425       fBIFactorA = 0.908854622177;
1426       fBIFactorC = 0.958851103977;
1427       break; 
1428     case 105256:
1429       fBIFactorA = 0.810076150206;
1430       fBIFactorC = 0.884663561883;
1431       break;
1432     case 105257:
1433       fBIFactorA = 0.80974912303;
1434       fBIFactorC = 0.878859123479;
1435       break;
1436     case 105268:
1437       fBIFactorA = 0.809052110679;
1438       fBIFactorC = 0.87233890989;
1439       break;
1440     default:
1441       fBIFactorA = 1;
1442       fBIFactorC = 1;
1443     }
1444   }
1445
1446 }