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