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