Extending the functionality of the physics selection. Instead of providing a yes...
[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
116 ClassImp(AliPhysicsSelection)
117
118 AliPhysicsSelection::AliPhysicsSelection() :
119   AliAnalysisCuts("AliPhysicsSelection", "AliPhysicsSelection"),
120   fCurrentRun(-1),
121   fMC(kFALSE),
122   fCollTrigClasses(),
123   fBGTrigClasses(),
124   fTriggerAnalysis(),
125   fBackgroundIdentification(0),
126   fHistBunchCrossing(0),
127   fHistTriggerPattern(0),
128   fSkipTriggerClassSelection(0),
129   fUsingCustomClasses(0),
130   fSkipV0(0),
131   fBIFactorA(1),
132   fBIFactorC(1),
133   fComputeBG(0),
134   fUseBXNumbers(1),
135   fUseMuonTriggers(0),
136   fFillingScheme(""),
137   fBin0CallBack(""),
138   fBin0CallBackPointer(0)
139 {
140   // constructor
141   
142   fCollTrigClasses.SetOwner(1);
143   fBGTrigClasses.SetOwner(1);
144   fTriggerAnalysis.SetOwner(1);
145   fHistStatistics[0] = 0;
146   fHistStatistics[1] = 0;
147   
148   AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kWarning);
149 }
150     
151 AliPhysicsSelection::~AliPhysicsSelection()
152 {
153   // destructor
154
155   fCollTrigClasses.Delete();
156   fBGTrigClasses.Delete();
157   fTriggerAnalysis.Delete();
158
159   if (fHistStatistics[0])
160   {
161     delete fHistStatistics[0];
162     fHistStatistics[0] = 0;
163   }
164   if (fHistStatistics[1])
165   {
166     delete fHistStatistics[1];
167     fHistStatistics[1] = 0;
168   }
169   
170   if (fHistBunchCrossing)
171   {
172     delete fHistBunchCrossing;
173     fHistBunchCrossing = 0;
174   }
175   if (fHistTriggerPattern)
176   {
177     delete fHistTriggerPattern;
178     fHistTriggerPattern = 0;
179   }
180
181 }
182
183 UInt_t AliPhysicsSelection::CheckTriggerClass(const AliESDEvent* aEsd, const char* trigger) const
184 {
185   // checks if the given trigger class(es) are found for the current event
186   // format of trigger: +TRIGGER1 -TRIGGER2 [#XXX] [&YY]
187   //   requires TRIGGER1 and rejects TRIGGER2
188   //   in bunch crossing XXX
189   //   if successful, a word with bit YY set is returned (for association between entry in fCollTrigClasses and AliVEvent::EOfflineTriggerTypes)
190   
191   Bool_t foundBCRequirement = kFALSE;
192   Bool_t foundCorrectBC = kFALSE;
193   
194   UInt_t returnCode = AliVEvent::kUserDefined;
195   
196   TString str(trigger);
197   TObjArray* tokens = str.Tokenize(" ");
198   
199   for (Int_t i=0; i < tokens->GetEntries(); i++)
200   {
201     TString str2(((TObjString*) tokens->At(i))->String());
202     
203     if (str2[0] == '+' || str2[0] == '-')
204     {
205       Bool_t flag = (str2[0] == '+');
206       
207       str2.Remove(0, 1);
208       
209       if (flag && !aEsd->IsTriggerClassFired(str2))
210       {
211         AliDebug(AliLog::kDebug, Form("Rejecting event because trigger class %s is not present", str2.Data()));
212         delete tokens;
213         return kFALSE;
214       }
215       if (!flag && aEsd->IsTriggerClassFired(str2))
216       {
217         AliDebug(AliLog::kDebug, Form("Rejecting event because trigger class %s is present", str2.Data()));
218         delete tokens;
219         return kFALSE;
220       }
221     }
222     else if (str2[0] == '#')
223     {
224       foundBCRequirement = kTRUE;
225     
226       str2.Remove(0, 1);
227       
228       Int_t bcNumber = str2.Atoi();
229       AliDebug(AliLog::kDebug, Form("Checking for bunch crossing number %d", bcNumber));
230       
231       if (aEsd->GetBunchCrossNumber() == bcNumber)
232       {
233         foundCorrectBC = kTRUE;
234         AliDebug(AliLog::kDebug, Form("Found correct bunch crossing %d", bcNumber));
235       }
236     }
237     else if (str2[0] == '&' && !fUsingCustomClasses)
238     {
239       str2.Remove(0, 1);
240       
241       returnCode = 1 << str2.Atoi();
242     }
243     else
244       AliFatal(Form("Invalid trigger syntax: %s", trigger));
245   }
246   
247   delete tokens;
248   
249   if (foundBCRequirement && !foundCorrectBC)
250     return kFALSE;
251   
252   return returnCode;
253 }
254     
255 UInt_t AliPhysicsSelection::IsCollisionCandidate(const AliESDEvent* aEsd)
256 {
257   // checks if the given event is a collision candidate
258   //
259   // returns a bit word describing the fired offline triggers (see AliVEvent::EOfflineTriggerTypes)
260   
261   if (fCurrentRun != aEsd->GetRunNumber())
262     if (!Initialize(aEsd->GetRunNumber()))
263       AliFatal(Form("Could not initialize for run %d", aEsd->GetRunNumber()));
264     
265   const AliESDHeader* esdHeader = aEsd->GetHeader();
266   if (!esdHeader)
267   {
268     AliError("ESD Header could not be retrieved");
269     return kFALSE;
270   }
271   
272   // check event type; should be PHYSICS = 7 for data and 0 for MC
273   if (!fMC)
274   {
275     if (esdHeader->GetEventType() != 7)
276       return kFALSE;
277   }
278   else
279   {
280     if (esdHeader->GetEventType() != 0)
281       AliFatal(Form("Invalid event type for MC: %d", esdHeader->GetEventType()));
282   }
283   
284   UInt_t accept = 0;
285     
286   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
287   for (Int_t i=0; i < count; i++)
288   {
289     const char* triggerClass = 0;
290     if (i < fCollTrigClasses.GetEntries())
291       triggerClass = ((TObjString*) fCollTrigClasses.At(i))->String();
292     else
293       triggerClass = ((TObjString*) fBGTrigClasses.At(i - fCollTrigClasses.GetEntries()))->String();
294   
295     AliDebug(AliLog::kDebug, Form("Processing trigger class %s", triggerClass));
296   
297     AliTriggerAnalysis* triggerAnalysis = static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i));
298   
299     triggerAnalysis->FillTriggerClasses(aEsd);
300     
301     UInt_t singleTriggerResult = CheckTriggerClass(aEsd, triggerClass);
302     if (singleTriggerResult)
303     {
304       triggerAnalysis->FillHistograms(aEsd);
305   
306       Bool_t isBin0 = kFALSE;
307       if (fBin0CallBack != "") {
308         AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
309         if (!mgr) {
310           AliError("Cannot get the analysis manager");
311         }  
312         else {
313           isBin0 = ((AliAnalysisTaskSE*)mgr->GetTask(fBin0CallBack.Data()))->IsEventInBinZero();
314         }
315       } else if (fBin0CallBackPointer) {
316           isBin0 = (*fBin0CallBackPointer)(aEsd);
317         
318       }
319       
320
321       
322       // hardware trigger (should only remove events for MC)
323       // replay CINT1B hardware trigger
324       // TODO this has to depend on the actual hardware trigger (and that depends on the run...)
325       Int_t fastORHW = triggerAnalysis->SPDFiredChips(aEsd, 1); // SPD number of chips from trigger bits (!)
326       Bool_t v0A       = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0A);
327       Bool_t v0C       = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0C);
328       Bool_t v0AHW     = fSkipV0 ? 0 :(triggerAnalysis->V0Trigger(aEsd, AliTriggerAnalysis::kASide, kTRUE) == AliTriggerAnalysis::kV0BB);// should replay hw trigger
329       Bool_t v0CHW     = fSkipV0 ? 0 :(triggerAnalysis->V0Trigger(aEsd, AliTriggerAnalysis::kCSide, kTRUE) == AliTriggerAnalysis::kV0BB);// should replay hw trigger
330       // offline trigger
331       Int_t fastOROffline = triggerAnalysis->SPDFiredChips(aEsd, 0); // SPD number of chips from clusters (!)
332       Bool_t v0ABG = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0ABG);
333       Bool_t v0CBG = fSkipV0 ? 0 :triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0CBG);
334       Bool_t v0BG = v0ABG || v0CBG;
335
336       // fmd
337       Bool_t fmdA =  triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kFMDA);
338       Bool_t fmdC =  triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kFMDC);
339       Bool_t fmd  = fmdA || fmdC;
340     
341       // SSD
342       Int_t ssdClusters = triggerAnalysis->SSDClusters(aEsd);
343
344       // Some "macros"
345       Bool_t mb1 = (fastOROffline > 0 || v0A || v0C) && (!v0BG);
346       Bool_t mb1prime = (fastOROffline > 1 || (fastOROffline > 0 && (v0A || v0C)) || (v0A && v0C) ) && (!v0BG);
347
348       // Background rejection
349       Bool_t bgID = kFALSE;
350       if (fBackgroundIdentification)
351         bgID = ! fBackgroundIdentification->IsSelected(const_cast<AliESDEvent*> (aEsd));
352       
353       Int_t ntrig = fastOROffline; // any 2 hits
354       if(v0A)              ntrig += 1;
355       if(v0C)              ntrig += 1; //v0C alone is enough
356       if(fmd)              ntrig += 1;
357       if(ssdClusters>1)    ntrig += 1;
358
359
360       Bool_t hwTrig = fastORHW > 0 || v0AHW || v0CHW;
361
362       // Fill trigger pattern histo
363       Int_t tpatt = 0;
364       if (fastORHW>0) tpatt+=1;
365       if (v0AHW)      tpatt+=2;
366       if (v0CHW)      tpatt+=4;
367       fHistTriggerPattern->Fill( tpatt );
368
369       // fill statistics and return decision
370       const Int_t nHistStat = 2;
371       for(Int_t iHistStat = 0; iHistStat < nHistStat; iHistStat++){
372         if (iHistStat == kStatIdxBin0 && !isBin0) continue; // skip the filling of bin0 stats if the event is not in the bin0
373       
374         fHistStatistics[iHistStat]->Fill(kStatTriggerClass, i);
375
376
377
378         // We fill the rest only if hw trigger is ok
379         if (!hwTrig)
380           {
381             AliDebug(AliLog::kDebug, "Rejecting event because hardware trigger is not fired");
382             continue;
383           } else {       
384           fHistStatistics[iHistStat]->Fill(kStatHWTrig, i);
385         }
386       
387
388         // v0 BG stats
389         if (v0ABG)
390           fHistStatistics[iHistStat]->Fill(kStatV0ABG, i);
391         if (v0CBG)
392           fHistStatistics[iHistStat]->Fill(kStatV0CBG, i);
393
394         // We fill the rest only if mb1 && ! v0BG
395         if (mb1)
396           fHistStatistics[iHistStat]->Fill(kStatMB1, i);
397         else continue;
398
399         if (mb1prime)
400           fHistStatistics[iHistStat]->Fill(kStatMB1Prime, i);
401
402         if (fmd)
403           fHistStatistics[iHistStat]->Fill(kStatFMD, i);
404
405         if(ssdClusters>1)
406           fHistStatistics[iHistStat]->Fill(kStatSSD1, i);
407
408         if(ntrig >= 2 && !v0BG) 
409           fHistStatistics[iHistStat]->Fill(kStatAny2Hits, i);
410
411         if (fastOROffline > 0)
412           fHistStatistics[iHistStat]->Fill(kStatFO1, i);
413         if (fastOROffline > 1)
414           fHistStatistics[iHistStat]->Fill(kStatFO2, i);
415         
416         if (v0A)
417           fHistStatistics[iHistStat]->Fill(kStatV0A, i);
418         if (v0C)
419           fHistStatistics[iHistStat]->Fill(kStatV0C, i);
420         
421         //       if (fastOROffline > 1 && !v0BG)
422         //         fHistStatistics[iHistStat]->Fill(kStatFO2NoBG, i);
423             
424         if (fastOROffline > 0 && (v0A || v0C) && !v0BG)
425           fHistStatistics[iHistStat]->Fill(kStatFO1AndV0, i);
426   
427         if (v0A && v0C && !v0BG && !bgID)
428           fHistStatistics[iHistStat]->Fill(kStatV0, i);
429
430
431         if ( mb1 )
432         
433           {
434             if (!v0BG || fSkipV0)
435               {
436                 if (!v0BG) fHistStatistics[iHistStat]->Fill(kStatOffline, i);
437       
438                 if (fBackgroundIdentification && bgID)
439                   {
440                     AliDebug(AliLog::kDebug, "Rejecting event because of background identification");
441                     fHistStatistics[iHistStat]->Fill(kStatBG, i);
442                   }
443                 else
444                   {
445                     AliDebug(AliLog::kDebug, "Accepted event for histograms");
446             
447                     fHistStatistics[iHistStat]->Fill(kStatAccepted, i);
448                     if(iHistStat == kStatIdxAll) fHistBunchCrossing->Fill(aEsd->GetBunchCrossNumber(), i); // Fill only for all (avoid double counting)
449                     if((i < fCollTrigClasses.GetEntries() || fSkipTriggerClassSelection) && (iHistStat==kStatIdxAll))
450                       accept |= singleTriggerResult; // only set for "all" (should not really matter)
451                   }
452               }
453             else
454               AliDebug(AliLog::kDebug, "Rejecting event because of V0 BG flag");
455           }
456         else
457           AliDebug(AliLog::kDebug, "Rejecting event because trigger condition is not fulfilled");
458       }
459     }
460   }
461  
462   if (accept)
463     AliDebug(AliLog::kDebug, Form("Accepted event as collision candidate with bit mask %d", accept));
464   
465   return accept;
466 }
467
468 Int_t AliPhysicsSelection::GetTriggerScheme(UInt_t runNumber) const
469 {
470   // returns the current trigger scheme (classes that are accepted/rejected)
471   
472   if (fMC)
473     return 0;
474     
475   // TODO dependent on run number
476   
477   switch (runNumber)
478   {
479     // CSMBB triggers
480     case 104044:
481     case 105054:
482     case 105057:
483       return 2;
484   }
485   
486   // default: CINT1 suite
487   return 1;
488 }  
489
490 const char * AliPhysicsSelection::GetFillingScheme(UInt_t runNumber)  {
491
492   if(fMC) return "MC";
493
494   if      (runNumber >= 104065 && runNumber <= 104160) {
495     return "4x4a";
496   } 
497   else if (runNumber >= 104315 && runNumber <= 104321) {
498     return "4x4a*";
499   }
500   else if (runNumber >= 104792 && runNumber <= 104803) {
501     return "4x4b";
502   }
503   else if (runNumber >= 104824 && runNumber <= 104892) {
504     return "4x4c";
505   }
506   else if (runNumber == 105143 || runNumber == 105160) {
507     return "16x16a";
508   }
509   else if (runNumber >= 105256 && runNumber <= 105268) {
510     return "4x4c";
511   } 
512   else if (runNumber >= 114786 && runNumber <= 116684) {
513     return "Single_2b_1_1_1";
514   }
515   else if (runNumber >= 117048 && runNumber <= 117120) {
516     return "Single_3b_2_2_2";
517   }
518   else if (runNumber >= 117220 && runNumber <= 119163) {
519     return "Single_2b_1_1_1";
520   }
521   else if (runNumber >= 119837 && runNumber <= 119862) {
522     return "Single_4b_2_2_2";
523   }
524   else if (runNumber >= 119902 && runNumber <= 120691) {
525     return "Single_6b_3_3_3";
526   }
527   else if (runNumber >= 120741 && runNumber <= 122375) {
528     return "Single_13b_8_8_8";
529   }
530   else {
531     AliError(Form("Unknown filling scheme (run %d)", runNumber));
532   }
533
534   return "Unknown";
535 }
536
537 const char * AliPhysicsSelection::GetBXIDs(UInt_t runNumber, const char * trigger)  {
538
539   if (!fUseBXNumbers || fMC) return "";
540
541   if      (runNumber >= 104065 && runNumber <= 104160) {
542     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #2128 #3019";
543     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #346 #3465";
544     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1234 #1680";
545     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
546     else AliError(Form("Unknown trigger: %s", trigger));
547   } 
548   else if (runNumber >= 104315 && runNumber <= 104321) {
549     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #2000 #2891";
550     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #218 #3337";
551     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1106 #1552";
552     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
553     else AliError(Form("Unknown trigger: %s", trigger));
554   }
555   else if (runNumber >= 104792 && runNumber <= 104803) {
556     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #2228 #3119";
557     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2554 #446";
558     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1334 #769";
559     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
560     else AliError(Form("Unknown trigger: %s", trigger));
561   }
562   else if (runNumber >= 104824 && runNumber <= 104892) {
563     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #3119 #769";
564     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2554 #446";
565     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1334 #2228";
566     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
567     else AliError(Form("Unknown trigger: %s", trigger));
568   }
569   else if (runNumber == 105143 || runNumber == 105160) {
570     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #1337 #1418 #2228 #2309 #3119 #3200 #446 #527";
571     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #1580  #1742  #1904  #2066  #2630  #2792  #2954  #3362";
572     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "  #845  #1007  #1169   #1577 #3359 #3521 #119  #281 ";
573     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
574     else AliError(Form("Unknown trigger: %s", trigger));
575   }
576   else if (runNumber >= 105256 && runNumber <= 105268) {
577     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #3019 #669";
578     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2454 #346";
579     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #1234 #2128";
580     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #790";
581     else AliError(Form("Unknown trigger: %s", trigger));
582   } else if (runNumber >= 114786 && runNumber <= 116684) { // 7 TeV 2010, assume always the same filling scheme
583     if     (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return " #346";
584     else if(!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return " #2131";
585     else if(!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return " #3019";
586     else if(!strcmp("CINT1-E-NOPF-ALL",trigger))     return " #1238";
587     else AliError(Form("Unknown trigger: %s", trigger));
588   }
589   else if (runNumber >= 117048 && runNumber <= 117120) {
590     //    return "Single_3b_2_2_2";
591    if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346  #1240 ";
592    else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #2131 ";
593    else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019 ";
594    else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1238";
595    else AliError(Form("Unknown trigger: %s", trigger));
596
597   }
598   else if (runNumber >= 117220 && runNumber <= 119163) {
599     //    return "Single_2b_1_1_1";
600     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346 ";
601     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #2131 ";
602     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019 ";
603     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1238 ";
604     else AliError(Form("Unknown trigger: %s", trigger));                                                    
605   }
606   else if (runNumber >= 119837 && runNumber <= 119862) {
607     //    return "Single_4b_2_2_2";
608     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #669  #3019 ";
609     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #346  #2454 ";
610     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #1234  #2128 ";
611     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1681 #3463";
612     else AliError(Form("Unknown trigger: %s", trigger));
613
614   }
615   else if (runNumber >= 119902 && runNumber <= 120691) {
616     //    return "Single_6b_3_3_3";
617     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346  #546  #746 ";
618     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #2131  #2331  #2531 ";
619     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019  #3219  #3419 ";
620     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1296 #1670";
621     else AliError(Form("Unknown trigger: %s", trigger));
622   }
623   else if (runNumber >= 120741 && runNumber <= 122375) {
624     //    return "Single_13b_8_8_8";
625     if      (!strcmp("CINT1B-ABCE-NOPF-ALL",trigger)) return "   #346  #446  #546  #646  #1240  #1340  #1440  #1540 ";
626     else if (!strcmp("CINT1A-ABCE-NOPF-ALL",trigger)) return "   #946  #2131  #2231  #2331  #2431 ";
627     else if (!strcmp("CINT1C-ABCE-NOPF-ALL",trigger)) return "   #3019  #3119  #3219  #3319  #3519 ";
628     else if (!strcmp("CINT1-E-NOPF-ALL",trigger)) return " #1835 #2726";
629     else AliError(Form("Unknown trigger: %s", trigger));
630     
631   }
632
633   else {
634     AliError(Form("Unknown run %d, using all BXs!",runNumber));
635   }
636
637   return "";
638 }
639     
640 Bool_t AliPhysicsSelection::Initialize(Int_t runNumber)
641 {
642   // initializes the object for the given run
643   // 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
644   
645   Bool_t oldStatus = TH1::AddDirectoryStatus();
646   TH1::AddDirectory(kFALSE);
647   
648   if(!fBin0CallBack) 
649     AliError("Bin0 Callback not set: will not fill the statistics for the bin 0");
650
651   if (fMC) {
652     // override BX and bg options in case of MC
653     fComputeBG    = kFALSE;
654     fUseBXNumbers = kFALSE;
655   }
656
657   Int_t triggerScheme = GetTriggerScheme(runNumber);
658   if (!fUsingCustomClasses && fCurrentRun != -1 && triggerScheme != GetTriggerScheme(fCurrentRun))
659     AliFatal("Processing several runs with different trigger schemes is not supported");
660   
661   if(fComputeBG && fCurrentRun != -1 && fCurrentRun != runNumber) 
662     AliFatal("Cannot process several runs because BG computation is requested");
663
664   if(fComputeBG && !fUseBXNumbers) 
665     AliFatal("Cannot compute BG id BX numbers are not used");
666   
667   if(fUseBXNumbers && fFillingScheme != "" && fFillingScheme != GetFillingScheme(runNumber))
668     AliFatal("Cannot process runs with different filling scheme if usage of BX numbers is requested");
669
670   fFillingScheme      = GetFillingScheme(runNumber);
671
672   if(fComputeBG) SetBIFactors(runNumber);
673
674   AliInfo(Form("Initializing for run %d", runNumber));
675   
676   // initialize first time?
677   if (fCurrentRun == -1)
678   {
679     if (fUsingCustomClasses) {
680       AliInfo("Using user-provided trigger classes");
681     } else {
682       switch (triggerScheme)
683       {
684       case 0:
685         fCollTrigClasses.Add(new TObjString("&0"));
686         break;
687         
688       case 1:
689         fCollTrigClasses.Add(new TObjString(Form("%s%s &0","+CINT1B-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1B-ABCE-NOPF-ALL"))));
690         fBGTrigClasses.Add  (new TObjString(Form("%s%s &0","+CINT1A-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1A-ABCE-NOPF-ALL"))));
691         fBGTrigClasses.Add  (new TObjString(Form("%s%s &0","+CINT1C-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1C-ABCE-NOPF-ALL"))));
692         fBGTrigClasses.Add  (new TObjString(Form("%s%s &0","+CINT1-E-NOPF-ALL",      GetBXIDs(runNumber,"CINT1-E-NOPF-ALL"))));
693
694         // Muon trigger have the same BXIDs of the corresponding CINT triggers
695         fCollTrigClasses.Add(new TObjString(Form("%s%s &2","+CMUS1B-ABCE-NOPF-MUON",  GetBXIDs(runNumber,"CINT1B-ABCE-NOPF-ALL"))));
696         fBGTrigClasses.Add  (new TObjString(Form("%s%s &2","+CMUS1A-ABCE-NOPF-MUON",  GetBXIDs(runNumber,"CINT1A-ABCE-NOPF-ALL"))));
697         fBGTrigClasses.Add  (new TObjString(Form("%s%s &2","+CMUS1C-ABCE-NOPF-MUON",  GetBXIDs(runNumber,"CINT1C-ABCE-NOPF-ALL"))));        
698         fBGTrigClasses.Add  (new TObjString(Form("%s%s &2","+CMUS1-E-NOPF-MUON"    ,  GetBXIDs(runNumber,"CINT1-E-NOPF-ALL"))));
699         break;
700         
701       case 2:
702         fCollTrigClasses.Add(new TObjString("+CSMBB-ABCE-NOPF-ALL &0"));
703         fBGTrigClasses.Add(new TObjString("+CSMBA-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL &0"));
704         fBGTrigClasses.Add(new TObjString("+CSMBC-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL &0"));
705         break;
706         
707       case 3:
708         // 
709         break;
710         
711       default:
712         AliFatal(Form("Unsupported trigger scheme %d", triggerScheme));
713       }
714     }
715     
716     Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
717     
718     for (Int_t i=0; i<count; i++)
719     {
720       AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis;
721       triggerAnalysis->SetAnalyzeMC(fMC);
722       triggerAnalysis->EnableHistograms();
723       triggerAnalysis->SetSPDGFOThreshhold(1);
724       fTriggerAnalysis.Add(triggerAnalysis);
725     }
726       
727     // TODO: shall I really delete this?
728     if (fHistStatistics[0])
729       delete fHistStatistics[0];
730     if (fHistStatistics[1])
731       delete fHistStatistics[1];
732   
733     fHistStatistics[kStatIdxBin0] = BookHistStatistics("_Bin0");
734     fHistStatistics[kStatIdxAll]  = BookHistStatistics("");
735     
736     if (fHistBunchCrossing)
737       delete fHistBunchCrossing;
738   
739     fHistBunchCrossing = new TH2F("fHistBunchCrossing", "fHistBunchCrossing;bunch crossing number;", 4000, -0.5, 3999.5,  count, -0.5, -0.5 + count);
740
741     if (fHistTriggerPattern)
742       delete fHistTriggerPattern;
743     
744     const int ntrig=3;
745     Int_t n = 1;
746     const Int_t nbinTrig = TMath::Nint(TMath::Power(2,ntrig));
747
748     fHistTriggerPattern = new TH1F("fHistTriggerPattern", "Trigger pattern: FO + 2*v0A + 4*v0C", 
749                                    nbinTrig, -0.5, nbinTrig-0.5);    
750     fHistTriggerPattern->GetXaxis()->SetBinLabel(1,"NO TRIG");
751     fHistTriggerPattern->GetXaxis()->SetBinLabel(2,"FO");
752     fHistTriggerPattern->GetXaxis()->SetBinLabel(3,"v0A");
753     fHistTriggerPattern->GetXaxis()->SetBinLabel(4,"FO & v0A");
754     fHistTriggerPattern->GetXaxis()->SetBinLabel(5,"v0C");
755     fHistTriggerPattern->GetXaxis()->SetBinLabel(6,"FO & v0C");
756     fHistTriggerPattern->GetXaxis()->SetBinLabel(7,"v0A & v0C");
757     fHistTriggerPattern->GetXaxis()->SetBinLabel(8,"FO & v0A & v0C");
758
759   
760     n = 1;
761     for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
762     {
763       fHistBunchCrossing->GetYaxis()->SetBinLabel(n, ((TObjString*) fCollTrigClasses.At(i))->String());
764       n++;
765     }
766     for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
767     {
768       fHistBunchCrossing->GetYaxis()->SetBinLabel(n, ((TObjString*) fBGTrigClasses.At(i))->String());
769       n++;
770     }
771
772     
773
774   }
775     
776   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
777   for (Int_t i=0; i<count; i++)
778   {
779     AliTriggerAnalysis* triggerAnalysis = static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i));
780   
781     switch (runNumber)
782     {
783       case 104315:
784       case 104316:
785       case 104320:
786       case 104321:
787       case 104439:
788         triggerAnalysis->SetV0TimeOffset(7.5);
789         break;
790       default:
791         triggerAnalysis->SetV0TimeOffset(0);
792     }
793   }
794     
795   fCurrentRun = runNumber;
796   
797   TH1::AddDirectory(oldStatus);
798   
799   return kTRUE;
800 }
801
802 TH2F * AliPhysicsSelection::BookHistStatistics(const char * tag) {
803
804     // add 6 rows to count for the estimate of good, accidentals and
805     // BG and the ratio of BG and accidentals to total +ratio goot to
806     // first col + 2 for error on good.
807
808   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
809 #ifdef VERBOSE_STAT
810   Int_t extrarows = fComputeBG ? 8 : 0;
811 #else
812   Int_t extrarows = fComputeBG ? 3 : 0;
813 #endif
814   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);
815
816   h->GetXaxis()->SetBinLabel(kStatTriggerClass,  "Trigger class");
817   h->GetXaxis()->SetBinLabel(kStatHWTrig,        "Hardware trigger");
818   h->GetXaxis()->SetBinLabel(kStatFO1,           "FO >= 1");
819   h->GetXaxis()->SetBinLabel(kStatFO2,           "FO >= 2");
820   h->GetXaxis()->SetBinLabel(kStatV0A,           "V0A");
821   h->GetXaxis()->SetBinLabel(kStatV0C,           "V0C");
822   h->GetXaxis()->SetBinLabel(kStatFMD,           "FMD");
823   h->GetXaxis()->SetBinLabel(kStatSSD1,          "SSD >= 2");
824   h->GetXaxis()->SetBinLabel(kStatV0ABG,         "V0A BG");
825   h->GetXaxis()->SetBinLabel(kStatV0CBG,         "V0C BG");
826   h->GetXaxis()->SetBinLabel(kStatMB1,           "(FO >= 1 | V0A | V0C) & !V0 BG");
827   h->GetXaxis()->SetBinLabel(kStatMB1Prime,      "(FO >= 2 | (FO >= 1 & (V0A | V0C)) | (V0A &v0C) ) & !V0 BG");
828   h->GetXaxis()->SetBinLabel(kStatFO1AndV0,      "FO >= 1 & (V0A | V0C) & !V0 BG");
829   h->GetXaxis()->SetBinLabel(kStatV0,            "V0A & V0C & !V0 BG & !BG ID");
830   h->GetXaxis()->SetBinLabel(kStatOffline,       "Offline Trigger");
831   h->GetXaxis()->SetBinLabel(kStatAny2Hits,      "2 Hits & !V0 BG");
832   h->GetXaxis()->SetBinLabel(kStatBG,            "Background identification");
833   h->GetXaxis()->SetBinLabel(kStatAccepted,      "Accepted");
834
835   Int_t n = 1;
836   for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
837     {
838       h->GetYaxis()->SetBinLabel(n, ((TObjString*) fCollTrigClasses.At(i))->String());
839       n++;
840     }
841   for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
842     {
843       h->GetYaxis()->SetBinLabel(n, ((TObjString*) fBGTrigClasses.At(i))->String());
844       n++;
845     }
846
847   if(fComputeBG) {
848     h->GetYaxis()->SetBinLabel(n++, "BG (A+C)");
849     h->GetYaxis()->SetBinLabel(n++, "ACC");
850 #ifdef VERBOSE_STAT
851     h->GetYaxis()->SetBinLabel(n++, "BG (A+C) %  (rel. to CINT1B)");
852     h->GetYaxis()->SetBinLabel(n++, "ACC % (rel. to CINT1B)");
853     h->GetYaxis()->SetBinLabel(n++, "ERR GOOD %");
854     h->GetYaxis()->SetBinLabel(n++, "GOOD % (rel. to 1st col)");
855     h->GetYaxis()->SetBinLabel(n++, "ERR GOOD");
856 #endif
857     h->GetYaxis()->SetBinLabel(n++, "GOOD");
858   }
859
860   return h;
861 }
862
863 void AliPhysicsSelection::Print(Option_t* /* option */) const
864 {
865   // print the configuration
866   
867   Printf("Configuration initialized for run %d (MC: %d):", fCurrentRun, fMC);
868   
869   Printf("Collision trigger classes:");
870   for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
871     Printf("%s", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
872   
873   Printf("Background trigger classes:");
874   for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
875     Printf("%s", ((TObjString*) fBGTrigClasses.At(i))->String().Data());
876
877   AliTriggerAnalysis* triggerAnalysis = dynamic_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(0));
878   
879   if (triggerAnalysis)
880   {
881     if (triggerAnalysis->GetV0TimeOffset() > 0)
882       Printf("V0 time offset active: %.2f ns", triggerAnalysis->GetV0TimeOffset());
883     
884     Printf("\nTotal available events:");
885     
886     triggerAnalysis->PrintTriggerClasses();
887   }
888   
889   if (fHistStatistics[kStatIdxAll])
890   {
891     for (Int_t i=0; i<fCollTrigClasses.GetEntries(); i++)
892     {
893       Printf("\nSelection statistics for collision trigger %s:", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
894       
895       Printf("Total events with correct trigger class: %d", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(1, i+1));
896       Printf("Selected collision candidates: %d", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(fHistStatistics[kStatIdxAll]->GetXaxis()->FindBin("Accepted"), i+1));
897     }
898   }
899   
900   if (fHistBunchCrossing)
901   {
902     Printf("\nBunch crossing statistics:");
903     
904     for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
905     {
906       TString str;
907       str.Form("Trigger %s has accepted events in the bunch crossings: ", fHistBunchCrossing->GetYaxis()->GetBinLabel(i));
908       
909       for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
910         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
911           str += Form("%d, ", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
912        
913       Printf("%s", str.Data());
914     }
915     
916     for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
917     {
918       Int_t countColl = 0;
919       Int_t countBG = 0;
920       for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
921       {
922         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
923         {
924           if (fCollTrigClasses.FindObject(fHistBunchCrossing->GetYaxis()->GetBinLabel(i)))
925             countColl++;
926           if (fBGTrigClasses.FindObject(fHistBunchCrossing->GetYaxis()->GetBinLabel(i)))
927             countBG++;
928         }
929       }
930       if (countColl > 0 && countBG > 0)
931         Printf("WARNING: Bunch crossing %d has collision and BG trigger classes active. Check BPTX functioning for this run!", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
932     }
933   }
934
935   if (fUsingCustomClasses)        
936     Printf("WARNING: Using custom trigger classes!");
937   if (fSkipTriggerClassSelection) 
938     Printf("WARNING: Skipping trigger class selection!");
939   if (fSkipV0) 
940     Printf("WARNING: Ignoring V0 information in selection");
941   if(!fBin0CallBack) 
942     Printf("WARNING: Callback not set: will not fill the statistics for the bin 0");
943
944 }
945
946 Long64_t AliPhysicsSelection::Merge(TCollection* list)
947 {
948   // Merge a list of AliMultiplicityCorrection objects with this (needed for
949   // PROOF).
950   // Returns the number of merged objects (including this).
951
952   if (!list)
953     return 0;
954
955   if (list->IsEmpty())
956     return 1;
957
958   TIterator* iter = list->MakeIterator();
959   TObject* obj;
960   
961   // collections of all histograms
962   const Int_t nHists = 9;
963   TList collections[nHists];
964
965   Int_t count = 0;
966   while ((obj = iter->Next())) {
967
968     AliPhysicsSelection* entry = dynamic_cast<AliPhysicsSelection*> (obj);
969     if (entry == 0) 
970       continue;
971       
972     collections[0].Add(&(entry->fTriggerAnalysis));
973     if (entry->fHistStatistics[0])
974       collections[1].Add(entry->fHistStatistics[0]);
975     if (entry->fHistStatistics[1])
976       collections[2].Add(entry->fHistStatistics[1]);
977     if (entry->fHistBunchCrossing)
978       collections[3].Add(entry->fHistBunchCrossing);
979     if (entry->fHistTriggerPattern)
980       collections[4].Add(entry->fHistTriggerPattern);
981     if (entry->fBackgroundIdentification)
982       collections[5].Add(entry->fBackgroundIdentification);
983
984     count++;
985   }
986
987   fTriggerAnalysis.Merge(&collections[0]);
988   if (fHistStatistics[0])
989     fHistStatistics[0]->Merge(&collections[1]);
990   if (fHistStatistics[1])
991     fHistStatistics[1]->Merge(&collections[2]);
992   if (fHistBunchCrossing)
993     fHistBunchCrossing->Merge(&collections[3]);
994   if (fHistTriggerPattern)
995     fHistTriggerPattern->Merge(&collections[4]);
996   if (fBackgroundIdentification)
997     fBackgroundIdentification->Merge(&collections[5]);
998   
999   delete iter;
1000
1001   return count+1;
1002 }
1003
1004 void AliPhysicsSelection::SaveHistograms(const char* folder) const
1005 {
1006   // write histograms to current directory
1007   
1008   if (!fHistStatistics[0] || !fHistStatistics[1])
1009     return;
1010     
1011   if (folder)
1012   {
1013     gDirectory->mkdir(folder);
1014     gDirectory->cd(folder);
1015   }
1016   
1017
1018   // Fill the last rows of fHistStatistics before saving
1019   if (fComputeBG) {
1020     Int_t triggerScheme = GetTriggerScheme(UInt_t(fCurrentRun));
1021     if(triggerScheme != 1){
1022       AliWarning("BG estimate only supported for trigger scheme \"1\" (CINT1 suite)");
1023     } else if (fUseMuonTriggers) {
1024       AliWarning("BG estimate with muon triggers to be implemented");
1025     } else {
1026       // 0. Determine the ratios of triggers E/B, A/B, C/B from the stat histogram
1027       // Those are used to rescale the different classes to the same number of bx ids
1028       Float_t nB = 0;
1029       Float_t nC = 0;
1030       Float_t nA = 0;
1031       Float_t nE = 0;
1032       for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++) {
1033         if (fHistBunchCrossing->GetBinContent(j, 1) > 0) nB++;
1034         if (fHistBunchCrossing->GetBinContent(j, 2) > 0) nA++;
1035         if (fHistBunchCrossing->GetBinContent(j, 3) > 0) nC++;
1036         if (fHistBunchCrossing->GetBinContent(j, 4) > 0) nE++;
1037       }
1038       Float_t ratioBE = nB/nE;
1039       Float_t ratioBA = nB/nA;
1040       Float_t ratioBC = nB/nC;
1041       Printf("Ratio between the BX ids in the different trigger classes:");
1042       Printf("  B/E = %f",ratioBE);
1043       Printf("  B/A = %f",ratioBA);
1044       Printf("  B/C = %f",ratioBC);
1045       Int_t nHistStat = 2;
1046       // TODO: get number of rows in a more flexible way
1047       // 1. loop over all cols
1048
1049       for(Int_t iHistStat = 0; iHistStat < nHistStat; iHistStat++){
1050     
1051         Int_t ncol = fHistStatistics[iHistStat]->GetNbinsX();
1052         Float_t good1 = 0;
1053         for(Int_t icol = 1; icol <= ncol; icol++) {
1054           Int_t cint1B = (Int_t) fHistStatistics[iHistStat]->GetBinContent(icol,1);     
1055           Int_t cint1A = (Int_t) fHistStatistics[iHistStat]->GetBinContent(icol,2);     
1056           Int_t cint1C = (Int_t) fHistStatistics[iHistStat]->GetBinContent(icol,3);     
1057           Int_t cint1E = (Int_t) fHistStatistics[iHistStat]->GetBinContent(icol,4);      
1058       
1059           if (cint1B>0) {
1060             Int_t acc  = ratioBE*cint1E; 
1061             Double_t acc_err = TMath::Sqrt(ratioBE*ratioBE*cint1E);
1062             //      Int_t bg   = cint1A + cint1C - 2*acc;
1063             Float_t bg   = fBIFactorA*(ratioBA*cint1A-acc) + fBIFactorC*(ratioBC*cint1C-acc) ;
1064             Float_t good = Float_t(cint1B) - bg - acc;
1065             if (icol ==1) good1 = good;
1066             //      Float_t errGood     = TMath::Sqrt(2*(cint1A+cint1C+cint1E));// Error on the number of goods assuming only bg fluctuates
1067             //      DeltaG^2 = B + FA^2 A + FC^2 C + Ratio^2 (FA+FC-1)^2 E.
1068             Float_t errGood     = TMath::Sqrt( cint1B + 
1069                                                fBIFactorA*fBIFactorA*ratioBA*ratioBA*cint1A +
1070                                                fBIFactorC*fBIFactorC*ratioBC*ratioBC*cint1C +
1071                                                ratioBE * ratioBE * 
1072                                                (fBIFactorA + fBIFactorC - 1)*(fBIFactorA + fBIFactorC - 1)*cint1E);
1073             Float_t errBG = TMath::Sqrt(fBIFactorA*fBIFactorA*ratioBA*ratioBA*cint1A+
1074                                         fBIFactorC*fBIFactorC*ratioBC*ratioBC*cint1C+
1075                                         ratioBE*ratioBE*(fBIFactorA+fBIFactorC)*(fBIFactorA+fBIFactorC)*cint1E);
1076         
1077         
1078             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowBG,bg);      
1079             fHistStatistics[iHistStat]->SetBinError  (icol,kStatRowBG,errBG);   
1080             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowAcc,acc);    
1081             fHistStatistics[iHistStat]->SetBinError  (icol,kStatRowAcc,acc_err);        
1082             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowGood,good);    
1083             fHistStatistics[iHistStat]->SetBinError  (icol,kStatRowGood,errGood);    
1084
1085 #ifdef VERBOSE_STAT
1086             //kStatRowBG=5,kStatRowAcc,kStatRowBGFrac,kStatRowAccFrac,kStatRowErrGoodFrac,kStatRowGoodFrac,kStatRowGood,kStatRowErrGood
1087             Float_t accFrac   = Float_t(acc) / cint1B  *100;
1088             Float_t errAccFrac= Float_t(acc_err) / cint1B  *100;
1089             Float_t bgFrac    = Float_t(bg)  / cint1B  *100;
1090             Float_t goodFrac  = Float_t(good)  / good1 *100;
1091             Float_t errGoodFrac = errGood/good1 * 100;
1092             Float_t errFracBG = bg > 0 ? TMath::Sqrt((errBG/bg)*(errBG/bg) + 1/cint1B)*bgFrac : 0;
1093             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowBGFrac,bgFrac);      
1094             fHistStatistics[iHistStat]->SetBinError  (icol,kStatRowBGFrac,errFracBG);   
1095             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowAccFrac,accFrac);    
1096             fHistStatistics[iHistStat]->SetBinError  (icol,kStatRowAccFrac,errAccFrac);    
1097             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowGoodFrac,goodFrac);    
1098             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowErrGoodFrac,errGoodFrac);    
1099             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowErrGood,errGood);    
1100 #endif
1101           }
1102         }
1103       }
1104     }  
1105   }
1106
1107   fHistStatistics[0]->Write();
1108   fHistStatistics[1]->Write();
1109   if(fHistBunchCrossing ) fHistBunchCrossing ->Write();
1110   if(fHistTriggerPattern) fHistTriggerPattern->Write();
1111   
1112   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
1113   for (Int_t i=0; i < count; i++)
1114   {
1115     TString triggerClass = "trigger_histograms_";
1116     if (i < fCollTrigClasses.GetEntries())
1117       triggerClass += ((TObjString*) fCollTrigClasses.At(i))->String();
1118     else
1119       triggerClass += ((TObjString*) fBGTrigClasses.At(i - fCollTrigClasses.GetEntries()))->String();
1120   
1121     gDirectory->mkdir(triggerClass);
1122     gDirectory->cd(triggerClass);
1123   
1124     static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i))->SaveHistograms();
1125     
1126     gDirectory->cd("..");
1127   }
1128  
1129   if (fBackgroundIdentification)
1130   {
1131     gDirectory->mkdir("background_identification");
1132     gDirectory->cd("background_identification");
1133       
1134     fBackgroundIdentification->GetOutput()->Write();
1135       
1136     gDirectory->cd("..");
1137   }
1138   
1139   if (folder)
1140     gDirectory->cd("..");
1141 }
1142
1143 void AliPhysicsSelection::SetBIFactors(Int_t run) {
1144
1145   switch(run) {
1146   case 104155:
1147     fBIFactorA = 0.961912722908;
1148     fBIFactorC = 1.04992336081;
1149     break;
1150   case 104157:
1151     fBIFactorA = 0.947312854998;
1152     fBIFactorC = 1.01599706417;
1153     break;
1154   case 104159:
1155     fBIFactorA = 0.93659320151;
1156     fBIFactorC = 0.98580804207;
1157     break;
1158   case 104160:
1159     fBIFactorA = 0.929664189926;
1160     fBIFactorC = 0.963467679851;
1161     break;
1162   case 104315:
1163     fBIFactorA = 1.08939104979;
1164     fBIFactorC = 0.931113921925;
1165     break;
1166   case 104316:
1167     fBIFactorA = 1.08351880974;
1168     fBIFactorC = 0.916068345845;
1169     break;
1170   case 104320:
1171     fBIFactorA = 1.07669281245;
1172     fBIFactorC = 0.876818744763;
1173     break;
1174   case 104321:
1175     fBIFactorA = 1.00971079602;
1176     fBIFactorC = 0.773781299076;
1177     break;
1178   case 104792:
1179     fBIFactorA = 0.787215863962;
1180     fBIFactorC = 0.778253173071;
1181     break;
1182   case 104793:
1183     fBIFactorA = 0.692211363661;
1184     fBIFactorC = 0.733152456667;
1185     break;
1186   case 104799:
1187     fBIFactorA = 1.04027825161;
1188     fBIFactorC = 1.00530825942;
1189     break;
1190   case 104800:
1191     fBIFactorA = 1.05309910671;
1192     fBIFactorC = 1.00376801855;
1193     break;
1194   case 104801:
1195     fBIFactorA = 1.0531231922;
1196     fBIFactorC = 0.992439666758;
1197     break;
1198   case 104802:
1199     fBIFactorA = 1.04191478134;
1200     fBIFactorC = 0.979368585208;
1201     break;
1202   case 104803:
1203     fBIFactorA = 1.03121314094;
1204     fBIFactorC = 0.973379962609;
1205     break;
1206   case 104824:
1207     fBIFactorA = 0.969945926722;
1208     fBIFactorC = 0.39549745806;
1209     break;
1210   case 104825:
1211     fBIFactorA = 0.968627213937;
1212     fBIFactorC = 0.310100412205;
1213     break;
1214   case 104841:
1215     fBIFactorA = 0.991601393212;
1216     fBIFactorC = 0.83762204722;
1217     break;
1218   case 104845:
1219     fBIFactorA = 0.98040863886;
1220     fBIFactorC = 0.694824205793;
1221     break;
1222   case 104867:
1223     fBIFactorA = 1.10646173412;
1224     fBIFactorC = 0.841407246916;
1225     break;
1226   case 104876:
1227     fBIFactorA = 1.12063452421;
1228     fBIFactorC = 0.78726542895;
1229     break;
1230   case 104890:
1231     fBIFactorA = 1.02346137453;
1232     fBIFactorC = 1.03355663595;
1233     break;
1234   case 104892:
1235     fBIFactorA = 1.05406025913;
1236     fBIFactorC = 1.00029166135;
1237     break;
1238   case 105143:
1239     fBIFactorA = 0.947343384349;
1240     fBIFactorC = 0.972637444408;
1241     break;
1242   case 105160:
1243     fBIFactorA = 0.908854622177;
1244     fBIFactorC = 0.958851103977;
1245     break; 
1246   case 105256:
1247     fBIFactorA = 0.810076150206;
1248     fBIFactorC = 0.884663561883;
1249     break;
1250   case 105257:
1251     fBIFactorA = 0.80974912303;
1252     fBIFactorC = 0.878859123479;
1253     break;
1254   case 105268:
1255     fBIFactorA = 0.809052110679;
1256     fBIFactorC = 0.87233890989;
1257     break;
1258   default:
1259     fBIFactorA = 1;
1260     fBIFactorC = 1;
1261   }
1262
1263
1264 }