Added loop for extraction of clusters really attached to its track.
[u/mrichter/AliRoot.git] / ANALYSIS / AliBackgroundSelection.cxx
CommitLineData
b27a54cf 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
91bea6e7 23#include "AliBackgroundSelection.h"
24#include "TH2F.h"
25#include "TList.h"
91bea6e7 26#include "TString.h"
27#include "AliESDInputHandlerRP.h"
28#include "AliAnalysisManager.h"
29#include "TTree.h"
b27a54cf 30#include "AliMultiplicity.h"
91bea6e7 31#ifdef PASS1RECO
c93255fe 32#include "AliITSRecPoint.h"
91bea6e7 33#endif
b27a54cf 34
35
91bea6e7 36
37ClassImp(AliBackgroundSelection)
38
39AliBackgroundSelection::AliBackgroundSelection():
40 AliAnalysisCuts(), fOutputHist(0), fACut(0), fBCut(0), fDeltaPhiCut(10)
41{
b27a54cf 42 // ctor
91bea6e7 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
50AliBackgroundSelection::AliBackgroundSelection(const char* name, const char* title):
51 AliAnalysisCuts(name,title), fOutputHist(0), fACut(0), fBCut(0), fDeltaPhiCut(10)
52{
b27a54cf 53 // ctor
91bea6e7 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
62AliBackgroundSelection::AliBackgroundSelection(const AliBackgroundSelection& obj) : AliAnalysisCuts(obj),
63fOutputHist(0), fACut(0), fBCut(0), fDeltaPhiCut(0)
64{
b27a54cf 65 // copy ctor
91bea6e7 66 fOutputHist = obj.fOutputHist;
67 fACut = obj.fACut;
68 fBCut = obj.fBCut;
69 fDeltaPhiCut = obj.fDeltaPhiCut;
70}
71
72AliBackgroundSelection::~AliBackgroundSelection() {
b27a54cf 73 // dtor
91bea6e7 74 if(fOutputHist) {
75 delete fOutputHist;
76 fOutputHist = 0;
77 }
78
79}
80
b27a54cf 81Bool_t AliBackgroundSelection::IsSelected(TObject* const obj)
91bea6e7 82{
b27a54cf 83 // returns false if the event is identifiead as beam background,
84 // true otherwise.
91bea6e7 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
215void AliBackgroundSelection::Init(){
216
b27a54cf 217 // Set default cut values
91bea6e7 218 fACut = 65;
219 fBCut = 4;
220
221}
222
223
224void AliBackgroundSelection::BookClusterVsTrackletsHisto(const char * trigger_name){
225
b27a54cf 226 // Book control histogram for the cut on the correlation cluster vs tracklets
227
91bea6e7 228 Bool_t oldStatus = TH1::AddDirectoryStatus();
229 TH1::AddDirectory(kFALSE);
230
ce08cb1f 231 TH2F * h1 = new TH2F(GetClusterVsTrackletsHistoName(trigger_name),trigger_name, 300, -0.5, 2999.5, 1000, -0.5, 9999.5);
91bea6e7 232 h1->SetXTitle("Tracklets");
233 h1->SetYTitle("SPD Clusters");
adf3c99c 234 // AliInfo(Form("Creating histos: %s, all and accepted", GetClusterVsTrackletsHistoName(trigger_name)));
91bea6e7 235
236 TH2F * h2 = new TH2F(GetClusterVsTrackletsHistoNameAccepted(trigger_name),TString(trigger_name)+ "(accepted)",
ce08cb1f 237 300, -0.5, 2999.5, 1000, -0.5, 9999.5);
91bea6e7 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
248void AliBackgroundSelection::BookDeltaPhiHisto(const char * trigger_name){
249
b27a54cf 250 // Book control histogram for the cut on the DeltaPhi window used by vertexer Z
251
91bea6e7 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");
adf3c99c 257 // AliInfo(Form("Creating histos: %s, all and accepted", GetDeltaPhiHistoName(trigger_name)));
91bea6e7 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
270TH2F * AliBackgroundSelection::GetClusterVsTrackletsHisto(const char * trigger_name){
b27a54cf 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
91bea6e7 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}
285TH1F * AliBackgroundSelection::GetDeltaPhiHisto(const char * trigger_name){
b27a54cf 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
91bea6e7 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
301TH2F * AliBackgroundSelection::GetClusterVsTrackletsHistoAccepted(const char * trigger_name){
302
b27a54cf 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
91bea6e7 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
318TH1F * AliBackgroundSelection::GetDeltaPhiHistoAccepted(const char * trigger_name){
319
b27a54cf 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
91bea6e7 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
335const char * AliBackgroundSelection::GetClusterVsTrackletsHistoName(const char * trigger_name){
b27a54cf 336
337 // build up the name of the cluster vs tracklets histo using the trigger class
338
91bea6e7 339 static TString str;
340 str = ("hCvsT");
341 str = str+GetName()+"_"+trigger_name;
342 return str.Data();
343}
344
345const char * AliBackgroundSelection::GetClusterVsTrackletsHistoNameAccepted(const char * trigger_name){
b27a54cf 346
347 // build up the name of the cluster vs tracklets histo using the trigger class (accepted events)
91bea6e7 348 static TString str;
349 str = ("hCvsT");
350 str = str+GetName()+"_"+trigger_name + "_accepted";
351 return str.Data();
352}
353
354const char * AliBackgroundSelection::GetDeltaPhiHistoName(const char * trigger_name){
b27a54cf 355
356 // build up the name of the delta phi histo using the trigger class
357
358
91bea6e7 359 static TString str;
360 str = ("hDeltaPhi");
361 str = str+GetName()+"_"+trigger_name;
362 return str.Data();
363}
364
365const char * AliBackgroundSelection::GetDeltaPhiHistoNameAccepted(const char * trigger_name){
b27a54cf 366
367 // build up the name of the delta phi histo using the trigger class (accepted events)
368
91bea6e7 369 static TString str;
370 str = ("hDeltaPhi");
371 str = str+GetName()+"_"+trigger_name + "_accepted";
372 return str.Data();
373}
374
b27a54cf 375Long64_t AliBackgroundSelection::Merge(TCollection* const list)
91bea6e7 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
0c6c629b 385 //AliInfo("Merging");
91bea6e7 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;
b27a54cf 418 Int_t maxLoops = hlist->GetSize() + fOutputHist->GetSize(); // In the worst case all of the histos will be different...
91bea6e7 419 while(areListsDifferent) {
b27a54cf 420 if(iloop>maxLoops) AliFatal("Infinite Loop?");
91bea6e7 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())){
0c6c629b 445 //AliInfo(Form("Adding object %s",hist->GetName()));
91bea6e7 446 foundDiffinThisIterStep = kTRUE;
447 TH1 * hclone = (TH1*) hist->Clone();
448 hclone->Reset();
449 otherlist->Add(hclone);
450 }
451 }
452
8663e267 453 delete iterlist;
91bea6e7 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