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