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