ec2bbc5d106177f117098523e9a9e37bb4e0470e
[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   return 1;
281 }  
282     
283 Bool_t AliPhysicsSelection::Initialize(UInt_t runNumber)
284 {
285   // initializes the object for the given run
286   // 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
287   
288   Bool_t oldStatus = TH1::AddDirectoryStatus();
289   TH1::AddDirectory(kFALSE);
290   
291   Int_t triggerScheme = GetTriggerScheme(runNumber);
292   
293   if (fCurrentRun != -1 && triggerScheme != GetTriggerScheme(fCurrentRun))
294     AliFatal("Processing several runs with different trigger schemes is not supported");
295   
296   AliInfo(Form("Initializing for run %d", runNumber));
297   
298   Bool_t firstTime = kFALSE;
299   if (fCurrentRun == -1)
300   {
301     // initialize first time
302     fCurrentRun = runNumber;
303     firstTime = kTRUE;
304   }
305   
306   if (firstTime)
307   {
308     switch (triggerScheme)
309     {
310       case 0:
311         fCollTrigClasses.Add(new TObjString(""));
312         break;
313       
314       case 1:
315         fCollTrigClasses.Add(new TObjString("+CINT1B-ABCE-NOPF-ALL"));
316         fBGTrigClasses.Add(new TObjString("+CINT1A-ABCE-NOPF-ALL"));
317         fBGTrigClasses.Add(new TObjString("+CINT1C-ABCE-NOPF-ALL"));
318         fBGTrigClasses.Add(new TObjString("+CINT1-E-NOPF-ALL"));
319         break;
320         
321       default:
322         AliFatal(Form("Unsupported trigger scheme %d", triggerScheme));
323     }
324   
325     Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
326     
327     for (Int_t i=0; i<count; i++)
328     {
329       AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis;
330       triggerAnalysis->SetAnalyzeMC(fMC);
331       triggerAnalysis->EnableHistograms();
332       triggerAnalysis->SetSPDGFOThreshhold(1);
333       fTriggerAnalysis.Add(triggerAnalysis);
334     }
335       
336     if (fHistStatistics)
337       delete fHistStatistics;
338   
339     fHistStatistics = new TH2F("fHistStatistics", "fHistStatistics;;", 15, 0.5, 15.5, count, -0.5, -0.5 + count);
340       
341     Int_t n = 1;
342     fHistStatistics->GetXaxis()->SetBinLabel(n++, "Trigger class");
343     fHistStatistics->GetXaxis()->SetBinLabel(n++, "Hardware trigger");
344     fHistStatistics->GetXaxis()->SetBinLabel(n++, "FO >= 1");
345     fHistStatistics->GetXaxis()->SetBinLabel(n++, "FO >= 2");
346     fHistStatistics->GetXaxis()->SetBinLabel(n++, "V0A");
347     fHistStatistics->GetXaxis()->SetBinLabel(n++, "V0C");
348     fHistStatistics->GetXaxis()->SetBinLabel(n++, "V0A BG");
349     fHistStatistics->GetXaxis()->SetBinLabel(n++, "V0C BG");
350     fHistStatistics->GetXaxis()->SetBinLabel(n++, "FO >= 2 &!V0 BG");
351     fHistStatistics->GetXaxis()->SetBinLabel(n++, "(FO >= 1 | V0A | V0C) & !V0 BG");
352     fHistStatistics->GetXaxis()->SetBinLabel(n++, "FO >= 1 & (V0A | V0C) & !V0 BG");
353     fHistStatistics->GetXaxis()->SetBinLabel(n++, "V0A & V0C & !V0 BG");
354     fHistStatistics->GetXaxis()->SetBinLabel(n++, "(FO >= 2 | (FO >= 1 & (V0A | V0C)) | (V0A & V0C)) & !V0 BG");
355     fHistStatistics->GetXaxis()->SetBinLabel(n++, "Background identification");
356     fHistStatistics->GetXaxis()->SetBinLabel(n++, "Accepted");
357     
358     if (fHistBunchCrossing)
359       delete fHistBunchCrossing;
360   
361     fHistBunchCrossing = new TH2F("fHistBunchCrossing", "fHistBunchCrossing;bunch crossing number;", 4000, -0.5, 3999.5,  count, -0.5, -0.5 + count);
362       
363     n = 1;
364     for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
365     {
366       fHistStatistics->GetYaxis()->SetBinLabel(n, ((TObjString*) fCollTrigClasses.At(i))->String());
367       fHistBunchCrossing->GetYaxis()->SetBinLabel(n, ((TObjString*) fCollTrigClasses.At(i))->String());
368       n++;
369     }
370     for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
371     {
372       fHistStatistics->GetYaxis()->SetBinLabel(n, ((TObjString*) fBGTrigClasses.At(i))->String());
373       fHistBunchCrossing->GetYaxis()->SetBinLabel(n, ((TObjString*) fBGTrigClasses.At(i))->String());
374       n++;
375     }
376   }
377     
378   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
379   for (Int_t i=0; i<count; i++)
380   {
381     AliTriggerAnalysis* triggerAnalysis = static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i));
382   
383     switch (runNumber)
384     {
385       case 104315:
386       case 104316:
387       case 104320:
388       case 104321:
389       case 104439:
390         triggerAnalysis->SetV0TimeOffset(7.5);
391         break;
392       default:
393         triggerAnalysis->SetV0TimeOffset(0);
394     }
395   }
396     
397   TH1::AddDirectory(oldStatus);
398   
399   return kTRUE;
400 }
401
402 void AliPhysicsSelection::Print(Option_t* /* option */) const
403 {
404   // print the configuration
405   
406   Printf("Configuration initialized for run %d (MC: %d):", fCurrentRun, fMC);
407   
408   Printf("Collision trigger classes:");
409   for (Int_t i=0; i < fCollTrigClasses.GetEntries(); i++)
410     Printf("%s", ((TObjString*) fCollTrigClasses.At(i))->String().Data());
411   
412   Printf("Background trigger classes:");
413   for (Int_t i=0; i < fBGTrigClasses.GetEntries(); i++)
414     Printf("%s", ((TObjString*) fBGTrigClasses.At(i))->String().Data());
415
416   AliTriggerAnalysis* triggerAnalysis = dynamic_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(0));
417   
418   if (triggerAnalysis)
419   {
420     if (triggerAnalysis->GetV0TimeOffset() > 0)
421       Printf("V0 time offset active: %.2f ns", triggerAnalysis->GetV0TimeOffset());
422     
423     Printf("\nTotal available events:");
424     
425     triggerAnalysis->PrintTriggerClasses();
426   }
427   
428   if (fHistStatistics && fCollTrigClasses.GetEntries() > 0)
429   {
430     Printf("\nSelection statistics for first collision trigger (%s):", ((TObjString*) fCollTrigClasses.First())->String().Data());
431     
432     Printf("Total events with correct trigger class: %d", (Int_t) fHistStatistics->GetBinContent(1, 1));
433     Printf("Selected collision candidates: %d", (Int_t) fHistStatistics->GetBinContent(fHistStatistics->GetXaxis()->FindBin("Accepted"), 1));
434   }
435   
436   if (fHistBunchCrossing)
437   {
438     Printf("\nBunch crossing statistics:");
439     
440     for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
441     {
442       TString str;
443       str.Form("Trigger %s has accepted events in the bunch crossings: ", fHistBunchCrossing->GetYaxis()->GetBinLabel(i));
444       
445       for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
446         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
447           str += Form("%d, ", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
448        
449       Printf("%s", str.Data());
450     }
451     
452     for (Int_t j=1; j<=fHistBunchCrossing->GetNbinsX(); j++)
453     {
454       Int_t count = 0;
455       for (Int_t i=1; i<=fHistBunchCrossing->GetNbinsY(); i++)
456       {
457         if (fHistBunchCrossing->GetBinContent(j, i) > 0)
458           count++;
459       }
460       if (count > 1)
461         Printf("WARNING: Bunch crossing %d has more than one trigger class active. Check BPTX functioning for this run!", (Int_t) fHistBunchCrossing->GetXaxis()->GetBinCenter(j));
462     }
463   }
464 }
465
466 Long64_t AliPhysicsSelection::Merge(TCollection* list)
467 {
468   // Merge a list of AliMultiplicityCorrection objects with this (needed for
469   // PROOF).
470   // Returns the number of merged objects (including this).
471
472   if (!list)
473     return 0;
474
475   if (list->IsEmpty())
476     return 1;
477
478   TIterator* iter = list->MakeIterator();
479   TObject* obj;
480   
481   // collections of all histograms
482   const Int_t nHists = 9;
483   TList collections[nHists];
484
485   Int_t count = 0;
486   while ((obj = iter->Next())) {
487
488     AliPhysicsSelection* entry = dynamic_cast<AliPhysicsSelection*> (obj);
489     if (entry == 0) 
490       continue;
491       
492     collections[0].Add(&(entry->fTriggerAnalysis));
493     if (entry->fHistStatistics)
494       collections[1].Add(entry->fHistStatistics);
495     if (entry->fHistBunchCrossing)
496       collections[2].Add(entry->fHistBunchCrossing);
497     if (entry->fBackgroundIdentification)
498       collections[3].Add(entry->fBackgroundIdentification);
499
500     count++;
501   }
502
503   fTriggerAnalysis.Merge(&collections[0]);
504   if (fHistStatistics)
505     fHistStatistics->Merge(&collections[1]);
506   if (fHistBunchCrossing)
507     fHistBunchCrossing->Merge(&collections[2]);
508   if (fBackgroundIdentification)
509     fBackgroundIdentification->Merge(&collections[3]);
510   
511   delete iter;
512
513   return count+1;
514 }
515
516 void AliPhysicsSelection::SaveHistograms(const char* folder) const
517 {
518   // write histograms to current directory
519   
520   if (!fHistStatistics)
521     return;
522     
523   if (folder)
524   {
525     gDirectory->mkdir(folder);
526     gDirectory->cd(folder);
527   }
528   
529   fHistStatistics->Write();
530   fHistBunchCrossing->Write();
531   
532   Int_t count = fCollTrigClasses.GetEntries() + fBGTrigClasses.GetEntries();
533   for (Int_t i=0; i < count; i++)
534   {
535     TString triggerClass = "trigger_histograms_";
536     if (i < fCollTrigClasses.GetEntries())
537       triggerClass += ((TObjString*) fCollTrigClasses.At(i))->String();
538     else
539       triggerClass += ((TObjString*) fBGTrigClasses.At(i - fCollTrigClasses.GetEntries()))->String();
540   
541     gDirectory->mkdir(triggerClass);
542     gDirectory->cd(triggerClass);
543   
544     static_cast<AliTriggerAnalysis*> (fTriggerAnalysis.At(i))->SaveHistograms();
545     
546     gDirectory->cd("..");
547   }
548  
549   if (fBackgroundIdentification)
550   {
551     gDirectory->mkdir("background_identification");
552     gDirectory->cd("background_identification");
553       
554     fBackgroundIdentification->GetOutput()->Write();
555       
556     gDirectory->cd("..");
557   }
558   
559   if (folder)
560     gDirectory->cd("..");
561 }