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