changing offline trigger definitions once more:
[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, YY 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 = str2.Atoll();
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 &%u","+CINT1B-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1B-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
691         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CINT1A-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1A-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
692         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CINT1C-ABCE-NOPF-ALL",  GetBXIDs(runNumber,"CINT1C-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
693         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CINT1-E-NOPF-ALL",      GetBXIDs(runNumber,"CINT1-E-NOPF-ALL"), (UInt_t) AliVEvent::kMB)));
694
695         // Muon trigger have the same BXIDs of the corresponding CINT triggers
696         fCollTrigClasses.Add(new TObjString(Form("%s%s &%u","+CMUS1B-ABCE-NOPF-MUON",  GetBXIDs(runNumber,"CINT1B-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
697         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CMUS1A-ABCE-NOPF-MUON",  GetBXIDs(runNumber,"CINT1A-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
698         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CMUS1C-ABCE-NOPF-MUON",  GetBXIDs(runNumber,"CINT1C-ABCE-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
699         fBGTrigClasses.Add  (new TObjString(Form("%s%s &%u","+CMUS1-E-NOPF-MUON"    ,  GetBXIDs(runNumber,"CINT1-E-NOPF-ALL"), (UInt_t) AliVEvent::kMUON)));
700         
701         // triggers classes used from August 2010
702         // MB
703         fCollTrigClasses.Add(new TObjString(Form("+CINT1-B-NOPF-ALLNOTRD &%u", (UInt_t) AliVEvent::kMB)));
704         fBGTrigClasses.Add  (new TObjString(Form("+CINT1-AC-NOPF-ALLNOTRD &%u", (UInt_t) AliVEvent::kMB)));
705         fBGTrigClasses.Add  (new TObjString(Form("+CINT1-E-NOPF-ALLNOTRD &%u", (UInt_t) AliVEvent::kMB)));
706         
707         // MUON
708         fCollTrigClasses.Add(new TObjString(Form("+CMUS1-B-NOPF-ALLNOTRD &%u", (UInt_t) AliVEvent::kMUON)));
709         fBGTrigClasses.Add  (new TObjString(Form("+CMUS1-AC-NOPF-ALLNOTRD &%u", (UInt_t) AliVEvent::kMUON)));
710         fBGTrigClasses.Add  (new TObjString(Form("+CMUS1-E-NOPF-ALLNOTRD &%u", (UInt_t) AliVEvent::kMUON)));
711
712         // High Multiplicity
713         fCollTrigClasses.Add(new TObjString(Form("+CSH1-B-NOPF-ALLNOTRD &%u", (UInt_t) AliVEvent::kHighMult)));
714         fBGTrigClasses.Add  (new TObjString(Form("+CSH1-AC-NOPF-ALLNOTRD &%u", (UInt_t) AliVEvent::kHighMult)));
715         fBGTrigClasses.Add  (new TObjString(Form("+CSH1-E-NOPF-ALLNOTRD &%u", (UInt_t) AliVEvent::kHighMult)));
716
717         break;
718         
719       case 2:
720         fCollTrigClasses.Add(new TObjString(Form("+CSMBB-ABCE-NOPF-ALL &%u", (UInt_t) AliVEvent::kMB)));
721         fBGTrigClasses.Add(new TObjString(Form("+CSMBA-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL &%u", (UInt_t) AliVEvent::kMB)));
722         fBGTrigClasses.Add(new TObjString(Form("+CSMBC-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL &%u", (UInt_t) AliVEvent::kMB)));
723         break;
724         
725       case 3:
726         // 
727         break;
728         
729       default:
730         AliFatal(Form("Unsupported trigger scheme %d", triggerScheme));
731       }
732     }
733     
734     Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
735     
736     for (Int_t i=0; i<count; i++)
737     {
738       AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis;
739       triggerAnalysis->SetAnalyzeMC(fMC);
740       triggerAnalysis->EnableHistograms();
741       triggerAnalysis->SetSPDGFOThreshhold(1);
742       fTriggerAnalysis.Add(triggerAnalysis);
743     }
744       
745     // TODO: shall I really delete this?
746     if (fHistStatistics[0])
747       delete fHistStatistics[0];
748     if (fHistStatistics[1])
749       delete fHistStatistics[1];
750   
751     fHistStatistics[kStatIdxBin0] = BookHistStatistics("_Bin0");
752     fHistStatistics[kStatIdxAll]  = BookHistStatistics("");
753     
754     if (fHistBunchCrossing)
755       delete fHistBunchCrossing;
756   
757     fHistBunchCrossing = new TH2F("fHistBunchCrossing", "fHistBunchCrossing;bunch crossing number;", 4000, -0.5, 3999.5,  count, -0.5, -0.5 + count);
758
759     if (fHistTriggerPattern)
760       delete fHistTriggerPattern;
761     
762     const int ntrig=3;
763     Int_t n = 1;
764     const Int_t nbinTrig = TMath::Nint(TMath::Power(2,ntrig));
765
766     fHistTriggerPattern = new TH1F("fHistTriggerPattern", "Trigger pattern: FO + 2*v0A + 4*v0C", 
767                                    nbinTrig, -0.5, nbinTrig-0.5);    
768     fHistTriggerPattern->GetXaxis()->SetBinLabel(1,"NO TRIG");
769     fHistTriggerPattern->GetXaxis()->SetBinLabel(2,"FO");
770     fHistTriggerPattern->GetXaxis()->SetBinLabel(3,"v0A");
771     fHistTriggerPattern->GetXaxis()->SetBinLabel(4,"FO & v0A");
772     fHistTriggerPattern->GetXaxis()->SetBinLabel(5,"v0C");
773     fHistTriggerPattern->GetXaxis()->SetBinLabel(6,"FO & v0C");
774     fHistTriggerPattern->GetXaxis()->SetBinLabel(7,"v0A & v0C");
775     fHistTriggerPattern->GetXaxis()->SetBinLabel(8,"FO & v0A & v0C");
776
777   
778     n = 1;
779     for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
780     {
781       fHistBunchCrossing->GetYaxis()->SetBinLabel(n, ((TObjString*) fCollTrigClasses.At(i))->String());
782       n++;
783     }
784     for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
785     {
786       fHistBunchCrossing->GetYaxis()->SetBinLabel(n, ((TObjString*) fBGTrigClasses.At(i))->String());
787       n++;
788     }
789
790     
791
792   }
793     
794   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
795   for (Int_t i=0; i<count; i++)
796   {
797     AliTriggerAnalysis* triggerAnalysis = static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i));
798   
799     switch (runNumber)
800     {
801       case 104315:
802       case 104316:
803       case 104320:
804       case 104321:
805       case 104439:
806         triggerAnalysis->SetV0TimeOffset(7.5);
807         break;
808       default:
809         triggerAnalysis->SetV0TimeOffset(0);
810     }
811   }
812     
813   fCurrentRun = runNumber;
814   
815   TH1::AddDirectory(oldStatus);
816   
817   return kTRUE;
818 }
819
820 TH2F * AliPhysicsSelection::BookHistStatistics(const char * tag) {
821
822     // add 6 rows to count for the estimate of good, accidentals and
823     // BG and the ratio of BG and accidentals to total +ratio goot to
824     // first col + 2 for error on good.
825
826   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
827 #ifdef VERBOSE_STAT
828   Int_t extrarows = fComputeBG ? 8 : 0;
829 #else
830   Int_t extrarows = fComputeBG ? 3 : 0;
831 #endif
832   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);
833
834   h->GetXaxis()->SetBinLabel(kStatTriggerClass,  "Trigger class");
835   h->GetXaxis()->SetBinLabel(kStatHWTrig,        "Hardware trigger");
836   h->GetXaxis()->SetBinLabel(kStatFO1,           "FO >= 1");
837   h->GetXaxis()->SetBinLabel(kStatFO2,           "FO >= 2");
838   h->GetXaxis()->SetBinLabel(kStatV0A,           "V0A");
839   h->GetXaxis()->SetBinLabel(kStatV0C,           "V0C");
840   h->GetXaxis()->SetBinLabel(kStatFMD,           "FMD");
841   h->GetXaxis()->SetBinLabel(kStatSSD1,          "SSD >= 2");
842   h->GetXaxis()->SetBinLabel(kStatV0ABG,         "V0A BG");
843   h->GetXaxis()->SetBinLabel(kStatV0CBG,         "V0C BG");
844   h->GetXaxis()->SetBinLabel(kStatMB1,           "(FO >= 1 | V0A | V0C) & !V0 BG");
845   h->GetXaxis()->SetBinLabel(kStatMB1Prime,      "(FO >= 2 | (FO >= 1 & (V0A | V0C)) | (V0A &v0C) ) & !V0 BG");
846   h->GetXaxis()->SetBinLabel(kStatFO1AndV0,      "FO >= 1 & (V0A | V0C) & !V0 BG");
847   h->GetXaxis()->SetBinLabel(kStatV0,            "V0A & V0C & !V0 BG & !BG ID");
848   h->GetXaxis()->SetBinLabel(kStatOffline,       "Offline Trigger");
849   h->GetXaxis()->SetBinLabel(kStatAny2Hits,      "2 Hits & !V0 BG");
850   h->GetXaxis()->SetBinLabel(kStatBG,            "Background identification");
851   h->GetXaxis()->SetBinLabel(kStatAccepted,      "Accepted");
852
853   Int_t n = 1;
854   for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
855     {
856       h->GetYaxis()->SetBinLabel(n, ((TObjString*) fCollTrigClasses.At(i))->String());
857       n++;
858     }
859   for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
860     {
861       h->GetYaxis()->SetBinLabel(n, ((TObjString*) fBGTrigClasses.At(i))->String());
862       n++;
863     }
864
865   if(fComputeBG) {
866     h->GetYaxis()->SetBinLabel(n++, "BG (A+C)");
867     h->GetYaxis()->SetBinLabel(n++, "ACC");
868 #ifdef VERBOSE_STAT
869     h->GetYaxis()->SetBinLabel(n++, "BG (A+C) %  (rel. to CINT1B)");
870     h->GetYaxis()->SetBinLabel(n++, "ACC % (rel. to CINT1B)");
871     h->GetYaxis()->SetBinLabel(n++, "ERR GOOD %");
872     h->GetYaxis()->SetBinLabel(n++, "GOOD % (rel. to 1st col)");
873     h->GetYaxis()->SetBinLabel(n++, "ERR GOOD");
874 #endif
875     h->GetYaxis()->SetBinLabel(n++, "GOOD");
876   }
877
878   return h;
879 }
880
881 void AliPhysicsSelection::Print(Option_t* /* option */) const
882 {
883   // print the configuration
884   
885   Printf("Configuration initialized for run %d (MC: %d):", fCurrentRun, fMC);
886   
887   Printf("Collision trigger classes:");
888   for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
889     Printf("%s", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
890   
891   Printf("Background trigger classes:");
892   for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
893     Printf("%s", ((TObjString*) fBGTrigClasses.At(i))->String().Data());
894
895   AliTriggerAnalysis* triggerAnalysis = dynamic_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(0));
896   
897   if (triggerAnalysis)
898   {
899     if (triggerAnalysis->GetV0TimeOffset() > 0)
900       Printf("V0 time offset active: %.2f ns", triggerAnalysis->GetV0TimeOffset());
901     
902     Printf("\nTotal available events:");
903     
904     triggerAnalysis->PrintTriggerClasses();
905   }
906   
907   if (fHistStatistics[kStatIdxAll])
908   {
909     for (Int_t i=0; i<fCollTrigClasses.GetEntries(); i++)
910     {
911       Printf("\nSelection statistics for collision trigger %s:", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
912       
913       Printf("Total events with correct trigger class: %d", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(1, i+1));
914       Printf("Selected collision candidates: %d", (Int_t) fHistStatistics[kStatIdxAll]->GetBinContent(fHistStatistics[kStatIdxAll]->GetXaxis()->FindBin("Accepted"), i+1));
915     }
916   }
917   
918   if (fHistBunchCrossing)
919   {
920     Printf("\nBunch crossing statistics:");
921     
922     for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
923     {
924       TString str;
925       str.Form("Trigger %s has accepted events in the bunch crossings: ", fHistBunchCrossing->GetYaxis()->GetBinLabel(i));
926       
927       for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
928         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
929           str += Form("%d, ", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
930        
931       Printf("%s", str.Data());
932     }
933     
934     for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
935     {
936       Int_t countColl = 0;
937       Int_t countBG = 0;
938       for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
939       {
940         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
941         {
942           if (fCollTrigClasses.FindObject(fHistBunchCrossing->GetYaxis()->GetBinLabel(i)))
943             countColl++;
944           if (fBGTrigClasses.FindObject(fHistBunchCrossing->GetYaxis()->GetBinLabel(i)))
945             countBG++;
946         }
947       }
948       if (countColl > 0 && countBG > 0)
949         Printf("WARNING: Bunch crossing %d has collision and BG trigger classes active. Check BPTX functioning for this run!", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
950     }
951   }
952
953   if (fUsingCustomClasses)        
954     Printf("WARNING: Using custom trigger classes!");
955   if (fSkipTriggerClassSelection) 
956     Printf("WARNING: Skipping trigger class selection!");
957   if (fSkipV0) 
958     Printf("WARNING: Ignoring V0 information in selection");
959   if(!fBin0CallBack) 
960     Printf("WARNING: Callback not set: will not fill the statistics for the bin 0");
961
962 }
963
964 Long64_t AliPhysicsSelection::Merge(TCollection* list)
965 {
966   // Merge a list of AliMultiplicityCorrection objects with this (needed for
967   // PROOF).
968   // Returns the number of merged objects (including this).
969
970   if (!list)
971     return 0;
972
973   if (list->IsEmpty())
974     return 1;
975
976   TIterator* iter = list->MakeIterator();
977   TObject* obj;
978   
979   // collections of all histograms
980   const Int_t nHists = 9;
981   TList collections[nHists];
982
983   Int_t count = 0;
984   while ((obj = iter->Next())) {
985
986     AliPhysicsSelection* entry = dynamic_cast<AliPhysicsSelection*> (obj);
987     if (entry == 0) 
988       continue;
989       
990     collections[0].Add(&(entry->fTriggerAnalysis));
991     if (entry->fHistStatistics[0])
992       collections[1].Add(entry->fHistStatistics[0]);
993     if (entry->fHistStatistics[1])
994       collections[2].Add(entry->fHistStatistics[1]);
995     if (entry->fHistBunchCrossing)
996       collections[3].Add(entry->fHistBunchCrossing);
997     if (entry->fHistTriggerPattern)
998       collections[4].Add(entry->fHistTriggerPattern);
999     if (entry->fBackgroundIdentification)
1000       collections[5].Add(entry->fBackgroundIdentification);
1001
1002     count++;
1003   }
1004
1005   fTriggerAnalysis.Merge(&collections[0]);
1006   if (fHistStatistics[0])
1007     fHistStatistics[0]->Merge(&collections[1]);
1008   if (fHistStatistics[1])
1009     fHistStatistics[1]->Merge(&collections[2]);
1010   if (fHistBunchCrossing)
1011     fHistBunchCrossing->Merge(&collections[3]);
1012   if (fHistTriggerPattern)
1013     fHistTriggerPattern->Merge(&collections[4]);
1014   if (fBackgroundIdentification)
1015     fBackgroundIdentification->Merge(&collections[5]);
1016   
1017   delete iter;
1018
1019   return count+1;
1020 }
1021
1022 void AliPhysicsSelection::SaveHistograms(const char* folder) const
1023 {
1024   // write histograms to current directory
1025   
1026   if (!fHistStatistics[0] || !fHistStatistics[1])
1027     return;
1028     
1029   if (folder)
1030   {
1031     gDirectory->mkdir(folder);
1032     gDirectory->cd(folder);
1033   }
1034   
1035
1036   // Fill the last rows of fHistStatistics before saving
1037   if (fComputeBG) {
1038     Int_t triggerScheme = GetTriggerScheme(UInt_t(fCurrentRun));
1039     if(triggerScheme != 1){
1040       AliWarning("BG estimate only supported for trigger scheme \"1\" (CINT1 suite)");
1041     } else if (fUseMuonTriggers) {
1042       AliWarning("BG estimate with muon triggers to be implemented");
1043     } else {
1044       // 0. Determine the ratios of triggers E/B, A/B, C/B from the stat histogram
1045       // Those are used to rescale the different classes to the same number of bx ids
1046       Float_t nB = 0;
1047       Float_t nC = 0;
1048       Float_t nA = 0;
1049       Float_t nE = 0;
1050       for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++) {
1051         if (fHistBunchCrossing->GetBinContent(j, 1) > 0) nB++;
1052         if (fHistBunchCrossing->GetBinContent(j, 2) > 0) nA++;
1053         if (fHistBunchCrossing->GetBinContent(j, 3) > 0) nC++;
1054         if (fHistBunchCrossing->GetBinContent(j, 4) > 0) nE++;
1055       }
1056       Float_t ratioBE = nB/nE;
1057       Float_t ratioBA = nB/nA;
1058       Float_t ratioBC = nB/nC;
1059       Printf("Ratio between the BX ids in the different trigger classes:");
1060       Printf("  B/E = %f",ratioBE);
1061       Printf("  B/A = %f",ratioBA);
1062       Printf("  B/C = %f",ratioBC);
1063       Int_t nHistStat = 2;
1064       // TODO: get number of rows in a more flexible way
1065       // 1. loop over all cols
1066
1067       for(Int_t iHistStat = 0; iHistStat < nHistStat; iHistStat++){
1068     
1069         Int_t ncol = fHistStatistics[iHistStat]->GetNbinsX();
1070         Float_t good1 = 0;
1071         for(Int_t icol = 1; icol <= ncol; icol++) {
1072           Int_t cint1B = (Int_t) fHistStatistics[iHistStat]->GetBinContent(icol,1);     
1073           Int_t cint1A = (Int_t) fHistStatistics[iHistStat]->GetBinContent(icol,2);     
1074           Int_t cint1C = (Int_t) fHistStatistics[iHistStat]->GetBinContent(icol,3);     
1075           Int_t cint1E = (Int_t) fHistStatistics[iHistStat]->GetBinContent(icol,4);      
1076       
1077           if (cint1B>0) {
1078             Int_t acc  = ratioBE*cint1E; 
1079             Double_t acc_err = TMath::Sqrt(ratioBE*ratioBE*cint1E);
1080             //      Int_t bg   = cint1A + cint1C - 2*acc;
1081             Float_t bg   = fBIFactorA*(ratioBA*cint1A-acc) + fBIFactorC*(ratioBC*cint1C-acc) ;
1082             Float_t good = Float_t(cint1B) - bg - acc;
1083             if (icol ==1) good1 = good;
1084             //      Float_t errGood     = TMath::Sqrt(2*(cint1A+cint1C+cint1E));// Error on the number of goods assuming only bg fluctuates
1085             //      DeltaG^2 = B + FA^2 A + FC^2 C + Ratio^2 (FA+FC-1)^2 E.
1086             Float_t errGood     = TMath::Sqrt( cint1B + 
1087                                                fBIFactorA*fBIFactorA*ratioBA*ratioBA*cint1A +
1088                                                fBIFactorC*fBIFactorC*ratioBC*ratioBC*cint1C +
1089                                                ratioBE * ratioBE * 
1090                                                (fBIFactorA + fBIFactorC - 1)*(fBIFactorA + fBIFactorC - 1)*cint1E);
1091             Float_t errBG = TMath::Sqrt(fBIFactorA*fBIFactorA*ratioBA*ratioBA*cint1A+
1092                                         fBIFactorC*fBIFactorC*ratioBC*ratioBC*cint1C+
1093                                         ratioBE*ratioBE*(fBIFactorA+fBIFactorC)*(fBIFactorA+fBIFactorC)*cint1E);
1094         
1095         
1096             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowBG,bg);      
1097             fHistStatistics[iHistStat]->SetBinError  (icol,kStatRowBG,errBG);   
1098             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowAcc,acc);    
1099             fHistStatistics[iHistStat]->SetBinError  (icol,kStatRowAcc,acc_err);        
1100             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowGood,good);    
1101             fHistStatistics[iHistStat]->SetBinError  (icol,kStatRowGood,errGood);    
1102
1103 #ifdef VERBOSE_STAT
1104             //kStatRowBG=5,kStatRowAcc,kStatRowBGFrac,kStatRowAccFrac,kStatRowErrGoodFrac,kStatRowGoodFrac,kStatRowGood,kStatRowErrGood
1105             Float_t accFrac   = Float_t(acc) / cint1B  *100;
1106             Float_t errAccFrac= Float_t(acc_err) / cint1B  *100;
1107             Float_t bgFrac    = Float_t(bg)  / cint1B  *100;
1108             Float_t goodFrac  = Float_t(good)  / good1 *100;
1109             Float_t errGoodFrac = errGood/good1 * 100;
1110             Float_t errFracBG = bg > 0 ? TMath::Sqrt((errBG/bg)*(errBG/bg) + 1/cint1B)*bgFrac : 0;
1111             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowBGFrac,bgFrac);      
1112             fHistStatistics[iHistStat]->SetBinError  (icol,kStatRowBGFrac,errFracBG);   
1113             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowAccFrac,accFrac);    
1114             fHistStatistics[iHistStat]->SetBinError  (icol,kStatRowAccFrac,errAccFrac);    
1115             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowGoodFrac,goodFrac);    
1116             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowErrGoodFrac,errGoodFrac);    
1117             fHistStatistics[iHistStat]->SetBinContent(icol,kStatRowErrGood,errGood);    
1118 #endif
1119           }
1120         }
1121       }
1122     }  
1123   }
1124
1125   fHistStatistics[0]->Write();
1126   fHistStatistics[1]->Write();
1127   if(fHistBunchCrossing ) fHistBunchCrossing ->Write();
1128   if(fHistTriggerPattern) fHistTriggerPattern->Write();
1129   
1130   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
1131   for (Int_t i=0; i < count; i++)
1132   {
1133     TString triggerClass = "trigger_histograms_";
1134     if (i < fCollTrigClasses.GetEntries())
1135       triggerClass += ((TObjString*) fCollTrigClasses.At(i))->String();
1136     else
1137       triggerClass += ((TObjString*) fBGTrigClasses.At(i - fCollTrigClasses.GetEntries()))->String();
1138   
1139     gDirectory->mkdir(triggerClass);
1140     gDirectory->cd(triggerClass);
1141   
1142     static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i))->SaveHistograms();
1143     
1144     gDirectory->cd("..");
1145   }
1146  
1147   if (fBackgroundIdentification)
1148   {
1149     gDirectory->mkdir("background_identification");
1150     gDirectory->cd("background_identification");
1151       
1152     fBackgroundIdentification->GetOutput()->Write();
1153       
1154     gDirectory->cd("..");
1155   }
1156   
1157   if (folder)
1158     gDirectory->cd("..");
1159 }
1160
1161 void AliPhysicsSelection::SetBIFactors(Int_t run) {
1162
1163   switch(run) {
1164   case 104155:
1165     fBIFactorA = 0.961912722908;
1166     fBIFactorC = 1.04992336081;
1167     break;
1168   case 104157:
1169     fBIFactorA = 0.947312854998;
1170     fBIFactorC = 1.01599706417;
1171     break;
1172   case 104159:
1173     fBIFactorA = 0.93659320151;
1174     fBIFactorC = 0.98580804207;
1175     break;
1176   case 104160:
1177     fBIFactorA = 0.929664189926;
1178     fBIFactorC = 0.963467679851;
1179     break;
1180   case 104315:
1181     fBIFactorA = 1.08939104979;
1182     fBIFactorC = 0.931113921925;
1183     break;
1184   case 104316:
1185     fBIFactorA = 1.08351880974;
1186     fBIFactorC = 0.916068345845;
1187     break;
1188   case 104320:
1189     fBIFactorA = 1.07669281245;
1190     fBIFactorC = 0.876818744763;
1191     break;
1192   case 104321:
1193     fBIFactorA = 1.00971079602;
1194     fBIFactorC = 0.773781299076;
1195     break;
1196   case 104792:
1197     fBIFactorA = 0.787215863962;
1198     fBIFactorC = 0.778253173071;
1199     break;
1200   case 104793:
1201     fBIFactorA = 0.692211363661;
1202     fBIFactorC = 0.733152456667;
1203     break;
1204   case 104799:
1205     fBIFactorA = 1.04027825161;
1206     fBIFactorC = 1.00530825942;
1207     break;
1208   case 104800:
1209     fBIFactorA = 1.05309910671;
1210     fBIFactorC = 1.00376801855;
1211     break;
1212   case 104801:
1213     fBIFactorA = 1.0531231922;
1214     fBIFactorC = 0.992439666758;
1215     break;
1216   case 104802:
1217     fBIFactorA = 1.04191478134;
1218     fBIFactorC = 0.979368585208;
1219     break;
1220   case 104803:
1221     fBIFactorA = 1.03121314094;
1222     fBIFactorC = 0.973379962609;
1223     break;
1224   case 104824:
1225     fBIFactorA = 0.969945926722;
1226     fBIFactorC = 0.39549745806;
1227     break;
1228   case 104825:
1229     fBIFactorA = 0.968627213937;
1230     fBIFactorC = 0.310100412205;
1231     break;
1232   case 104841:
1233     fBIFactorA = 0.991601393212;
1234     fBIFactorC = 0.83762204722;
1235     break;
1236   case 104845:
1237     fBIFactorA = 0.98040863886;
1238     fBIFactorC = 0.694824205793;
1239     break;
1240   case 104867:
1241     fBIFactorA = 1.10646173412;
1242     fBIFactorC = 0.841407246916;
1243     break;
1244   case 104876:
1245     fBIFactorA = 1.12063452421;
1246     fBIFactorC = 0.78726542895;
1247     break;
1248   case 104890:
1249     fBIFactorA = 1.02346137453;
1250     fBIFactorC = 1.03355663595;
1251     break;
1252   case 104892:
1253     fBIFactorA = 1.05406025913;
1254     fBIFactorC = 1.00029166135;
1255     break;
1256   case 105143:
1257     fBIFactorA = 0.947343384349;
1258     fBIFactorC = 0.972637444408;
1259     break;
1260   case 105160:
1261     fBIFactorA = 0.908854622177;
1262     fBIFactorC = 0.958851103977;
1263     break; 
1264   case 105256:
1265     fBIFactorA = 0.810076150206;
1266     fBIFactorC = 0.884663561883;
1267     break;
1268   case 105257:
1269     fBIFactorA = 0.80974912303;
1270     fBIFactorC = 0.878859123479;
1271     break;
1272   case 105268:
1273     fBIFactorA = 0.809052110679;
1274     fBIFactorC = 0.87233890989;
1275     break;
1276   default:
1277     fBIFactorA = 1;
1278     fBIFactorC = 1;
1279   }
1280
1281
1282 }