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