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