Changing to centrality dependent corrections
[u/mrichter/AliRoot.git] / ANALYSIS / AliBackgroundSelection.cxx
1 // ----------------------------------------------------------------
2 // AliBackgroundSelection
3 //
4 // This class implements to cuts to reject background events from the
5 // samples to be used in the physics analysis:
6 // 1. A linear cut on the correlation cluster vs tracklets
7 // 2. A cut on the delta phi window used by the vertexer Z
8 // The parameters used in both cuts can be set
9 // 
10 // The class also produces control histograms for all and accepted
11 // events, for each trigger class present in the data independently.
12 // Histograms are booked on the fly in the UserExec, whenever a new
13 // trigger class is found.
14 //
15 // After the first implementation, it was realized that the deltaphi
16 // cut is more a quality selection cut than an event selection cut, so
17 // it is effectively disabled by default.
18 //
19 // Author: Michele Floris, CERN
20 // ----------------------------------------------------------------
21
22
23 #include "AliBackgroundSelection.h"
24 #include "TH2F.h"
25 #include "TList.h"
26 #include "TString.h"
27 #include "AliESDInputHandlerRP.h"
28 #include "AliAnalysisManager.h"
29 #include "TTree.h"
30 #include "AliMultiplicity.h"
31 #ifdef PASS1RECO
32 #include "AliITSRecPoint.h"
33 #endif
34
35
36
37 ClassImp(AliBackgroundSelection)
38
39 AliBackgroundSelection::AliBackgroundSelection():
40   AliAnalysisCuts(), fOutputHist(0), fACut(0), fBCut(0), fDeltaPhiCut(10)
41 {
42   // ctor
43   fOutputHist = new TList();
44   fOutputHist->SetOwner();
45   fACut = 65;
46   fBCut = 4;
47   fDeltaPhiCut = 10; // effectively disabling delta phi cut by default
48 }
49
50 AliBackgroundSelection::AliBackgroundSelection(const char* name, const char* title):
51   AliAnalysisCuts(name,title), fOutputHist(0), fACut(0), fBCut(0), fDeltaPhiCut(10)
52 {
53   // ctor
54   fOutputHist = new TList();
55   fOutputHist->SetOwner();
56   fACut = 65;
57   fBCut = 4;
58   fDeltaPhiCut = 10; //  effectively disabling delta phi cut by default
59
60 }
61
62 AliBackgroundSelection::AliBackgroundSelection(const AliBackgroundSelection& obj) : AliAnalysisCuts(obj),
63 fOutputHist(0), fACut(0), fBCut(0), fDeltaPhiCut(0)
64 {
65   // copy ctor
66   fOutputHist  = obj.fOutputHist;
67   fACut        = obj.fACut;
68   fBCut        = obj.fBCut;
69   fDeltaPhiCut = obj.fDeltaPhiCut;
70 }
71
72 AliBackgroundSelection::~AliBackgroundSelection() {
73   // dtor
74   if(fOutputHist) {
75     delete fOutputHist;
76     fOutputHist = 0;
77   }
78
79 }
80
81 Bool_t AliBackgroundSelection::IsSelected(TObject* const obj) 
82 {
83   // returns false if the event is identifiead as beam background,
84   // true otherwise.
85
86   // reset fSelected
87   SetSelected(kFALSE);
88 #ifdef PASS1RECO
89   // Get rec points
90   AliESDInputHandlerRP* handlerRP = dynamic_cast<AliESDInputHandlerRP*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
91   if (!handlerRP)
92     AliFatal("Cannot get the AliESDInputHandlerRP");
93
94   TTree* itsClusterTree = handlerRP->GetTreeR("ITS");
95   if (!itsClusterTree){
96     AliError("Cannot get the ITS Cluster tree");
97     return kFALSE;
98   }
99   //    AliFatal("Cannot get the ITS Cluster tree");
100
101   TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
102   TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
103
104   itsClusterBranch->SetAddress(&itsClusters);
105
106   Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
107 #endif
108
109   AliESDEvent * esdEv = (AliESDEvent*) obj;
110
111 #ifdef PASS1RECO
112   Float_t deltaPhi = 0.0; // deltaPhi is not available in pass1
113
114   // Get # spd clusters and of tracklets
115   Int_t spdClusters=0;
116
117
118   // loop over the its subdetectors
119   for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
120
121     if (!itsClusterTree->GetEvent(iIts))
122       continue;
123
124     Int_t nClusters = itsClusters->GetEntriesFast();
125
126     // loop over clusters
127     while (nClusters--) {
128       AliITSRecPoint* cluster = (AliITSRecPoint*) itsClusters->UncheckedAt(nClusters);
129
130       Int_t layer = cluster->GetLayer();
131
132       if (layer < 3) { // SPD
133         spdClusters++;
134       }
135     }
136   }
137 #endif
138
139   const AliMultiplicity* mult = esdEv->GetMultiplicity();
140   if (!mult){
141     AliFatal("No multiplicity object"); // TODO: Should this be fatal?
142   }
143   Int_t ntracklet = mult->GetNumberOfTracklets();
144
145 #ifndef PASS1RECO
146   // get deltaphi if vertexer z
147   Float_t deltaPhi = 0.0;
148   // Get Vertex
149   const AliESDVertex * vtxESD = esdEv->GetPrimaryVertexSPD();
150   if(vtxESD) {
151     if (vtxESD->IsFromVertexerZ()) deltaPhi = vtxESD->GetDispersion(); // dispersion contains deltaphi in case of vertexer Z
152   }
153   else {
154     AliWarning("No Vertex");
155   }
156
157   
158
159   // compute number of spd clusters
160   Float_t spdClusters = 0;
161   for(Int_t ilayer = 0; ilayer < 2; ilayer++){
162     spdClusters += mult->GetNumberOfITSClusters(ilayer);
163   }
164 #endif
165
166   // Check cuts
167   Bool_t isCvsTOk     = kFALSE;
168   Bool_t isDeltaPhiOk = kFALSE;
169
170   Float_t limit = fACut + ntracklet * fBCut;  
171   if (spdClusters > limit)        isCvsTOk = kFALSE;
172   else                            isCvsTOk = kTRUE ;
173
174   if(deltaPhi > fDeltaPhiCut)     isDeltaPhiOk = kFALSE;
175   else                            isDeltaPhiOk = kTRUE ;
176
177   if (!isCvsTOk || !isDeltaPhiOk) SetSelected(kFALSE);
178   else                            SetSelected(kTRUE );
179
180   // Fill control histos for all trigger classes
181   TString trgClasses = esdEv->GetFiredTriggerClasses();
182   TObjArray * tokens = trgClasses.Tokenize(" ");
183   TIter iter(tokens);
184   while(TObjString * tok = (TObjString*) iter.Next()){
185     // clean up trigger name
186     TString trg = tok->GetString();
187     trg.Strip(TString::kTrailing, ' ');
188     trg.Strip(TString::kLeading, ' ');
189     
190     // cluster vs tracklets
191     TH2F * hCvsT = GetClusterVsTrackletsHisto(trg.Data());
192     TH2F * hCvsTa = GetClusterVsTrackletsHistoAccepted(trg.Data());
193     hCvsT->Fill(ntracklet,spdClusters);
194     if(isCvsTOk) hCvsTa->Fill(ntracklet,spdClusters);
195
196     // Delta phi
197     TH1F * hDeltaPhi = GetDeltaPhiHisto(trg.Data());
198     TH1F * hDeltaPhia = GetDeltaPhiHistoAccepted(trg.Data());
199     hDeltaPhi->Fill(deltaPhi);
200     if(isDeltaPhiOk) hDeltaPhia->Fill(deltaPhi);
201   }
202   if(tokens) delete tokens;
203   // return decision
204
205 #ifdef PASS1RECO
206   if(itsClusters) {
207     itsClusters->Delete();
208     delete itsClusters;
209   }
210 #endif 
211   return Selected();
212 }
213
214
215 void   AliBackgroundSelection::Init(){
216
217   // Set default cut values
218   fACut = 65;
219   fBCut = 4;
220
221 }
222
223
224 void AliBackgroundSelection::BookClusterVsTrackletsHisto(const char * trigger_name){
225
226   // Book control histogram for the cut on the correlation cluster vs tracklets
227
228   Bool_t oldStatus = TH1::AddDirectoryStatus();
229   TH1::AddDirectory(kFALSE);
230
231   TH2F * h1 = new TH2F(GetClusterVsTrackletsHistoName(trigger_name),trigger_name, 300, -0.5, 2999.5, 1000, -0.5, 9999.5);
232   h1->SetXTitle("Tracklets");
233   h1->SetYTitle("SPD Clusters");
234   // AliInfo(Form("Creating histos: %s, all and accepted", GetClusterVsTrackletsHistoName(trigger_name)));
235
236   TH2F * h2 = new TH2F(GetClusterVsTrackletsHistoNameAccepted(trigger_name),TString(trigger_name)+ "(accepted)", 
237                        300, -0.5, 2999.5, 1000, -0.5, 9999.5);
238   h2->SetXTitle("Tracklets");
239   h2->SetYTitle("SPD Clusters");
240
241   fOutputHist->Add(h1);
242   fOutputHist->Add(h2);
243
244   TH1::AddDirectory(oldStatus);
245
246 }
247
248 void AliBackgroundSelection::BookDeltaPhiHisto(const char * trigger_name){
249
250   // Book control histogram for the cut on the DeltaPhi window used by vertexer Z
251
252   Bool_t oldStatus = TH1::AddDirectoryStatus();
253   TH1::AddDirectory(kFALSE);
254
255   TH1F * h1 = new TH1F(GetDeltaPhiHistoName(trigger_name),trigger_name, 100,0,0.5);
256   h1->SetXTitle("#Delta #phi");
257   // AliInfo(Form("Creating histos: %s, all and accepted", GetDeltaPhiHistoName(trigger_name)));
258
259   TH1F * h2 = new TH1F(GetDeltaPhiHistoNameAccepted(trigger_name),TString(trigger_name)+ "(accepted)", 100,0,0.5);
260   h2->SetXTitle("#Delta #phi");
261
262
263   fOutputHist->Add(h1);
264   fOutputHist->Add(h2);
265
266   TH1::AddDirectory(oldStatus);
267
268 }
269
270 TH2F * AliBackgroundSelection::GetClusterVsTrackletsHisto(const char * trigger_name){
271
272   // Returns the control histogram corresponding to a given trigger
273   // class. If it does not exist, it creates it and adds it to the
274   // output list
275   // All Events
276
277   if(!fOutputHist) {AliError("List of histos not initialized");return 0;}
278   TH2F * h = (TH2F*) fOutputHist->FindObject(GetClusterVsTrackletsHistoName(trigger_name));  
279   if(!h) {
280     BookClusterVsTrackletsHisto(trigger_name);
281     h = (TH2F*) fOutputHist->FindObject(GetClusterVsTrackletsHistoName(trigger_name));  
282   }
283   return h;
284 }
285 TH1F * AliBackgroundSelection::GetDeltaPhiHisto(const char * trigger_name){
286
287   // Returns the control histogram corresponding to a given trigger
288   // class. If it does not exist, it creates it and adds it to the
289   // output list
290   // All Events
291
292   if(!fOutputHist) {AliError("List of histos not initialized");return 0;}
293   TH1F * h = (TH1F*) fOutputHist->FindObject(GetDeltaPhiHistoName(trigger_name));  
294   if(!h) {
295     BookDeltaPhiHisto(trigger_name);
296     h  = (TH1F*) fOutputHist->FindObject(GetDeltaPhiHistoName(trigger_name));  
297   }
298   return h;
299 }
300
301 TH2F * AliBackgroundSelection::GetClusterVsTrackletsHistoAccepted(const char * trigger_name){
302
303   // Returns the control histogram corresponding to a given trigger
304   // class. If it does not exist, it creates it and adds it to the
305   // output list
306   // Events passing the cut only
307
308   if(!fOutputHist) {AliError("List of histos not initialized");return 0;}
309   TH2F * h = (TH2F*) fOutputHist->FindObject(GetClusterVsTrackletsHistoNameAccepted(trigger_name));
310   if(!h) {
311     BookClusterVsTrackletsHisto(trigger_name);
312     h = (TH2F*) fOutputHist->FindObject(GetClusterVsTrackletsHistoNameAccepted(trigger_name));  
313   }
314   return h;
315   
316 }
317
318 TH1F * AliBackgroundSelection::GetDeltaPhiHistoAccepted(const char * trigger_name){
319
320   // Returns the control histogram corresponding to a given trigger
321   // class. If it does not exist, it creates it and adds it to the
322   // output list
323   // Events passing the cut only
324
325   if(!fOutputHist) {AliError("List of histos not initialized");return 0;}
326   TH1F * h = (TH1F*) fOutputHist->FindObject(GetDeltaPhiHistoNameAccepted(trigger_name));  
327   if(!h) {
328     BookDeltaPhiHisto(trigger_name);
329     h  = (TH1F*) fOutputHist->FindObject(GetDeltaPhiHistoNameAccepted(trigger_name));  
330   }
331   return h;
332   
333 }
334
335 const char * AliBackgroundSelection::GetClusterVsTrackletsHistoName(const char * trigger_name){
336
337   // build up the name of the cluster vs tracklets histo using the trigger class
338
339     static TString str;
340     str = ("hCvsT");
341     str = str+GetName()+"_"+trigger_name;
342     return str.Data();
343 }
344
345 const char * AliBackgroundSelection::GetClusterVsTrackletsHistoNameAccepted(const char * trigger_name){
346
347   // build up the name of the cluster vs tracklets histo using the trigger class (accepted events)
348     static TString str;
349     str = ("hCvsT");
350     str = str+GetName()+"_"+trigger_name + "_accepted";
351     return str.Data();
352 }
353
354 const char * AliBackgroundSelection::GetDeltaPhiHistoName(const char * trigger_name){
355
356   // build up the name of the delta phi histo using the trigger class
357
358
359     static TString str;
360     str = ("hDeltaPhi");
361     str = str+GetName()+"_"+trigger_name;
362     return str.Data();
363 }
364
365 const char * AliBackgroundSelection::GetDeltaPhiHistoNameAccepted(const char * trigger_name){
366
367   // build up the name of the delta phi histo using the trigger class (accepted events)
368
369     static TString str;
370     str = ("hDeltaPhi");
371     str = str+GetName()+"_"+trigger_name + "_accepted";
372     return str.Data();
373 }
374
375 Long64_t AliBackgroundSelection::Merge(TCollection* const list)
376 {
377   // Merge a list of AliBackgroundSelection objects with this (needed for
378   // PROOF).
379   // Returns the number of merged objects (including this).
380
381   // We have to make sure that all the list contain the same histos in
382   // the same order. We thus also have to sort the list (sorting is
383   // done by name in TList).
384
385   //AliInfo("Merging");
386
387   if (!list)
388     return 0;
389
390   if (list->IsEmpty())
391     return 1;
392
393   TIterator* iter = list->MakeIterator();
394   TObject* obj;
395
396   // collections of all histograms
397   const Int_t nHists = 1;
398   TList collections[nHists];
399
400   Int_t count = 0;
401   // 1. Sort this list
402   fOutputHist->Sort();
403   
404   while ((obj = iter->Next())) {
405     Bool_t foundDiffinThisIterStep = kFALSE;
406     //    Printf("%d - %s",count, obj->GetName());
407     AliBackgroundSelection* entry = dynamic_cast<AliBackgroundSelection*> (obj);
408     if (entry == 0) 
409       continue;
410
411     TList * hlist = entry->fOutputHist;
412
413     // Check if all histos in this fOutputHist are also in the one from entry and viceversa
414     // Use getters to automatically book non defined histos    
415
416     Bool_t areListsDifferent=kTRUE;
417     Int_t iloop = 0;
418     Int_t maxLoops = hlist->GetSize() + fOutputHist->GetSize(); // In the worst case all of the histos will be different...    
419     while(areListsDifferent) {
420       if(iloop>maxLoops) AliFatal("Infinite Loop?");
421       iloop++;
422       // sort
423       hlist->Sort();
424       fOutputHist->Sort();
425       // loop over the largest 
426
427       // loop over the largest 
428       TObject * hist =0;
429       TIterator * iterlist = 0;
430       TList * thislist  = 0; // the list over which I'm iterating (i.e. the largest)
431       TList * otherlist = 0; // the other list
432
433       if (hlist->GetSize() >= fOutputHist->GetSize()) { 
434         thislist  = hlist;
435         otherlist = fOutputHist;
436       }
437       else{
438         thislist  = fOutputHist;
439         otherlist = hlist;      
440       }
441       iterlist = thislist->MakeIterator();
442
443       while ((hist= iterlist->Next())){ 
444         if(!otherlist->FindObject(hist->GetName())){
445           //AliInfo(Form("Adding object %s",hist->GetName()));
446           foundDiffinThisIterStep = kTRUE;
447           TH1 * hclone =  (TH1*) hist->Clone();
448           hclone->Reset();
449           otherlist->Add(hclone);
450         }
451       }
452
453       delete iterlist;
454       // re-sort before checking
455       hlist->Sort();
456       fOutputHist->Sort();
457
458       // check if everything is fine    
459       areListsDifferent=kFALSE;
460       if (hlist->GetSize() == fOutputHist->GetSize()) { 
461         Int_t nhist =  fOutputHist->GetSize();
462         for(Int_t ihist = 0; ihist < nhist; ihist++){
463           if(strcmp(fOutputHist->At(ihist)->GetName(),hlist->At(ihist)->GetName())) areListsDifferent = kTRUE;
464         }
465       } else {
466         areListsDifferent=kTRUE;
467       }
468     }
469
470     // last check: if something is not ok die loudly 
471     if (hlist->GetSize() != fOutputHist->GetSize()) {
472       AliFatal("Mismatching size!");
473     }
474     Int_t nhist =  fOutputHist->GetSize();
475     for(Int_t ihist = 0; ihist < nhist; ihist++){
476       if(strcmp(fOutputHist->At(ihist)->GetName(),hlist->At(ihist)->GetName())){
477         AliFatal(Form("Mismatching histos: %s -> %s", fOutputHist->At(ihist)->GetName(),hlist->At(ihist)->GetName()));
478       }
479     }
480
481     if (foundDiffinThisIterStep){
482       iter->Reset(); // We found a difference: previous lists could
483                      // also be affected... We start from scratch
484       Int_t n = 0;
485       collections[n++].Clear();
486       count = 0;
487     }
488     else {
489 //       AliInfo("hlist");
490 //       hlist->Print();
491 //       AliInfo("fOutputHist");
492 //       fOutputHist->Print();
493       
494       Int_t n = 0;
495       collections[n++].Add(hlist);
496       
497       count++;
498     }
499   }
500
501   Int_t n = 0;
502   fOutputHist->Merge(&collections[n++]);
503   
504   delete iter;
505
506   return count+1;
507 }
508
509