PWGHF/vertexingHF/AliNormalizationCounter.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //*************************************************************************
  19 // Class AliNormalizationCounter
  20 // Class to store the informations relevant for the normalization in the
  21 // barrel for each run
  22 // Authors: G. Ortona, ortona@to.infn.it
  23 // D. Caffarri, davide.caffarri@pd.to.infn.it
  24 // with many thanks to P. Pillot
  25 /////////////////////////////////////////////////////////////
  26
  27 #include "AliLog.h"
  28 #include "AliNormalizationCounter.h"
  29 #include <AliESDEvent.h>
  30 #include <AliESDtrack.h>
  31 #include <AliAODEvent.h>
  32 #include <AliVParticle.h>
  33 #include <AliTriggerAnalysis.h>
  34 #include <TH1F.h>
  35 #include <TH2F.h>
  36 #include <TList.h>
  37 #include <TString.h>
  38 #include <TCanvas.h>
  39 #include <AliPhysicsSelection.h>
  40 #include <AliMultiplicity.h>
  41
  42 ClassImp(AliNormalizationCounter)
  43
  44 //____________________________________________
  45 AliNormalizationCounter::AliNormalizationCounter():
  46 TNamed(),
  47 fCounters(),
  48 fESD(kFALSE),
  49 fMultiplicity(kFALSE),
  50 fMultiplicityEtaRange(1.0),
  51 fHistTrackFilterEvMult(0),
  52 fHistTrackAnaEvMult(0),
  53 fHistTrackFilterSpdMult(0),
  54 fHistTrackAnaSpdMult(0)
  55 {
  56   // empty constructor
  57 }
  58
  59 //__________________________________________________
  60 AliNormalizationCounter::AliNormalizationCounter(const char *name):
  61 TNamed(name,name),
  62 fCounters(name),
  63 fESD(kFALSE),
  64 fMultiplicity(kFALSE),
  65 fMultiplicityEtaRange(1.0),
  66 fHistTrackFilterEvMult(0),
  67 fHistTrackAnaEvMult(0),
  68 fHistTrackFilterSpdMult(0),
  69 fHistTrackAnaSpdMult(0)
  70 {
  71   ;
  72 }
  73
  74 //______________________________________________
  75 AliNormalizationCounter::~AliNormalizationCounter()
  76 {
  77   //destructor
  78   if(fHistTrackFilterEvMult){
  79     delete fHistTrackFilterEvMult;
  80     fHistTrackFilterEvMult =0;
  81   }
  82   if(fHistTrackAnaEvMult){
  83     delete fHistTrackAnaEvMult;
  84     fHistTrackAnaEvMult=0;
  85   }
  86   if(fHistTrackFilterSpdMult){
  87     delete fHistTrackFilterSpdMult;
  88     fHistTrackFilterSpdMult=0;
  89   }
  90   if(fHistTrackAnaSpdMult){
  91     delete fHistTrackAnaSpdMult;
  92     fHistTrackAnaSpdMult=0;
  93   }
  94 }
  95
  96 //______________________________________________
  97 void AliNormalizationCounter::Init()
  98 {
  99   //variables initialization
 100   fCounters.AddRubric("Event","triggered/V0AND/PileUp/PbPbC0SMH-B-NOPF-ALLNOTRD/Candles0.3/PrimaryV/countForNorm/noPrimaryV/zvtxGT10/!V0A&Candle03/!V0A&PrimaryV/Candid(Filter)/Candid(Analysis)/NCandid(Filter)/NCandid(Analysis)");
 101   if(fMultiplicity)  fCounters.AddRubric("Multiplicity", 5000);
 102   fCounters.AddRubric("Run", 1000000);
 103   fCounters.Init();
 104   fHistTrackFilterEvMult=new TH2F("FiltCandidvsTracksinEv","FiltCandidvsTracksinEv",10000,-0.5,9999.5,200,-0.5,199.5);
 105   fHistTrackFilterEvMult->GetYaxis()->SetTitle("NCandidates");
 106   fHistTrackFilterEvMult->GetXaxis()->SetTitle("NTracksinEvent");
 107   fHistTrackAnaEvMult=new TH2F("AnaCandidvsTracksinEv","AnaCandidvsTracksinEv",10000,-0.5,9999.5,100,-0.5,99.5);
 108   fHistTrackAnaEvMult->GetYaxis()->SetTitle("NCandidates");
 109   fHistTrackAnaEvMult->GetXaxis()->SetTitle("NTracksinEvent");
 110   fHistTrackFilterSpdMult=new TH2F("FilterCandidvsSpdMult","FilterCandidvsSpdMult",5000,-0.5,4999.5,200,-0.5,199.5);
 111   fHistTrackFilterSpdMult->GetYaxis()->SetTitle("NCandidates");
 112   fHistTrackFilterSpdMult->GetXaxis()->SetTitle("NSPDTracklets");
 113   fHistTrackAnaSpdMult=new TH2F("AnaCandidvsSpdMult","AnaCandidvsSpdMult",5000,-0.5,4999.5,100,-0.5,99.5);
 114   fHistTrackAnaSpdMult->GetYaxis()->SetTitle("NCandidates");
 115   fHistTrackAnaSpdMult->GetXaxis()->SetTitle("NSPDTracklets");
 116 }
 117
 118 //______________________________________________
 119 Long64_t AliNormalizationCounter::Merge(TCollection* list){
 120   if (!list) return 0;
 121   if (list->IsEmpty()) return 0;//(Long64_t)fCounters.Merge(list);
 122
 123   TIter next(list);
 124   const TObject* obj = 0x0;
 125   while ((obj = next())) {
 126
 127     // check that "obj" is an object of the class AliNormalizationCounter
 128     const AliNormalizationCounter* counter = dynamic_cast<const AliNormalizationCounter*>(obj);
 129     if (!counter) {
 130       AliError(Form("object named %s is not AliNormalizationCounter! Skipping it.", counter->GetName()));
 131       continue;
 132     }
 133
 134     Add(counter);
 135
 136   }
 137
 138   return (Long64_t)1;//(Long64_t)fCounters->GetEntries();
 139 }
 140 //_______________________________________
 141 void AliNormalizationCounter::Add(const AliNormalizationCounter *norm){
 142   fCounters.Add(&(norm->fCounters));
 143   fHistTrackFilterEvMult->Add(norm->fHistTrackFilterEvMult);
 144   fHistTrackAnaEvMult->Add(norm->fHistTrackAnaEvMult);
 145   fHistTrackFilterSpdMult->Add(norm->fHistTrackFilterSpdMult);
 146   fHistTrackAnaSpdMult->Add(norm->fHistTrackAnaSpdMult);
 147 }
 148 //_______________________________________
 149 /*
 150 Stores the variables used for normalization as function of run number
 151 returns kTRUE if the event is to be counted for normalization
 152 (pass event selection cuts OR has no primary vertex)
 153  */
 154 void AliNormalizationCounter::StoreEvent(AliVEvent *event,AliRDHFCuts *rdCut,Bool_t mc, Int_t multiplicity){
 155   //
 156
 157   Bool_t isEventSelected = rdCut->IsEventSelected(event);
 158
 159   // events not passing physics selection. do nothing
 160   if(rdCut->IsEventRejectedDuePhysicsSelection()) return;
 161
 162   Bool_t v0A=kFALSE;
 163   Bool_t v0B=kFALSE;
 164   Bool_t flag03=kFALSE;
 165   Bool_t flagPV=kFALSE;
 166
 167   //Run Number
 168   Int_t runNumber = event->GetRunNumber();
 169
 170   // Evaluate the multiplicity
 171   if(multiplicity==-9999) Multiplicity(event);
 172
 173   //Find CINT1B
 174   AliESDEvent *eventESD = (AliESDEvent*)event;
 175   if(!eventESD){AliError("ESD event not available");return;}
 176   if(mc&&event->GetEventType() != 0)return;
 177   //event must be either physics or MC
 178   if(!(event->GetEventType() == 7||event->GetEventType() == 0))return;
 179
 180   if(fMultiplicity)
 181     fCounters.Count(Form("Event:triggered/Run:%d/Multiplicity:%d",runNumber,multiplicity));
 182   else
 183     fCounters.Count(Form("Event:triggered/Run:%d",runNumber));
 184
 185   //Find V0AND
 186   AliTriggerAnalysis trAn; /// Trigger Analysis
 187   v0B = trAn.IsOfflineTriggerFired(eventESD , AliTriggerAnalysis::kV0C);
 188   v0A = trAn.IsOfflineTriggerFired(eventESD , AliTriggerAnalysis::kV0A);
 189   if(v0A&&v0B){
 190     if(fMultiplicity)
 191       fCounters.Count(Form("Event:V0AND/Run:%d/Multiplicity:%d",runNumber,multiplicity));
 192     else
 193       fCounters.Count(Form("Event:V0AND/Run:%d",runNumber));
 194   }
 195
 196   //FindPrimary vertex
 197   // AliVVertex *vtrc =  (AliVVertex*)event->GetPrimaryVertex();
 198   // if(vtrc && vtrc->GetNContributors()>0){
 199   //   fCounters.Count(Form("Event:PrimaryV/Run:%d",runNumber));
 200   //   flagPV=kTRUE;
 201   // }
 202
 203   //trigger
 204   AliAODEvent *eventAOD = (AliAODEvent*)event;
 205   TString trigclass=eventAOD->GetFiredTriggerClasses();
 206   if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD")||trigclass.Contains("C0SMH-B-NOPF-ALL")){
 207     if(fMultiplicity)
 208       fCounters.Count(Form("Event:PbPbC0SMH-B-NOPF-ALLNOTRD/Run:%d/Multiplicity:%d",runNumber,multiplicity));
 209     else
 210       fCounters.Count(Form("Event:PbPbC0SMH-B-NOPF-ALLNOTRD/Run:%d",runNumber));
 211   }
 212
 213   //FindPrimary vertex
 214   if(isEventSelected){
 215     if(fMultiplicity)
 216       fCounters.Count(Form("Event:PrimaryV/Run:%d/Multiplicity:%d",runNumber,multiplicity));
 217     else
 218       fCounters.Count(Form("Event:PrimaryV/Run:%d",runNumber));
 219     flagPV=kTRUE;
 220   }else{
 221     if(rdCut->GetWhyRejection()==0){
 222       if(fMultiplicity)
 223         fCounters.Count(Form("Event:noPrimaryV/Run:%d/Multiplicity:%d",runNumber,multiplicity));
 224       else
 225         fCounters.Count(Form("Event:noPrimaryV/Run:%d",runNumber));
 226     }
 227     //find good vtx outside range
 228     if(rdCut->GetWhyRejection()==6){
 229       if(fMultiplicity) {
 230         fCounters.Count(Form("Event:zvtxGT10/Run:%d/Multiplicity:%d",runNumber,multiplicity));
 231         fCounters.Count(Form("Event:PrimaryV/Run:%d/Multiplicity:%d",runNumber,multiplicity));
 232       } else {
 233         fCounters.Count(Form("Event:zvtxGT10/Run:%d",runNumber));
 234         fCounters.Count(Form("Event:PrimaryV/Run:%d",runNumber));
 235       }
 236       flagPV=kTRUE;
 237     }
 238     if(rdCut->GetWhyRejection()==1){
 239       if(fMultiplicity)
 240         fCounters.Count(Form("Event:PileUp/Run:%d/Multiplicity:%d",runNumber,multiplicity));
 241       else
 242         fCounters.Count(Form("Event:PileUp/Run:%d",runNumber));
 243     }
 244   }
 245   //to be counted for normalization
 246   if(rdCut->CountEventForNormalization()){
 247     if(fMultiplicity)
 248       fCounters.Count(Form("Event:countForNorm/Run:%d/Multiplicity:%d",runNumber,multiplicity));
 249     else
 250       fCounters.Count(Form("Event:countForNorm/Run:%d",runNumber));
 251   }
 252
 253
 254   //Find Candle
 255   Int_t trkEntries = (Int_t)event->GetNumberOfTracks();
 256   for(Int_t i=0;i<trkEntries&&!flag03;i++){
 257     AliAODTrack *track=(AliAODTrack*)event->GetTrack(i);
 258     if((track->Pt()>0.3)&&(!flag03)){
 259       if(fMultiplicity)
 260         fCounters.Count(Form("Event:Candles0.3/Run:%d/Multiplicity:%d",runNumber,multiplicity));
 261       else
 262         fCounters.Count(Form("Event:Candles0.3/Run:%d",runNumber));
 263       flag03=kTRUE;
 264       break;
 265     }
 266   }
 267
 268   if(!(v0A&&v0B)&&(flag03)){
 269     if(fMultiplicity)
 270       fCounters.Count(Form("Event:!V0A&Candle03/Run:%d/Multiplicity:%d",runNumber,multiplicity));
 271     else
 272       fCounters.Count(Form("Event:!V0A&Candle03/Run:%d",runNumber));
 273   }
 274   if(!(v0A&&v0B)&&flagPV){
 275     if(fMultiplicity)
 276       fCounters.Count(Form("Event:!V0A&PrimaryV/Run:%d/Multiplicity:%d",runNumber,multiplicity));
 277     else
 278       fCounters.Count(Form("Event:!V0A&PrimaryV/Run:%d",runNumber));
 279   }
 280
 281   return;
 282 }
 283 //_____________________________________________________________________
 284 void AliNormalizationCounter::StoreCandidates(AliVEvent *event,Int_t nCand,Bool_t flagFilter){
 285
 286   Int_t ntracks=event->GetNumberOfTracks();
 287   if(flagFilter)fHistTrackFilterEvMult->Fill(ntracks,nCand);
 288   else fHistTrackAnaEvMult->Fill(ntracks,nCand);
 289   Int_t nSPD=0;
 290   if(fESD){
 291     AliESDEvent *ESDevent=(AliESDEvent*)event;
 292     const AliMultiplicity *alimult = ESDevent->GetMultiplicity();
 293     nSPD = alimult->GetNumberOfTracklets();
 294
 295   }else{
 296     AliAODEvent *aodEvent =(AliAODEvent*)event;
 297     AliAODTracklets *trklets=aodEvent->GetTracklets();
 298     nSPD = trklets->GetNumberOfTracklets();
 299   }
 300   if(flagFilter)fHistTrackFilterSpdMult->Fill(nSPD,nCand);
 301   else fHistTrackAnaSpdMult->Fill(nSPD,nCand);
 302
 303   Int_t runNumber = event->GetRunNumber();
 304   Int_t multiplicity = Multiplicity(event);
 305   if(nCand==0)return;
 306   if(flagFilter){
 307     if(fMultiplicity)
 308       fCounters.Count(Form("Event:Candid(Filter)/Run:%d/Multiplicity:%d",runNumber,multiplicity));
 309     else
 310       fCounters.Count(Form("Event:Candid(Filter)/Run:%d",runNumber));
 311     for(Int_t i=0;i<nCand;i++){
 312       if(fMultiplicity)
 313         fCounters.Count(Form("Event:NCandid(Filter)/Run:%d/Multiplicity:%d",runNumber,multiplicity));
 314       else
 315         fCounters.Count(Form("Event:NCandid(Filter)/Run:%d",runNumber));
 316     }
 317   }else{
 318     if(fMultiplicity)
 319       fCounters.Count(Form("Event:Candid(Analysis)/Run:%d/Multiplicity:%d",runNumber,multiplicity));
 320     else
 321       fCounters.Count(Form("Event:Candid(Analysis)/Run:%d",runNumber));
 322     for(Int_t i=0;i<nCand;i++){
 323       if(fMultiplicity)
 324         fCounters.Count(Form("Event:NCandid(Analysis)/Run:%d/Multiplicity:%d",runNumber,multiplicity));
 325       else
 326         fCounters.Count(Form("Event:NCandid(Analysis)/Run:%d",runNumber));
 327     }
 328   }
 329   return;
 330 }
 331 //_______________________________________________________________________
 332 TH1D* AliNormalizationCounter::DrawAgainstRuns(TString candle,Bool_t drawHist){
 333   //
 334   fCounters.SortRubric("Run");
 335   TString selection;
 336   selection.Form("event:%s",candle.Data());
 337   TH1D* histoneD = fCounters.Get("run",selection.Data());
 338
 339   histoneD->Sumw2();
 340   if(drawHist)histoneD->DrawClone();
 341   return histoneD;
 342 }
 343 //___________________________________________________________________________
 344 TH1D* AliNormalizationCounter::DrawRatio(TString candle1,TString candle2){
 345   //
 346   fCounters.SortRubric("Run");
 347   TString name;
 348
 349   name.Form("%s/%s",candle1.Data(),candle2.Data());
 350   TH1D* num=DrawAgainstRuns(candle1.Data(),kFALSE);
 351   TH1D* den=DrawAgainstRuns(candle2.Data(),kFALSE);
 352
 353   den->SetTitle(candle2.Data());
 354   den->SetName(candle2.Data());
 355   num->Divide(num,den,1,1,"B");
 356   num->SetTitle(name.Data());
 357   num->SetName(name.Data());
 358   num->DrawClone();
 359   return num;
 360 }
 361 //___________________________________________________________________________
 362 void AliNormalizationCounter::PrintRubrics(){
 363   fCounters.PrintKeyWords();
 364 }
 365 //___________________________________________________________________________
 366 Double_t AliNormalizationCounter::GetSum(TString candle){
 367   TString selection="event:";
 368   selection.Append(candle);
 369   return fCounters.GetSum(selection.Data());
 370 }
 371 //___________________________________________________________________________
 372 TH2F* AliNormalizationCounter::GetHist(Bool_t filtercuts,Bool_t spdtracklets,Bool_t drawHist){
 373   if(filtercuts){
 374     if(spdtracklets){
 375       if(drawHist)fHistTrackFilterSpdMult->DrawCopy("LEGO2Z 0");
 376       return fHistTrackFilterSpdMult;
 377     }else{
 378       if(drawHist)fHistTrackFilterEvMult->DrawCopy("LEGO2Z 0");
 379       return fHistTrackFilterEvMult;
 380     }
 381   }else{
 382     if(spdtracklets){
 383       if(drawHist)fHistTrackAnaSpdMult->DrawCopy("LEGO2Z 0");
 384       return fHistTrackAnaSpdMult;
 385     }else{
 386       if(drawHist)fHistTrackAnaEvMult->DrawCopy("LEGO2Z 0");
 387       return fHistTrackAnaEvMult;
 388     }
 389   }
 390 }
 391 //___________________________________________________________________________
 392 Double_t AliNormalizationCounter::GetNEventsForNorm(){
 393   Double_t noVtxzGT10=GetSum("noPrimaryV")*GetSum("zvtxGT10")/GetSum("PrimaryV");
 394   return GetSum("countForNorm")-noVtxzGT10;
 395 }
 396 //___________________________________________________________________________
 397 Double_t AliNormalizationCounter::GetNEventsForNorm(Int_t runnumber){
 398   TString listofruns = fCounters.GetKeyWords("RUN");
 399   if(!listofruns.Contains(Form("%d",runnumber))){
 400     printf("WARNING: %d is not a valid run number\n",runnumber);
 401     fCounters.Print("Run","",kTRUE);
 402     return 0.;
 403   }
 404   TString suffix;suffix.Form("/RUN:%d",runnumber);
 405   TString zvtx;zvtx.Form("zvtxGT10%s",suffix.Data());
 406   TString noPV;noPV.Form("noPrimaryV%s",suffix.Data());
 407   TString pV;pV.Form("PrimaryV%s",suffix.Data());
 408   TString tbc;tbc.Form("countForNorm%s",suffix.Data());
 409   Double_t noVtxzGT10=GetSum(noPV.Data())*GetSum(zvtx.Data())/GetSum(pV.Data());
 410   return GetSum(tbc.Data())-noVtxzGT10;
 411 }
 412
 413 //___________________________________________________________________________
 414 Double_t AliNormalizationCounter::GetNEventsForNorm(Int_t minmultiplicity, Int_t maxmultiplicity){
 415
 416   if(!fMultiplicity) {
 417     AliInfo("Sorry, you didn't activate the multiplicity in the counter!");
 418     return 0.;
 419   }
 420
 421   TString listofruns = fCounters.GetKeyWords("Multiplicity");
 422
 423   Int_t nmultbins = maxmultiplicity - minmultiplicity;
 424   Double_t sumnoPV=0., sumZvtx=0., sumPv=0., sumEvtNorm=0.;
 425   for (Int_t ibin=0; ibin<=nmultbins; ibin++) {
 426     //    cout << " Looking at bin "<< ibin+minmultiplicity<<endl;
 427     if(!listofruns.Contains(Form("%d",ibin+minmultiplicity))){
 428       AliInfo(Form("WARNING: %d is not a valid multiplicity number. \n",ibin+minmultiplicity));
 429       continue;
 430     }
 431     TString suffix;suffix.Form("/Multiplicity:%d",ibin+minmultiplicity);
 432     TString zvtx;zvtx.Form("zvtxGT10%s",suffix.Data());
 433     TString noPV;noPV.Form("noPrimaryV%s",suffix.Data());
 434     TString pV;pV.Form("PrimaryV%s",suffix.Data());
 435     TString tbc;tbc.Form("countForNorm%s",suffix.Data());
 436     sumnoPV += GetSum(noPV.Data());
 437     sumZvtx += GetSum(zvtx.Data());
 438     sumPv += GetSum(pV.Data());
 439     sumEvtNorm += GetSum(tbc.Data());
 440   }
 441   Double_t noVtxzGT10 = sumPv>0. ? sumnoPV * sumZvtx / sumPv : 0.;
 442   return sumEvtNorm - noVtxzGT10;
 443 }
 444
 445 //___________________________________________________________________________
 446 TH1D* AliNormalizationCounter::DrawNEventsForNorm(Bool_t drawRatio){
 447   //usare algebra histos
 448   fCounters.SortRubric("Run");
 449   TString selection;
 450
 451   selection.Form("event:noPrimaryV");
 452   TH1D* hnoPrimV = fCounters.Get("run",selection.Data());
 453   hnoPrimV->Sumw2();
 454
 455   selection.Form("event:zvtxGT10");
 456   TH1D*  hzvtx= fCounters.Get("run",selection.Data());
 457   hzvtx->Sumw2();
 458
 459   selection.Form("event:PrimaryV");
 460   TH1D* hPrimV = fCounters.Get("run",selection.Data());
 461   hPrimV->Sumw2();
 462
 463   hzvtx->Multiply(hnoPrimV);
 464   hzvtx->Divide(hPrimV);
 465
 466   selection.Form("event:countForNorm");
 467   TH1D* hCountForNorm = fCounters.Get("run",selection.Data());
 468   hCountForNorm->Sumw2();
 469
 470   hCountForNorm->Add(hzvtx,-1.);
 471
 472   if(drawRatio){
 473     selection.Form("event:triggered");
 474     TH1D* htriggered = fCounters.Get("run",selection.Data());
 475     htriggered->Sumw2();
 476     hCountForNorm->Divide(htriggered);
 477   }
 478
 479   hCountForNorm->DrawClone();
 480   return hCountForNorm;
 481 }
 482
 483 //___________________________________________________________________________
 484 Int_t AliNormalizationCounter::Multiplicity(AliVEvent* event){
 485
 486   Int_t multiplicity = 0;
 487   AliAODEvent *eventAOD = (AliAODEvent*)event;
 488   AliAODTracklets * aodTracklets = (AliAODTracklets*)eventAOD->GetTracklets();
 489   Int_t ntracklets = (Int_t)aodTracklets->GetNumberOfTracklets();
 490   for(Int_t i=0;i<ntracklets; i++){
 491     Double_t theta = aodTracklets->GetTheta(i);
 492     Double_t eta = -TMath::Log( TMath::Tan(theta/2.) ); // check the formula
 493     if(TMath::Abs(eta)<fMultiplicityEtaRange){ // set the proper cut on eta
 494       multiplicity++;
 495     }
 496   }
 497
 498   return multiplicity;
 499 }