e425a92b778ad2e0fffd6741c1d7c9a087c0bade
[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 //   Origin: Jan Fiete Grosse-Oetringhaus, CERN
42 //-------------------------------------------------------------------------
43
44 #include <Riostream.h>
45 #include <TH1F.h>
46 #include <TH2F.h>
47 #include <TList.h>
48 #include <TIterator.h>
49 #include <TDirectory.h>
50 #include <TObjArray.h>
51
52 #include <AliPhysicsSelection.h>
53
54 #include <AliTriggerAnalysis.h>
55 #include <AliLog.h>
56
57 #include <AliESDEvent.h>
58
59 ClassImp(AliPhysicsSelection)
60
61 AliPhysicsSelection::AliPhysicsSelection() :
62   AliAnalysisCuts("AliPhysicsSelection", "AliPhysicsSelection"),
63   fCurrentRun(-1),
64   fMC(kFALSE),
65   fCollTrigClasses(),
66   fBGTrigClasses(),
67   fTriggerAnalysis(),
68   fBackgroundIdentification(0),
69   fHistStatistics(0),
70   fHistBunchCrossing(0)
71 {
72   // constructor
73   
74   fCollTrigClasses.SetOwner(1);
75   fBGTrigClasses.SetOwner(1);
76   fTriggerAnalysis.SetOwner(1);
77   
78   AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kWarning);
79 }
80     
81 AliPhysicsSelection::~AliPhysicsSelection()
82 {
83   // destructor
84   
85   fCollTrigClasses.Delete();
86   fBGTrigClasses.Delete();
87   fTriggerAnalysis.Delete();
88
89   if (fHistStatistics)
90   {
91     delete fHistStatistics;
92     fHistStatistics = 0;
93   }
94   
95   if (fHistBunchCrossing)
96   {
97     delete fHistBunchCrossing;
98     fHistBunchCrossing = 0;
99   }
100 }
101
102 Bool_t AliPhysicsSelection::CheckTriggerClass(const AliESDEvent* aEsd, const char* trigger) const
103 {
104   // checks if the given trigger class(es) are found for the current event
105   // format of trigger: +TRIGGER1 -TRIGGER2
106   //   requires TRIGGER1 and rejects TRIGGER2
107   
108   TString str(trigger);
109   TObjArray* tokens = str.Tokenize(" ");
110   
111   for (Int_t i=0; i < tokens->GetEntries(); i++)
112   {
113     TString str2(((TObjString*) tokens->At(i))->String());
114     
115     if (str2[0] != '+' && str2[0] != '-')
116       AliFatal(Form("Invalid trigger syntax: %s", trigger));
117       
118     Bool_t flag = (str2[0] == '+');
119     
120     str2.Remove(0, 1);
121     
122     if (flag && !aEsd->IsTriggerClassFired(str2))
123     {
124       AliDebug(AliLog::kDebug, Form("Rejecting event because trigger class %s is not present", str2.Data()));
125       delete tokens;
126       return kFALSE;
127     }
128     if (!flag && aEsd->IsTriggerClassFired(str2))
129     {
130       AliDebug(AliLog::kDebug, Form("Rejecting event because trigger class %s is present", str2.Data()));
131       delete tokens;
132       return kFALSE;
133     }
134   }
135   
136   delete tokens;
137   return kTRUE;
138 }
139     
140 Bool_t AliPhysicsSelection::IsCollisionCandidate(const AliESDEvent* aEsd)
141 {
142   // checks if the given event is a collision candidate
143   
144   if (fCurrentRun != aEsd->GetRunNumber())
145     if (!Initialize(aEsd->GetRunNumber()))
146       AliFatal(Form("Could not initialize for run %d", aEsd->GetRunNumber()));
147     
148   const AliESDHeader* esdHeader = aEsd->GetHeader();
149   if (!esdHeader)
150   {
151     AliError("ESD Header could not be retrieved");
152     return kFALSE;
153   }
154   
155   // check event type; should be PHYSICS = 7 for data and 0 for MC
156   if (!fMC)
157   {
158     if (esdHeader->GetEventType() != 7)
159       return kFALSE;
160   }
161   else
162   {
163     if (esdHeader->GetEventType() != 0)
164       AliFatal(Form("Invalid event type for MC: %d", esdHeader->GetEventType()));
165   }
166   
167   Bool_t accept = kFALSE;
168     
169   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
170   for (Int_t i=0; i < count; i++)
171   {
172     const char* triggerClass = 0;
173     if (i < fCollTrigClasses.GetEntries())
174       triggerClass = ((TObjString*) fCollTrigClasses.At(i))->String();
175     else
176       triggerClass = ((TObjString*) fBGTrigClasses.At(i - fCollTrigClasses.GetEntries()))->String();
177   
178     AliDebug(AliLog::kDebug, Form("Processing trigger class %s", triggerClass));
179   
180     AliTriggerAnalysis* triggerAnalysis = static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i));
181   
182     triggerAnalysis->FillTriggerClasses(aEsd);
183     
184     if (CheckTriggerClass(aEsd, triggerClass))
185     {
186       triggerAnalysis->FillHistograms(aEsd);
187   
188       fHistStatistics->Fill(1, i);
189       
190       // hardware trigger (should only remove events for MC)
191       // replay CINT1B hardware trigger
192       // TODO this has to depend on the actual hardware trigger (and that depends on the run...)
193       Int_t fastORHW = triggerAnalysis->SPDFiredChips(aEsd, 1); // SPD number of chips from trigger bits (!)
194       Bool_t v0A = triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0A);
195       Bool_t v0C = triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0C);
196         
197       if (fastORHW == 0 && !v0A && !v0C)
198       {
199         AliDebug(AliLog::kDebug, "Rejecting event because hardware trigger is not fired");
200         continue;
201       }
202       
203       fHistStatistics->Fill(2, i);
204     
205       // offline trigger
206       Int_t fastOROffline = triggerAnalysis->SPDFiredChips(aEsd, 0); // SPD number of chips from clusters (!)
207       Bool_t v0ABG = triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0ABG);
208       Bool_t v0CBG = triggerAnalysis->IsOfflineTriggerFired(aEsd, AliTriggerAnalysis::kV0CBG);
209       Bool_t v0BG = v0ABG || v0CBG;
210     
211       if (fastOROffline > 0)
212         fHistStatistics->Fill(3, i);
213       if (fastOROffline > 1)
214         fHistStatistics->Fill(4, i);
215         
216       if (v0A)
217         fHistStatistics->Fill(5, i);
218       if (v0C)
219         fHistStatistics->Fill(6, i);
220       if (v0ABG)
221         fHistStatistics->Fill(7, i);
222       if (v0CBG)
223         fHistStatistics->Fill(8, i);
224         
225       if (fastOROffline > 1 && !v0BG)
226         fHistStatistics->Fill(9, i);
227         
228       if ((fastOROffline > 0 || v0A || v0C) && !v0BG)
229         fHistStatistics->Fill(10, i);
230     
231       if (fastOROffline > 0 && (v0A || v0C) && !v0BG)
232         fHistStatistics->Fill(11, i);
233   
234       if (v0A && v0C && !v0BG)
235         fHistStatistics->Fill(12, i);
236         
237       if (fastOROffline > 1 || (fastOROffline > 0 && (v0A || v0C)) || (v0A && v0C))
238       {
239         if (!v0BG)
240         {
241           fHistStatistics->Fill(13, i);
242       
243           if (fBackgroundIdentification && !fBackgroundIdentification->IsSelected(const_cast<AliESDEvent*> (aEsd)))
244           {
245             AliDebug(AliLog::kDebug, "Rejecting event because of background identification");
246             fHistStatistics->Fill(14, i);
247           }
248           else
249           {
250             AliDebug(AliLog::kDebug, "Accepted event for histograms");
251             
252             fHistStatistics->Fill(15, i);
253             fHistBunchCrossing->Fill(aEsd->GetBunchCrossNumber(), i);
254             if (i < fCollTrigClasses.GetEntries())
255               accept = kTRUE;
256           }
257         }
258         else
259           AliDebug(AliLog::kDebug, "Rejecting event because of V0 BG flag");
260       }
261       else
262         AliDebug(AliLog::kDebug, "Rejecting event because trigger condition is not fulfilled");
263     }
264   }
265  
266   if (accept)
267     AliDebug(AliLog::kDebug, "Accepted event as collision candidate");
268   
269   return accept;
270 }
271
272 Int_t AliPhysicsSelection::GetTriggerScheme(UInt_t runNumber)
273 {
274   // returns the current trigger scheme (classes that are accepted/rejected)
275   
276   if (fMC)
277     return 0;
278     
279   // TODO dependent on run number
280   
281   switch (runNumber)
282   {
283     // CSMBB triggers
284     case 104044:
285     case 105054:
286     case 105057:
287       return 2;
288   }
289   
290   // default: CINT1 suite
291   return 1;
292 }  
293     
294 Bool_t AliPhysicsSelection::Initialize(UInt_t runNumber)
295 {
296   // initializes the object for the given run
297   // 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
298   
299   Bool_t oldStatus = TH1::AddDirectoryStatus();
300   TH1::AddDirectory(kFALSE);
301   
302   Int_t triggerScheme = GetTriggerScheme(runNumber);
303   
304   if (fCurrentRun != -1 && triggerScheme != GetTriggerScheme(fCurrentRun))
305     AliFatal("Processing several runs with different trigger schemes is not supported");
306   
307   AliInfo(Form("Initializing for run %d", runNumber));
308   
309   // initialize first time?
310   if (fCurrentRun == -1)
311   {
312     switch (triggerScheme)
313     {
314       case 0:
315         fCollTrigClasses.Add(new TObjString(""));
316         break;
317       
318       case 1:
319         fCollTrigClasses.Add(new TObjString("+CINT1B-ABCE-NOPF-ALL"));
320         fBGTrigClasses.Add(new TObjString("+CINT1A-ABCE-NOPF-ALL"));
321         fBGTrigClasses.Add(new TObjString("+CINT1C-ABCE-NOPF-ALL"));
322         fBGTrigClasses.Add(new TObjString("+CINT1-E-NOPF-ALL"));
323         break;
324         
325       case 2:
326         fCollTrigClasses.Add(new TObjString("+CSMBB-ABCE-NOPF-ALL"));
327         fBGTrigClasses.Add(new TObjString("+CSMBA-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL"));
328         fBGTrigClasses.Add(new TObjString("+CSMBC-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL"));
329         break;
330   
331       default:
332         AliFatal(Form("Unsupported trigger scheme %d", triggerScheme));
333     }
334   
335     Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
336     
337     for (Int_t i=0; i<count; i++)
338     {
339       AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis;
340       triggerAnalysis->SetAnalyzeMC(fMC);
341       triggerAnalysis->EnableHistograms();
342       triggerAnalysis->SetSPDGFOThreshhold(1);
343       fTriggerAnalysis.Add(triggerAnalysis);
344     }
345       
346     if (fHistStatistics)
347       delete fHistStatistics;
348   
349     fHistStatistics = new TH2F("fHistStatistics", "fHistStatistics;;", 15, 0.5, 15.5, count, -0.5, -0.5 + count);
350       
351     Int_t n = 1;
352     fHistStatistics->GetXaxis()->SetBinLabel(n++, "Trigger class");
353     fHistStatistics->GetXaxis()->SetBinLabel(n++, "Hardware trigger");
354     fHistStatistics->GetXaxis()->SetBinLabel(n++, "FO >= 1");
355     fHistStatistics->GetXaxis()->SetBinLabel(n++, "FO >= 2");
356     fHistStatistics->GetXaxis()->SetBinLabel(n++, "V0A");
357     fHistStatistics->GetXaxis()->SetBinLabel(n++, "V0C");
358     fHistStatistics->GetXaxis()->SetBinLabel(n++, "V0A BG");
359     fHistStatistics->GetXaxis()->SetBinLabel(n++, "V0C BG");
360     fHistStatistics->GetXaxis()->SetBinLabel(n++, "FO >= 2 &!V0 BG");
361     fHistStatistics->GetXaxis()->SetBinLabel(n++, "(FO >= 1 | V0A | V0C) & !V0 BG");
362     fHistStatistics->GetXaxis()->SetBinLabel(n++, "FO >= 1 & (V0A | V0C) & !V0 BG");
363     fHistStatistics->GetXaxis()->SetBinLabel(n++, "V0A & V0C & !V0 BG");
364     fHistStatistics->GetXaxis()->SetBinLabel(n++, "(FO >= 2 | (FO >= 1 & (V0A | V0C)) | (V0A & V0C)) & !V0 BG");
365     fHistStatistics->GetXaxis()->SetBinLabel(n++, "Background identification");
366     fHistStatistics->GetXaxis()->SetBinLabel(n++, "Accepted");
367     
368     if (fHistBunchCrossing)
369       delete fHistBunchCrossing;
370   
371     fHistBunchCrossing = new TH2F("fHistBunchCrossing", "fHistBunchCrossing;bunch crossing number;", 4000, -0.5, 3999.5,  count, -0.5, -0.5 + count);
372       
373     n = 1;
374     for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
375     {
376       fHistStatistics->GetYaxis()->SetBinLabel(n, ((TObjString*) fCollTrigClasses.At(i))->String());
377       fHistBunchCrossing->GetYaxis()->SetBinLabel(n, ((TObjString*) fCollTrigClasses.At(i))->String());
378       n++;
379     }
380     for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
381     {
382       fHistStatistics->GetYaxis()->SetBinLabel(n, ((TObjString*) fBGTrigClasses.At(i))->String());
383       fHistBunchCrossing->GetYaxis()->SetBinLabel(n, ((TObjString*) fBGTrigClasses.At(i))->String());
384       n++;
385     }
386   }
387     
388   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
389   for (Int_t i=0; i<count; i++)
390   {
391     AliTriggerAnalysis* triggerAnalysis = static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i));
392   
393     switch (runNumber)
394     {
395       case 104315:
396       case 104316:
397       case 104320:
398       case 104321:
399       case 104439:
400         triggerAnalysis->SetV0TimeOffset(7.5);
401         break;
402       default:
403         triggerAnalysis->SetV0TimeOffset(0);
404     }
405   }
406     
407   fCurrentRun = runNumber;
408   
409   TH1::AddDirectory(oldStatus);
410   
411   return kTRUE;
412 }
413
414 void AliPhysicsSelection::Print(Option_t* /* option */) const
415 {
416   // print the configuration
417   
418   Printf("Configuration initialized for run %d (MC: %d):", fCurrentRun, fMC);
419   
420   Printf("Collision trigger classes:");
421   for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
422     Printf("%s", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
423   
424   Printf("Background trigger classes:");
425   for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
426     Printf("%s", ((TObjString*) fBGTrigClasses.At(i))->String().Data());
427
428   AliTriggerAnalysis* triggerAnalysis = dynamic_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(0));
429   
430   if (triggerAnalysis)
431   {
432     if (triggerAnalysis->GetV0TimeOffset() > 0)
433       Printf("V0 time offset active: %.2f ns", triggerAnalysis->GetV0TimeOffset());
434     
435     Printf("\nTotal available events:");
436     
437     triggerAnalysis->PrintTriggerClasses();
438   }
439   
440   if (fHistStatistics && fCollTrigClasses.GetEntries() > 0)
441   {
442     Printf("\nSelection statistics for first collision trigger (%s):", ((TObjString*) fCollTrigClasses.First())->String().Data());
443     
444     Printf("Total events with correct trigger class: %d", (Int_t) fHistStatistics->GetBinContent(1, 1));
445     Printf("Selected collision candidates: %d", (Int_t) fHistStatistics->GetBinContent(fHistStatistics->GetXaxis()->FindBin("Accepted"), 1));
446   }
447   
448   if (fHistBunchCrossing)
449   {
450     Printf("\nBunch crossing statistics:");
451     
452     for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
453     {
454       TString str;
455       str.Form("Trigger %s has accepted events in the bunch crossings: ", fHistBunchCrossing->GetYaxis()->GetBinLabel(i));
456       
457       for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
458         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
459           str += Form("%d, ", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
460        
461       Printf("%s", str.Data());
462     }
463     
464     for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
465     {
466       Int_t count = 0;
467       for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
468       {
469         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
470           count++;
471       }
472       if (count > 1)
473         Printf("WARNING: Bunch crossing %d has more than one trigger class active. Check BPTX functioning for this run!", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
474     }
475   }
476 }
477
478 Long64_t AliPhysicsSelection::Merge(TCollection* list)
479 {
480   // Merge a list of AliMultiplicityCorrection objects with this (needed for
481   // PROOF).
482   // Returns the number of merged objects (including this).
483
484   if (!list)
485     return 0;
486
487   if (list->IsEmpty())
488     return 1;
489
490   TIterator* iter = list->MakeIterator();
491   TObject* obj;
492   
493   // collections of all histograms
494   const Int_t nHists = 9;
495   TList collections[nHists];
496
497   Int_t count = 0;
498   while ((obj = iter->Next())) {
499
500     AliPhysicsSelection* entry = dynamic_cast<AliPhysicsSelection*> (obj);
501     if (entry == 0) 
502       continue;
503       
504     collections[0].Add(&(entry->fTriggerAnalysis));
505     if (entry->fHistStatistics)
506       collections[1].Add(entry->fHistStatistics);
507     if (entry->fHistBunchCrossing)
508       collections[2].Add(entry->fHistBunchCrossing);
509     if (entry->fBackgroundIdentification)
510       collections[3].Add(entry->fBackgroundIdentification);
511
512     count++;
513   }
514
515   fTriggerAnalysis.Merge(&collections[0]);
516   if (fHistStatistics)
517     fHistStatistics->Merge(&collections[1]);
518   if (fHistBunchCrossing)
519     fHistBunchCrossing->Merge(&collections[2]);
520   if (fBackgroundIdentification)
521     fBackgroundIdentification->Merge(&collections[3]);
522   
523   delete iter;
524
525   return count+1;
526 }
527
528 void AliPhysicsSelection::SaveHistograms(const char* folder) const
529 {
530   // write histograms to current directory
531   
532   if (!fHistStatistics)
533     return;
534     
535   if (folder)
536   {
537     gDirectory->mkdir(folder);
538     gDirectory->cd(folder);
539   }
540   
541   fHistStatistics->Write();
542   fHistBunchCrossing->Write();
543   
544   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
545   for (Int_t i=0; i < count; i++)
546   {
547     TString triggerClass = "trigger_histograms_";
548     if (i < fCollTrigClasses.GetEntries())
549       triggerClass += ((TObjString*) fCollTrigClasses.At(i))->String();
550     else
551       triggerClass += ((TObjString*) fBGTrigClasses.At(i - fCollTrigClasses.GetEntries()))->String();
552   
553     gDirectory->mkdir(triggerClass);
554     gDirectory->cd(triggerClass);
555   
556     static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i))->SaveHistograms();
557     
558     gDirectory->cd("..");
559   }
560  
561   if (fBackgroundIdentification)
562   {
563     gDirectory->mkdir("background_identification");
564     gDirectory->cd("background_identification");
565       
566     fBackgroundIdentification->GetOutput()->Write();
567       
568     gDirectory->cd("..");
569   }
570   
571   if (folder)
572     gDirectory->cd("..");
573 }