Fixed warning
[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 if (runNumber >= 130601 && runNumber <= 130640) {
538     return "1000ns_50b_35_14_35";
539   }
540   else {
541     AliError(Form("Unknown filling scheme (run %d)", runNumber));
542   }
543
544   return "Unknown";
545 }
546
547 const char * AliPhysicsSelection::GetBXIDs(UInt_t runNumber, const char * trigger)  {
548
549   if (!fUseBXNumbers || fMC) return "";
550
551   if      (runNumber >= 104065 && runNumber <= 104160) {
552     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #2128 #3019";
553     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #346 #3465";
554     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1234 #1680";
555     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
556     //    else AliError(Form("Unknown trigger: %s", trigger));
557   } 
558   else if (runNumber >= 104315 && runNumber <= 104321) {
559     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #2000 #2891";
560     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #218 #3337";
561     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1106 #1552";
562     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
563     //    else AliError(Form("Unknown trigger: %s", trigger));
564   }
565   else if (runNumber >= 104792 && runNumber <= 104803) {
566     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #2228 #3119";
567     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2554 #446";
568     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1334 #769";
569     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
570     //    else AliError(Form("Unknown trigger: %s", trigger));
571   }
572   else if (runNumber >= 104824 && runNumber <= 104892) {
573     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #3119 #769";
574     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2554 #446";
575     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1334 #2228";
576     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
577     //    else AliError(Form("Unknown trigger: %s", trigger));
578   }
579   else if (runNumber == 105143 || runNumber == 105160) {
580     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #1337 #1418 #2228 #2309 #3119 #3200 #446 #527";
581     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #1580  #1742  #1904  #2066  #2630  #2792  #2954  #3362";
582     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "  #845  #1007  #1169   #1577 #3359 #3521 #119  #281 ";
583     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
584     //    else AliError(Form("Unknown trigger: %s", trigger));
585   }
586   else if (runNumber >= 105256 && runNumber <= 105268) {
587     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #3019 #669";
588     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2454 #346";
589     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1234 #2128";
590     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
591     //    else AliError(Form("Unknown trigger: %s", trigger));
592   } else if (runNumber >= 114786 && runNumber <= 116684) { // 7 TeV 2010, assume always the same filling scheme
593     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #346";
594     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2131";
595     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #3019";
596     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #1238";
597     //    else AliError(Form("Unknown trigger: %s", trigger));
598   }
599   else if (runNumber >= 117048 && runNumber <= 117120) {
600     //    return "Single_3b_2_2_2";
601    if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346  #1240 ";
602    else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #2131 ";
603    else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019 ";
604    else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1238";
605    //   else AliError(Form("Unknown trigger: %s", trigger));
606
607   }
608   else if (runNumber >= 117220 && runNumber <= 119163) {
609     //    return "Single_2b_1_1_1";
610     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346 ";
611     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #2131 ";
612     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019 ";
613     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1238 ";
614     //    else AliError(Form("Unknown trigger: %s", trigger));                                              
615   }
616   else if (runNumber >= 119837 && runNumber <= 119862) {
617     //    return "Single_4b_2_2_2";
618     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #669  #3019 ";
619     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #346  #2454 ";
620     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #1234  #2128 ";
621     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1681 #3463";
622     //    else AliError(Form("Unknown trigger: %s", trigger));
623
624   }
625   else if (runNumber >= 119902 && runNumber <= 120691) {
626     //    return "Single_6b_3_3_3";
627     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346  #546  #746 ";
628     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #2131  #2331  #2531 ";
629     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019  #3219  #3419 ";
630     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1296 #1670";
631     //    else AliError(Form("Unknown trigger: %s", trigger));
632   }
633   else if (runNumber >= 120741 && runNumber <= 122375) {
634     //    return "Single_13b_8_8_8";
635     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346  #446  #546  #646  #1240  #1340  #1440  #1540 ";
636     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #946  #2131  #2231  #2331  #2431 ";
637     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019  #3119  #3219  #3319  #3519 ";
638     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1835 #2726";
639     //    else AliError(Form("Unknown trigger: %s", trigger));
640     
641   } 
642   else if (runNumber >= 130148 && runNumber <= 130375) {
643     TString triggerString = trigger;
644     static TString returnString = " ";
645     returnString = "";
646     if (triggerString.Contains("B")) returnString += "   #346  #396  #446  #496  #546  #596  #646  #696  #1240  #1290  #1340  #1390  #1440  #1490  #1540  #1590 ";
647     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 ";
648     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 ";
649     // Printf("0x%x",returnString.Data());
650     // Printf("%s",returnString.Data());
651     return returnString.Data();
652   } 
653   else if (runNumber >= 130601 && runNumber <= 130640) {
654     TString triggerString = trigger;
655     static TString returnString = " ";
656     returnString = "";
657     if (triggerString.Contains("B")) returnString += "  #346  #386  #426  #466  #506  #546  #586  #1240  #1280  #1320  #1360  #1400  #1440  #1480 ";
658     if (triggerString.Contains("A")) returnString += "  #626  #666  #706  #746  #786  #826  #866  #1520  #1560  #1600  #1640  #1680  #1720  #1760  #2076  #2131  #2171  #2211  #2251  #2291  #2331  #2371  #2414  #2454  #2494  #2534  #2574  #2614  #2654  #2694  #2734  #2774  #2814 "; //#2854  #2894  #2934 not present in this run
659     if (triggerString.Contains("C")) returnString += "  #3019  #3059  #3099  #3139  #3179  #3219  #3259  #3299  #3339  #3379  #3419  #3459  #3499  #3539  #115  #629  #669  #709  #749  #789  #829  #869  #909  #949  #989  #1029  #1069  #1109  #1149  #1523  #1563  #1603  #1643 "; //#1683  #1723  #1763 not present in this run
660     return returnString.Data();
661   }
662
663   else {
664     AliWarning(Form("Unknown run %d, using all BXs!",runNumber));
665   }
666
667   return "";
668 }
669     
670 Bool_t AliPhysicsSelection::Initialize(Int_t runNumber)
671 {
672   // initializes the object for the given run
673   // 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
674   
675   Bool_t oldStatus = TH1::AddDirectoryStatus();
676   TH1::AddDirectory(kFALSE);
677   
678   if(!fBin0CallBack) 
679     AliError("Bin0 Callback not set: will not fill the statistics for the bin 0");
680
681   if (fMC) {
682     // override BX and bg options in case of MC
683     fComputeBG    = kFALSE;
684     fUseBXNumbers = kFALSE;
685   }
686
687   Int_t triggerScheme = GetTriggerScheme(runNumber);
688   if (!fUsingCustomClasses && fCurrentRun != -1 && triggerScheme != GetTriggerScheme(fCurrentRun))
689     AliFatal("Processing several runs with different trigger schemes is not supported");
690   
691   if(fComputeBG && fCurrentRun != -1 && fCurrentRun != runNumber) 
692     AliFatal("Cannot process several runs because BG computation is requested");
693
694   if(fComputeBG && !fUseBXNumbers) 
695     AliFatal("Cannot compute BG if BX numbers are not used");
696   
697   if(fUseBXNumbers && fFillingScheme != "" && fFillingScheme != GetFillingScheme(runNumber))
698     AliFatal("Cannot process runs with different filling scheme if usage of BX numbers is requested");
699
700   fFillingScheme      = GetFillingScheme(runNumber);
701
702   AliInfo(Form("Initializing for run %d", runNumber));
703   
704   // initialize first time?
705   if (fCurrentRun == -1)
706   {
707     if (fUsingCustomClasses) {
708       AliInfo("Using user-provided trigger classes");
709     } else {
710       switch (triggerScheme)
711       {
712       case 0:
713         fCollTrigClasses.Add(new TObjString(Form("&%u", (UInt_t) AliVEvent::kMB)));
714         break;
715         
716       case 1:
717         // trigger classes used before August 2010
718         fCollTrigClasses.Add(new TObjString(Form("%s%s &%u","+CINT1B-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1B-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
719         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CINT1A-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1A-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
720         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CINT1C-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1C-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
721         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CINT1-E-NOPF-ALL",      GetBXIDs(runNumber,"CINT1-E-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
722
723         // Muon trigger have the same BXIDs of the corresponding CINT triggers
724         fCollTrigClasses.Add(new TObjString(Form("%s%s &%u","+CMUS1B-ABCE-NOPF-MUON",  GetBXIDs(runNumber,"CINT1B-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
725         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CMUS1A-ABCE-NOPF-MUON",  GetBXIDs(runNumber,"CINT1A-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
726         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CMUS1C-ABCE-NOPF-MUON",  GetBXIDs(runNumber,"CINT1C-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
727         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CMUS1-E-NOPF-MUON"    ,  GetBXIDs(runNumber,"CINT1-E-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
728         
729         // triggers classes used from August 2010
730         // MB
731         fCollTrigClasses.Add(new TObjString(Form("%s%s &%u", "+CINT1-B-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "B"),  (UInt_t) AliVEvent::kMB)));
732         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CINT1-AC-NOPF-ALLNOTRD", GetBXIDs(runNumber, "AC"), (UInt_t) AliVEvent::kMB)));
733         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CINT1-E-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "E"),  (UInt_t) AliVEvent::kMB)));
734                                                                                       
735         // MUON                                                                       
736         fCollTrigClasses.Add(new TObjString(Form("%s%s &%u", "+CMUS1-B-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "B"),  (UInt_t) AliVEvent::kMUON)));
737         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CMUS1-AC-NOPF-ALLNOTRD", GetBXIDs(runNumber, "AC"), (UInt_t) AliVEvent::kMUON)));
738         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CMUS1-E-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "E"),  (UInt_t) AliVEvent::kMUON)));
739                                                                                      
740         // High Multiplicity                                                         
741         fCollTrigClasses.Add(new TObjString(Form("%s%s &%u", "+CSH1-B-NOPF-ALLNOTRD"  , GetBXIDs(runNumber, "B"),  (UInt_t) AliVEvent::kHighMult)));
742         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CSH1-AC-NOPF-ALLNOTRD" , GetBXIDs(runNumber, "AC"), (UInt_t) AliVEvent::kHighMult)));
743         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u", "+CSH1-E-NOPF-ALLNOTRD"  , GetBXIDs(runNumber, "E"),  (UInt_t) AliVEvent::kHighMult)));
744
745         // WARNING: IF YOU ADD MORE TRIGGER CLASSES, PLEASE CHECK THAT THE REGULAR EXPRESSION IN GetStatRow IS STILL VALID
746
747         break;
748         
749       case 2:
750         fCollTrigClasses.Add(new TObjString(Form("+CSMBB-ABCE-NOPF-ALL &%u", (UInt_t) AliVEvent::kMB)));
751         fBGTrigClasses.Add(new TObjString(Form("+CSMBA-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL &%u", (UInt_t) AliVEvent::kMB)));
752         fBGTrigClasses.Add(new TObjString(Form("+CSMBC-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL &%u", (UInt_t) AliVEvent::kMB)));
753         break;
754         
755       case 3:
756         // 
757         break;
758         
759       default:
760         AliFatal(Form("Unsupported trigger scheme %d", triggerScheme));
761       }
762     }
763     
764     Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
765     
766     for (Int_t i=0; i<count; i++)
767     {
768       AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis;
769       triggerAnalysis->SetAnalyzeMC(fMC);
770       triggerAnalysis->EnableHistograms();
771       triggerAnalysis->SetSPDGFOThreshhold(1);
772       fTriggerAnalysis.Add(triggerAnalysis);
773     }
774       
775     // TODO: shall I really delete this?
776     if (fHistStatistics[0])
777       delete fHistStatistics[0];
778     if (fHistStatistics[1])
779       delete fHistStatistics[1];
780   
781     fHistStatistics[kStatIdxBin0] = BookHistStatistics("_Bin0");
782     fHistStatistics[kStatIdxAll]  = BookHistStatistics("");
783     
784     if (fHistBunchCrossing)
785       delete fHistBunchCrossing;
786   
787     fHistBunchCrossing = new TH2F("fHistBunchCrossing", "fHistBunchCrossing;bunch crossing number;", 4000, -0.5, 3999.5,  count, -0.5, -0.5 + count);
788
789     if (fHistTriggerPattern)
790       delete fHistTriggerPattern;
791     
792     const int ntrig=3;
793     Int_t n = 1;
794     const Int_t nbinTrig = TMath::Nint(TMath::Power(2,ntrig));
795
796     fHistTriggerPattern = new TH1F("fHistTriggerPattern", "Trigger pattern: FO + 2*v0A + 4*v0C", 
797                                    nbinTrig, -0.5, nbinTrig-0.5);    
798     fHistTriggerPattern->GetXaxis()->SetBinLabel(1,"NO TRIG");
799     fHistTriggerPattern->GetXaxis()->SetBinLabel(2,"FO");
800     fHistTriggerPattern->GetXaxis()->SetBinLabel(3,"v0A");
801     fHistTriggerPattern->GetXaxis()->SetBinLabel(4,"FO & v0A");
802     fHistTriggerPattern->GetXaxis()->SetBinLabel(5,"v0C");
803     fHistTriggerPattern->GetXaxis()->SetBinLabel(6,"FO & v0C");
804     fHistTriggerPattern->GetXaxis()->SetBinLabel(7,"v0A & v0C");
805     fHistTriggerPattern->GetXaxis()->SetBinLabel(8,"FO & v0A & v0C");
806
807   
808     n = 1;
809     for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
810     {
811       fHistBunchCrossing->GetYaxis()->SetBinLabel(n, ((TObjString*) fCollTrigClasses.At(i))->String());
812       n++;
813     }
814     for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
815     {
816       fHistBunchCrossing->GetYaxis()->SetBinLabel(n, ((TObjString*) fBGTrigClasses.At(i))->String());
817       n++;
818     }
819
820     
821
822   }
823     
824   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
825   for (Int_t i=0; i<count; i++)
826   {
827     AliTriggerAnalysis* triggerAnalysis = static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i));
828   
829     switch (runNumber)
830     {
831       case 104315:
832       case 104316:
833       case 104320:
834       case 104321:
835       case 104439:
836         triggerAnalysis->SetV0TimeOffset(7.5);
837         break;
838       default:
839         triggerAnalysis->SetV0TimeOffset(0);
840     }
841   }
842     
843   fCurrentRun = runNumber;
844   
845   TH1::AddDirectory(oldStatus);
846   
847   return kTRUE;
848 }
849
850 TH2F * AliPhysicsSelection::BookHistStatistics(const char * tag) {
851   // add 6 rows to count for the estimate of good, accidentals and
852   // BG and the ratio of BG and accidentals to total +ratio goot to
853   // first col + 2 for error on good.
854   // TODO: Remember the the indexes of rows for the BG selection. Add new member fBGRows[] and use kStat as indexes
855
856   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
857 #ifdef VERBOSE_STAT
858   Int_t extrarows = fComputeBG != 0 ? 8 : 0;
859 #else
860   Int_t extrarows = fComputeBG != 0 ? 3 : 0;
861 #endif
862   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);
863
864   h->GetXaxis()->SetBinLabel(kStatTriggerClass,  "Trigger class");
865   h->GetXaxis()->SetBinLabel(kStatHWTrig,        "Hardware trigger");
866   h->GetXaxis()->SetBinLabel(kStatFO1,           "FO >= 1");
867   h->GetXaxis()->SetBinLabel(kStatFO2,           "FO >= 2");
868   h->GetXaxis()->SetBinLabel(kStatV0A,           "V0A");
869   h->GetXaxis()->SetBinLabel(kStatV0C,           "V0C");
870   h->GetXaxis()->SetBinLabel(kStatFMD,           "FMD");
871   h->GetXaxis()->SetBinLabel(kStatSSD1,          "SSD >= 2");
872   h->GetXaxis()->SetBinLabel(kStatV0ABG,         "V0A BG");
873   h->GetXaxis()->SetBinLabel(kStatV0CBG,         "V0C BG");
874   h->GetXaxis()->SetBinLabel(kStatMB1,           "(FO >= 1 | V0A | V0C) & !V0 BG");
875   h->GetXaxis()->SetBinLabel(kStatMB1Prime,      "(FO >= 2 | (FO >= 1 & (V0A | V0C)) | (V0A &v0C) ) & !V0 BG");
876   h->GetXaxis()->SetBinLabel(kStatFO1AndV0,      "FO >= 1 & (V0A | V0C) & !V0 BG");
877   h->GetXaxis()->SetBinLabel(kStatV0,            "V0A & V0C & !V0 BG & !BG ID");
878   h->GetXaxis()->SetBinLabel(kStatOffline,       "Offline Trigger");
879   h->GetXaxis()->SetBinLabel(kStatAny2Hits,      "2 Hits & !V0 BG");
880   h->GetXaxis()->SetBinLabel(kStatBG,            "Background identification");
881   h->GetXaxis()->SetBinLabel(kStatAccepted,      "Accepted");
882
883   Int_t n = 1;
884   for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
885     {
886       h->GetYaxis()->SetBinLabel(n, ((TObjString*) fCollTrigClasses.At(i))->String());
887       n++;
888     }
889   for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
890     {
891       h->GetYaxis()->SetBinLabel(n, ((TObjString*) fBGTrigClasses.At(i))->String());
892       n++;
893     }
894
895   if(fComputeBG) {
896     fBGStatOffset = n;
897     h->GetYaxis()->SetBinLabel(n++, Form("BG (A+C) (Mask [0x%x])", fComputeBG));
898     h->GetYaxis()->SetBinLabel(n++, "ACC");
899 #ifdef VERBOSE_STAT
900     h->GetYaxis()->SetBinLabel(n++, "BG (A+C) %  (rel. to CINT1B)");
901     h->GetYaxis()->SetBinLabel(n++, "ACC % (rel. to CINT1B)");
902     h->GetYaxis()->SetBinLabel(n++, "ERR GOOD %");
903     h->GetYaxis()->SetBinLabel(n++, "GOOD % (rel. to 1st col)");
904     h->GetYaxis()->SetBinLabel(n++, "ERR GOOD");
905 #endif
906     h->GetYaxis()->SetBinLabel(n++, "GOOD");
907   }
908
909   return h;
910 }
911
912 void AliPhysicsSelection::Print(Option_t* /* option */) const
913 {
914   // print the configuration
915   
916   Printf("Configuration initialized for run %d (MC: %d):", fCurrentRun, fMC);
917   
918   Printf("Collision trigger classes:");
919   for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
920     Printf("%s", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
921   
922   Printf("Background trigger classes:");
923   for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
924     Printf("%s", ((TObjString*) fBGTrigClasses.At(i))->String().Data());
925
926   AliTriggerAnalysis* triggerAnalysis = dynamic_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(0));
927   
928   if (triggerAnalysis)
929   {
930     if (triggerAnalysis->GetV0TimeOffset() > 0)
931       Printf("V0 time offset active: %.2f ns", triggerAnalysis->GetV0TimeOffset());
932     
933     Printf("\nTotal available events:");
934     
935     triggerAnalysis->PrintTriggerClasses();
936   }
937   
938   if (fHistStatistics[kStatIdxAll])
939   {
940     for (Int_t i=0; i<fCollTrigClasses.GetEntries(); i++)
941     {
942       Printf("\nSelection statistics for collision trigger %s:", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
943       
944       Printf("Total events with correct trigger class: %d", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(1, i+1));
945       Printf("Selected collision candidates: %d", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(fHistStatistics[kStatIdxAll]->GetXaxis()->FindBin("Accepted"), i+1));
946     }
947   }
948   
949   if (fHistBunchCrossing)
950   {
951     Printf("\nBunch crossing statistics:");
952     
953     for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
954     {
955       TString str;
956       str.Form("Trigger %s has accepted events in the bunch crossings: ", fHistBunchCrossing->GetYaxis()->GetBinLabel(i));
957       
958       for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
959         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
960           str += Form("%d, ", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
961        
962       Printf("%s", str.Data());
963     }
964     
965     for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
966     {
967       Int_t countColl = 0;
968       Int_t countBG = 0;
969       for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
970       {
971         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
972         {
973           if (fCollTrigClasses.FindObject(fHistBunchCrossing->GetYaxis()->GetBinLabel(i)))
974             countColl++;
975           if (fBGTrigClasses.FindObject(fHistBunchCrossing->GetYaxis()->GetBinLabel(i)))
976             countBG++;
977         }
978       }
979       if (countColl > 0 && countBG > 0)
980         Printf("WARNING: Bunch crossing %d has collision and BG trigger classes active. Check BPTX functioning for this run!", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
981     }
982   }
983
984   if (fUsingCustomClasses)        
985     Printf("WARNING: Using custom trigger classes!");
986   if (fSkipTriggerClassSelection) 
987     Printf("WARNING: Skipping trigger class selection!");
988   if (fSkipV0) 
989     Printf("WARNING: Ignoring V0 information in selection");
990   if(!fBin0CallBack) 
991     Printf("WARNING: Callback not set: will not fill the statistics for the bin 0");
992
993 }
994
995 Long64_t AliPhysicsSelection::Merge(TCollection* list)
996 {
997   // Merge a list of AliMultiplicityCorrection objects with this (needed for
998   // PROOF).
999   // Returns the number of merged objects (including this).
1000
1001   if (!list)
1002     return 0;
1003
1004   if (list->IsEmpty())
1005     return 1;
1006
1007   TIterator* iter = list->MakeIterator();
1008   TObject* obj;
1009   
1010   // collections of all histograms
1011   const Int_t nHists = 9;
1012   TList collections[nHists];
1013
1014   Int_t count = 0;
1015   while ((obj = iter->Next())) {
1016
1017     AliPhysicsSelection* entry = dynamic_cast<AliPhysicsSelection*> (obj);
1018     if (entry == 0) 
1019       continue;
1020     // Update run number. If this one is not initialized (-1) take the one from 
1021     // the next physics selection to be merged with. In case of 2 different run
1022     // numbers issue a warning (should physics selections from different runs be 
1023     // merged together) A.G.
1024     Int_t currentRun = entry->GetCurrentRun();
1025     // Nothing to merge with since run number was not initialized.
1026     if (currentRun < 0) continue;
1027     if (fCurrentRun < 0) fCurrentRun = currentRun;
1028     if (fCurrentRun != currentRun)
1029        AliWarning(Form("Current run %d not matching the one to be merged with %d", fCurrentRun, currentRun));
1030     
1031     collections[0].Add(&(entry->fTriggerAnalysis));
1032     if (entry->fHistStatistics[0])
1033       collections[1].Add(entry->fHistStatistics[0]);
1034     if (entry->fHistStatistics[1])
1035       collections[2].Add(entry->fHistStatistics[1]);
1036     if (entry->fHistBunchCrossing)
1037       collections[3].Add(entry->fHistBunchCrossing);
1038     if (entry->fHistTriggerPattern)
1039       collections[4].Add(entry->fHistTriggerPattern);
1040     if (entry->fBackgroundIdentification)
1041       collections[5].Add(entry->fBackgroundIdentification);
1042
1043     count++;
1044   }
1045
1046   fTriggerAnalysis.Merge(&collections[0]);
1047   if (fHistStatistics[0])
1048     fHistStatistics[0]->Merge(&collections[1]);
1049   if (fHistStatistics[1])
1050     fHistStatistics[1]->Merge(&collections[2]);
1051   if (fHistBunchCrossing)
1052     fHistBunchCrossing->Merge(&collections[3]);
1053   if (fHistTriggerPattern)
1054     fHistTriggerPattern->Merge(&collections[4]);
1055   if (fBackgroundIdentification)
1056     fBackgroundIdentification->Merge(&collections[5]);
1057   
1058   delete iter;
1059
1060   return count+1;
1061 }
1062
1063 void AliPhysicsSelection::SaveHistograms(const char* folder) const
1064 {
1065   // write histograms to current directory
1066   
1067   if (!fHistStatistics[0] || !fHistStatistics[1])
1068     return;
1069     
1070   if (folder)
1071   {
1072     gDirectory->mkdir(folder);
1073     gDirectory->cd(folder);
1074   }
1075   
1076
1077   // Fill the last rows of fHistStatistics before saving
1078   if (fComputeBG) {
1079     Int_t triggerScheme = GetTriggerScheme(UInt_t(fCurrentRun));
1080     if(triggerScheme != 1){
1081       AliWarning("BG estimate only supported for trigger scheme \"1\" (CINT1 suite)");
1082     } else if (fUseMuonTriggers) {
1083       AliWarning("BG estimate with muon triggers to be implemented");
1084     } else {      
1085       AliInfo("BG estimate assumes that for a given run you only have A and C triggers separately or"
1086               " toghether as a AC class! Make sure this assumption holds in your case"); 
1087       
1088       // use an anum for the different trigger classes, to make loops easier to read
1089       enum {kClassB =0 , kClassA, kClassC, kClassAC, kClassE, kNClasses};
1090       const char * classFlags[] = {"B", "A", "C", "AC", "E"}; // labels
1091
1092       UInt_t * rows[kNClasses] = {0}; // Array of matching rows
1093       Int_t nrows[kNClasses] = {0};
1094       // Get rows matching the requested trigger bits for all trigger classes
1095       for(Int_t iTrigClass = 0; iTrigClass < kNClasses; iTrigClass++){
1096         nrows[iTrigClass] = GetStatRow(classFlags[iTrigClass],fComputeBG,&rows[iTrigClass]);    
1097       }
1098
1099       // 0. Determine the ratios of triggers E/B, A/B, C/B from the stat histogram
1100       // Those are used to rescale the different classes to the same number of bx ids
1101       // 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?         
1102       Int_t nBXIds[kNClasses] = {0};
1103       cout <<"Computing BG:" << endl;
1104       
1105       for(Int_t iTrigClass = 0; iTrigClass < kNClasses; iTrigClass++){
1106         for(Int_t irow = 0; irow < nrows[iTrigClass]; irow++) {
1107           if(irow==0) cout << "- Class " << classFlags[iTrigClass] << endl;   
1108           for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++) {
1109             if (fHistBunchCrossing->GetBinContent(j, rows[iTrigClass][irow]) > 0) {
1110               nBXIds[iTrigClass]++;      
1111             }
1112           }
1113           if(nBXIds[iTrigClass]>0) cout << "   Using row " << rows[iTrigClass][irow] <<  ": " 
1114                                         << fHistBunchCrossing->GetYaxis()->GetBinLabel(rows[iTrigClass][irow]) 
1115                                         << " (nBXID "<< nBXIds[iTrigClass] << ")"<< endl;
1116
1117         }
1118
1119       }
1120
1121       Float_t ratioToB[kNClasses];
1122       ratioToB[kClassE]  = nBXIds[kClassE]  >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassE]   : 0;
1123       ratioToB[kClassA]  = nBXIds[kClassA]  >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassA]   : 0;
1124       ratioToB[kClassC]  = nBXIds[kClassC]  >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassC]   : 0;
1125       ratioToB[kClassAC] = nBXIds[kClassAC] >0 ? Float_t(nBXIds[kClassB])/nBXIds[kClassAC]  : 0;
1126       Printf("Ratio between the BX ids in the different trigger classes:");
1127       Printf("  B/E  = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassE], ratioToB[kClassE] );
1128       Printf("  B/A  = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassA], ratioToB[kClassA] );
1129       Printf("  B/C  = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassC], ratioToB[kClassC] );
1130       Printf("  B/AC = %d/%d = %f", nBXIds[kClassB],nBXIds[kClassAC],ratioToB[kClassAC]);
1131       Int_t nHistStat = 2;
1132
1133       // 1. loop over all cols
1134       for(Int_t iHistStat = 0; iHistStat < nHistStat; iHistStat++){
1135         Int_t ncol = fHistStatistics[iHistStat]->GetNbinsX();
1136         Float_t good1 = 0;      
1137         for(Int_t icol = 1; icol <= ncol; icol++) {
1138           Int_t nEvents[kNClasses] = {0}; // number of events should be reset at every column
1139           // For all trigger classes, add up over row matching trigger mask (as selected before)
1140           for(Int_t iTrigClass = 0; iTrigClass < kNClasses; iTrigClass++){
1141             for(Int_t irow = 0; irow < nrows[iTrigClass]; irow++) {       
1142               nEvents[iTrigClass] += (Int_t) fHistStatistics[iHistStat]->GetBinContent(icol,rows[iTrigClass][irow]);    
1143             }
1144             //      cout << "Events " << classFlags[iTrigClass] << " ("<<icol<<") " << nEvents[iTrigClass] << endl;         
1145           }
1146           if (nEvents[kClassB]>0) {
1147             Float_t acc  = ratioToB[kClassE]*nEvents[kClassE]; 
1148             Double_t acc_err = TMath::Sqrt(ratioToB[kClassE]*ratioToB[kClassE]*nEvents[kClassE]);
1149             //      Int_t bg   = cint1A + cint1C - 2*acc;
1150             
1151              // Assuming that for a given class the triggers are either recorded as A+C or AC
1152             Float_t bg   = nEvents[kClassAC] > 0 ?
1153               fBIFactorAC*(ratioToB[kClassAC]*nEvents[kClassAC] - 2*acc):
1154               fBIFactorA* (ratioToB[kClassA]*nEvents[kClassA]-acc) + 
1155               fBIFactorC* (ratioToB[kClassC]*nEvents[kClassC]-acc) ;
1156
1157             // cout << "-----------------------" << endl;
1158             // cout << "Factors: " << fBIFactorA << " " << fBIFactorC << " " << fBIFactorAC << endl;
1159             // cout << "Ratios: "  << ratioToB[kClassA] << " " << ratioToB[kClassC] << " " << ratioToB[kClassAC] << endl;
1160             // cout << "Evts:   "  << nEvents[kClassA] << " " << nEvents[kClassC] << " " << nEvents[kClassAC] << " " <<  nEvents[kClassB] << endl;
1161             // cout << "Acc: " << acc << endl;
1162             // cout << "BG: " << bg << endl;
1163             // cout  << "  " <<   fBIFactorA* (ratioToB[kClassA]*nEvents[kClassA]-acc) <<endl;
1164             // cout  << "  " <<   fBIFactorC* (ratioToB[kClassC]*nEvents[kClassC]-acc) <<endl;
1165             // cout  << "  " <<   fBIFactorAC*(ratioToB[kClassAC]*nEvents[kClassAC] - 2*acc) << endl;
1166             // cout << "-----------------------" << endl;
1167
1168             Float_t good = Float_t(nEvents[kClassB]) - bg - acc;
1169             if (icol ==1) good1 = good;
1170             //      Float_t errGood     = TMath::Sqrt(2*(nEvents[kClassA]+nEvents[kClassC]+nEvents[kClassE]));// Error on the number of goods assuming only bg fluctuates
1171             //      DeltaG^2 = B + FA^2 A + FC^2 C + Ratio^2 (FA+FC-1)^2 E.
1172             Float_t errGood     = nEvents[kClassAC] > 0 ? 
1173               TMath::Sqrt( nEvents[kClassB] +
1174                            fBIFactorAC*fBIFactorAC*ratioToB[kClassAC]*ratioToB[kClassAC]*nEvents[kClassAC] +
1175                            ratioToB[kClassE] * ratioToB[kClassE] * 
1176                            (fBIFactorAC - 1)*(fBIFactorAC - 1)*nEvents[kClassE]) :  
1177               TMath::Sqrt( nEvents[kClassB] + 
1178                            fBIFactorA*fBIFactorA*ratioToB[kClassA]*ratioToB[kClassA]*nEvents[kClassA] +
1179                            fBIFactorC*fBIFactorC*ratioToB[kClassC]*ratioToB[kClassC]*nEvents[kClassC] +
1180                            ratioToB[kClassE] * ratioToB[kClassE] * 
1181                            (fBIFactorA + fBIFactorC - 1)*(fBIFactorA + fBIFactorC - 1)*nEvents[kClassE]);
1182
1183             Float_t errBG = nEvents[kClassAC] > 0 ? 
1184               TMath::Sqrt(fBIFactorAC*fBIFactorAC*ratioToB[kClassAC]*ratioToB[kClassAC]*nEvents[kClassAC]+
1185                           4*ratioToB[kClassE]*ratioToB[kClassE]*(fBIFactorAC*fBIFactorAC)*nEvents[kClassE]) :
1186               TMath::Sqrt(fBIFactorA*fBIFactorA*ratioToB[kClassA]*ratioToB[kClassA]*nEvents[kClassA]+
1187                           fBIFactorC*fBIFactorC*ratioToB[kClassC]*ratioToB[kClassC]*nEvents[kClassC]+
1188                           ratioToB[kClassE]*ratioToB[kClassE]*(fBIFactorA+fBIFactorC)*(fBIFactorA+fBIFactorC)*nEvents[kClassE]);
1189         
1190         
1191             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowBG,bg);        
1192             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowBG,errBG);     
1193             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowAcc,acc);      
1194             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowAcc,acc_err);  
1195             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowGood,good);    
1196             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowGood,errGood);    
1197
1198 #ifdef VERBOSE_STAT
1199             //kStatRowBG=0,kStatRowAcc,kStatRowBGFrac,kStatRowAccFrac,kStatRowErrGoodFrac,kStatRowGoodFrac,kStatRowGood,kStatRowErrGood
1200             Float_t accFrac   = Float_t(acc) / nEvents[kClassB]  *100;
1201             Float_t errAccFrac= Float_t(acc_err) / nEvents[kClassB]  *100;
1202             Float_t bgFrac    = Float_t(bg)  / nEvents[kClassB]  *100;
1203             Float_t goodFrac  = Float_t(good)  / good1 *100;
1204             Float_t errGoodFrac = errGood/good1 * 100;
1205             Float_t errFracBG = bg > 0 ? TMath::Sqrt((errBG/bg)*(errBG/bg) + 1/nEvents[kClassB])*bgFrac : 0;
1206             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowBGFrac,bgFrac);        
1207             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowBGFrac,errFracBG);     
1208             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowAccFrac,accFrac);    
1209             fHistStatistics[iHistStat]->SetBinError  (icol,fBGStatOffset+kStatRowAccFrac,errAccFrac);    
1210             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowGoodFrac,goodFrac);    
1211             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowErrGoodFrac,errGoodFrac);    
1212             fHistStatistics[iHistStat]->SetBinContent(icol,fBGStatOffset+kStatRowErrGood,errGood);    
1213 #endif
1214           }
1215         }
1216       }
1217     }  
1218   }
1219
1220   fHistStatistics[0]->Write();
1221   fHistStatistics[1]->Write();
1222   if(fHistBunchCrossing ) fHistBunchCrossing ->Write();
1223   if(fHistTriggerPattern) fHistTriggerPattern->Write();
1224   
1225   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
1226   for (Int_t i=0; i < count; i++)
1227   {
1228     TString triggerClass = "trigger_histograms_";
1229     if (i < fCollTrigClasses.GetEntries())
1230       triggerClass += ((TObjString*) fCollTrigClasses.At(i))->String();
1231     else
1232       triggerClass += ((TObjString*) fBGTrigClasses.At(i - fCollTrigClasses.GetEntries()))->String();
1233   
1234     gDirectory->mkdir(triggerClass);
1235     gDirectory->cd(triggerClass);
1236   
1237     static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i))->SaveHistograms();
1238     
1239     gDirectory->cd("..");
1240   }
1241  
1242   if (fBackgroundIdentification)
1243   {
1244     gDirectory->mkdir("background_identification");
1245     gDirectory->cd("background_identification");
1246       
1247     fBackgroundIdentification->GetOutput()->Write();
1248       
1249     gDirectory->cd("..");
1250   }
1251   
1252   if (folder)
1253     gDirectory->cd("..");
1254 }
1255
1256 Int_t AliPhysicsSelection::GetStatRow(const char * triggerBXClass, UInt_t offlineTriggerType, UInt_t ** rowIDs) const {
1257   // Puts inside the array rowIDs the row number for a given offline
1258   // trigger in a given bx class. Returns the total number of lines
1259   // matching the selection
1260   // triggerBXClass can be either "A", "AC", "B" or "E"
1261   // offlineTriggerType is one of the types defined in AliVEvent
1262   // User should delete rowIDs if no longer needed
1263
1264   if(!fHistStatistics[0]) {
1265     AliWarning("Not initialized, returning 0");
1266     return 0;
1267   }
1268   const Int_t nrows = fHistStatistics[0]->GetNbinsY();
1269
1270   // allocate memory for at maximum nrows
1271   Int_t nMatches = 0;
1272   (*rowIDs) = new UInt_t[nrows];
1273
1274   // Build regular expression. look for a +, followed by the beginning
1275   // of a word. Within the word, look for the class id after any
1276   // number of any char, but at most one dash ("[^-]*-?"), followed by
1277   // a - and then any char (".*") and at the class id ("\\&(\\d)")
1278   // The class id is stored.
1279   // WARNING: please check this if the trigger classes change
1280   TPRegexp re(Form("\\+\\b[^-]*-?%s-.*\\&(\\d)",triggerBXClass));  
1281   // Loop over rows and find matching ones:
1282   for(Int_t irow = 1; irow <= nrows; irow++){
1283     TObjArray * matches = re.MatchS(fHistStatistics[0]->GetYaxis()->GetBinLabel(irow));
1284     if (matches->GetEntries()) {
1285       TString s = ((TObjString*)matches->At(1))->GetString();      
1286       if(UInt_t(s.Atoi()) & offlineTriggerType) { // bitwise comparison with the requested mask
1287         //      cout << "Marching " << s.Data() << " " << offlineTriggerType << " " << fHistStatistics[0]->GetYaxis()->GetBinLabel(irow) << endl;       
1288         (*rowIDs)[nMatches] = irow;
1289         nMatches++;     
1290       }
1291     }
1292     delete matches;
1293   }
1294   
1295   return nMatches;
1296 }
1297
1298 void AliPhysicsSelection::SetBIFactors(const AliESDEvent * aESD) {
1299   // Set factors for realtive bunch intesities
1300   if(!aESD) { 
1301     AliFatal("ESD not given");
1302   }
1303   Int_t run = aESD->GetRunNumber();
1304   if (run > 105268) {
1305     // intensities stored in the ESDs
1306     const AliESDRun* esdRun = aESD->GetESDRun();
1307     Double_t intAB = esdRun->GetMeanIntensityIntecting(0);
1308     Double_t intCB = esdRun->GetMeanIntensityIntecting(1);
1309     Double_t intAA = esdRun->GetMeanIntensityNonIntecting(0);
1310     Double_t intCC = esdRun->GetMeanIntensityNonIntecting(1);
1311
1312     // cout << "INT " <<intAB <<endl;
1313     // cout << "INT " <<intCB <<endl;
1314     // cout << "INT " <<intAA <<endl;
1315     // cout << "INT " <<intCC <<endl;
1316
1317     if (intAB > -1 && intAA > -1) {
1318       fBIFactorA = intAB/intAA;
1319     } else {
1320       AliWarning("Cannot set fBIFactorA, assuming 1");
1321     }
1322     
1323     if (intCB > -1 && intCC > -1) {
1324       fBIFactorC = intCB/intCC;
1325     } else {
1326       AliWarning("Cannot set fBIFactorC, assuming 1");
1327     }
1328       
1329     if (intAB > -1 && intAA > -1 &&
1330         intCB > -1 && intCC > -1) {
1331       fBIFactorAC = (intAB+intCB)/(intAA+intCC);
1332     } else {
1333       AliWarning("Cannot set fBIFactorAC, assuming 1");
1334     }
1335         
1336   }  
1337   else {
1338     // First runs. Intensities hardcoded
1339     switch(run) {
1340     case 104155:
1341       fBIFactorA = 0.961912722908;
1342       fBIFactorC = 1.04992336081;
1343       break;
1344     case 104157:
1345       fBIFactorA = 0.947312854998;
1346       fBIFactorC = 1.01599706417;
1347       break;
1348     case 104159:
1349       fBIFactorA = 0.93659320151;
1350       fBIFactorC = 0.98580804207;
1351       break;
1352     case 104160:
1353       fBIFactorA = 0.929664189926;
1354       fBIFactorC = 0.963467679851;
1355       break;
1356     case 104315:
1357       fBIFactorA = 1.08939104979;
1358       fBIFactorC = 0.931113921925;
1359       break;
1360     case 104316:
1361       fBIFactorA = 1.08351880974;
1362       fBIFactorC = 0.916068345845;
1363       break;
1364     case 104320:
1365       fBIFactorA = 1.07669281245;
1366       fBIFactorC = 0.876818744763;
1367       break;
1368     case 104321:
1369       fBIFactorA = 1.00971079602;
1370       fBIFactorC = 0.773781299076;
1371       break;
1372     case 104792:
1373       fBIFactorA = 0.787215863962;
1374       fBIFactorC = 0.778253173071;
1375       break;
1376     case 104793:
1377       fBIFactorA = 0.692211363661;
1378       fBIFactorC = 0.733152456667;
1379       break;
1380     case 104799:
1381       fBIFactorA = 1.04027825161;
1382       fBIFactorC = 1.00530825942;
1383       break;
1384     case 104800:
1385       fBIFactorA = 1.05309910671;
1386       fBIFactorC = 1.00376801855;
1387       break;
1388     case 104801:
1389       fBIFactorA = 1.0531231922;
1390       fBIFactorC = 0.992439666758;
1391       break;
1392     case 104802:
1393       fBIFactorA = 1.04191478134;
1394       fBIFactorC = 0.979368585208;
1395       break;
1396     case 104803:
1397       fBIFactorA = 1.03121314094;
1398       fBIFactorC = 0.973379962609;
1399       break;
1400     case 104824:
1401       fBIFactorA = 0.969945926722;
1402       fBIFactorC = 0.39549745806;
1403       break;
1404     case 104825:
1405       fBIFactorA = 0.968627213937;
1406       fBIFactorC = 0.310100412205;
1407       break;
1408     case 104841:
1409       fBIFactorA = 0.991601393212;
1410       fBIFactorC = 0.83762204722;
1411       break;
1412     case 104845:
1413       fBIFactorA = 0.98040863886;
1414       fBIFactorC = 0.694824205793;
1415       break;
1416     case 104867:
1417       fBIFactorA = 1.10646173412;
1418       fBIFactorC = 0.841407246916;
1419       break;
1420     case 104876:
1421       fBIFactorA = 1.12063452421;
1422       fBIFactorC = 0.78726542895;
1423       break;
1424     case 104890:
1425       fBIFactorA = 1.02346137453;
1426       fBIFactorC = 1.03355663595;
1427       break;
1428     case 104892:
1429       fBIFactorA = 1.05406025913;
1430       fBIFactorC = 1.00029166135;
1431       break;
1432     case 105143:
1433       fBIFactorA = 0.947343384349;
1434       fBIFactorC = 0.972637444408;
1435       break;
1436     case 105160:
1437       fBIFactorA = 0.908854622177;
1438       fBIFactorC = 0.958851103977;
1439       break; 
1440     case 105256:
1441       fBIFactorA = 0.810076150206;
1442       fBIFactorC = 0.884663561883;
1443       break;
1444     case 105257:
1445       fBIFactorA = 0.80974912303;
1446       fBIFactorC = 0.878859123479;
1447       break;
1448     case 105268:
1449       fBIFactorA = 0.809052110679;
1450       fBIFactorC = 0.87233890989;
1451       break;
1452     default:
1453       fBIFactorA = 1;
1454       fBIFactorC = 1;
1455     }
1456   }
1457
1458 }