ANALYSIS/AliBackgroundSelection.cxx

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