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