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