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