PWG/CaloTrackCorrBase/AliAnalysisTaskCounter.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 //_________________________________________________________________________
  17 // Count events with different selections
  18 //
  19 // It produces a histogram with the number of events with 9 bins:
  20 // 0: all events (that passed the physics selection if it was on)
  21 // 1: same but cross check that event pointer did exist
  22 // 2: passes vertex cut
  23 // 3: passes track number cut, tracks for eta < 0.8
  24 // 4: 3 && 2
  25 // 5: pass VAND
  26 // 6: 5 && 2
  27 // 7: 5 && 3
  28 // 8: 5 && 3 && 2
  29 // 9: not pileup from SPD
  30 // 10: Good vertex
  31 // 11: 10 && 5
  32 // 12: 10 && 3
  33 // 13: 10 && 2
  34 // 14: 10 && 2 && 3 && 5
  35 // 15: 10 && 9
  36 // 16: 9  && 5
  37 //
  38 // Author: Gustavo Conesa Balbastre (LPSC)
  39 //
  40 //_________________________________________________________________________
  41
  42 #include "TH2F.h"
  43 #include "AliAODHeader.h"
  44 #include "AliTriggerAnalysis.h"
  45 #include "AliESDEvent.h"
  46 #include "AliAODEvent.h"
  47 #include "AliESDtrackCuts.h"
  48 #include "AliAnalysisManager.h"
  49 #include "AliInputEventHandler.h"
  50
  51 #include "AliAnalysisTaskCounter.h"
  52 ClassImp(AliAnalysisTaskCounter)
  53
  54 //______________________________________________________________
  55 AliAnalysisTaskCounter::AliAnalysisTaskCounter(const char *name)
  56 : AliAnalysisTaskSE(name),
  57   fAcceptFastCluster(kTRUE),
  58   fZVertexCut(10.),
  59   fTrackMultEtaCut(0.8),
  60   fCaloFilterPatch(kFALSE),
  61   fOutputContainer(0x0),
  62   fESDtrackCuts(AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()),
  63   fTriggerAnalysis (new AliTriggerAnalysis),
  64   fhNEvents(0),
  65   fhXVertex(0),    fhYVertex(0),    fhZVertex(0),
  66   fhXGoodVertex(0),fhYGoodVertex(0),fhZGoodVertex(0),
  67   fhCentrality(0), fhEventPlaneAngle(0)
  68 {
  69   //ctor
  70   DefineOutput(1, TList::Class());
  71 }
  72
  73 //______________________________________________
  74 AliAnalysisTaskCounter::AliAnalysisTaskCounter()
  75   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskCounter"),
  76     fAcceptFastCluster(kTRUE),
  77     fZVertexCut(10.),
  78     fTrackMultEtaCut(0.8),
  79     fCaloFilterPatch(kFALSE),
  80     fOutputContainer(0x0),
  81     fESDtrackCuts(AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()),
  82     fTriggerAnalysis (new AliTriggerAnalysis),
  83     fhNEvents(0),
  84     fhXVertex(0),    fhYVertex(0),    fhZVertex(0),
  85     fhXGoodVertex(0),fhYGoodVertex(0),fhZGoodVertex(0),
  86     fhCentrality(0), fhEventPlaneAngle(0)
  87 {
  88   // ctor
  89   DefineOutput(1, TList::Class());
  90 }
  91
  92 //__________________________________________________
  93 AliAnalysisTaskCounter::~AliAnalysisTaskCounter()
  94 {
  95   //Destructor
  96
  97   if (AliAnalysisManager::GetAnalysisManager()->IsProofMode()) return;
  98
  99   if(fOutputContainer)
 100   {
 101     fOutputContainer->Delete() ;
 102     delete fOutputContainer ;
 103   }
 104
 105   if(fESDtrackCuts)    delete fESDtrackCuts;
 106   if(fTriggerAnalysis) delete fTriggerAnalysis;
 107
 108 }
 109
 110
 111 //____________________________________________________
 112 void AliAnalysisTaskCounter::UserCreateOutputObjects()
 113 {
 114   // Init histograms
 115
 116   fOutputContainer = new TList();
 117
 118   fhZVertex     = new TH1F("hZVertex", " Z vertex distribution"   , 200 , -50 , 50  ) ;
 119   fhZVertex->SetXTitle("v_{z} (cm)");
 120   fOutputContainer->Add(fhZVertex);
 121
 122   fhZGoodVertex     = new TH1F("hZGoodVertex", " Good Z vertex distribution"   , 200 , -50 , 50  ) ;
 123   fhZGoodVertex->SetXTitle("v_{z} (cm)");
 124   fOutputContainer->Add(fhZGoodVertex);
 125
 126   fhXVertex     = new TH1F("hXVertex", " X vertex distribution"   , 200 , -2 , 2  ) ;
 127   fhXVertex->SetXTitle("v_{x} (cm)");
 128   fOutputContainer->Add(fhXVertex);
 129
 130   fhXGoodVertex     = new TH1F("hXGoodVertex", " Good X vertex distribution"   , 200 , -2 , 2  ) ;
 131   fhXGoodVertex->SetXTitle("v_{x} (cm)");
 132   fOutputContainer->Add(fhXGoodVertex);
 133
 134   fhYVertex     = new TH1F("hYVertex", " Y vertex distribution"   , 200 , -2 , 2  ) ;
 135   fhYVertex->SetXTitle("v_{y} (cm)");
 136   fOutputContainer->Add(fhYVertex);
 137
 138   fhYGoodVertex     = new TH1F("hYGoodVertex", " Good Y vertex distribution"   , 200 , -2 , 2  ) ;
 139   fhYGoodVertex->SetXTitle("v_{y} (cm)");
 140   fOutputContainer->Add(fhYGoodVertex);
 141
 142   fhCentrality   = new TH1F("hCentrality","Number of events in centrality bin, |vz|<10 cm, method <V0M> ",100,0.,100.) ;
 143   fhCentrality->SetXTitle("Centrality bin");
 144   fOutputContainer->Add(fhCentrality) ;
 145
 146   fhEventPlaneAngle=new TH1F("hEventPlaneAngle","Number of events in event plane, |vz|<10 cm, method <V0> ",100,0.,TMath::Pi()) ;
 147   fhEventPlaneAngle->SetXTitle("EP angle (rad)");
 148   fOutputContainer->Add(fhEventPlaneAngle) ;
 149
 150   fhNEvents = new TH1I("hNEvents", "Number of analyzed events", 21, 0, 21) ;
 151   fhNEvents->SetXTitle("Selection");
 152   fhNEvents->SetYTitle("# events");
 153   fhNEvents->GetXaxis()->SetBinLabel(1 ,"1  = PS");
 154   fhNEvents->GetXaxis()->SetBinLabel(2 ,"2  = 1  & ESD");
 155   fhNEvents->GetXaxis()->SetBinLabel(3 ,"3  = 2  & |Z|<10");
 156   fhNEvents->GetXaxis()->SetBinLabel(4 ,"4  = 2  & |eta|<0.8");
 157   fhNEvents->GetXaxis()->SetBinLabel(5 ,"5  = 3  & 4");
 158   fhNEvents->GetXaxis()->SetBinLabel(6 ,"6  = 2  & V0AND");
 159   fhNEvents->GetXaxis()->SetBinLabel(7 ,"7  = 3  & 6");
 160   fhNEvents->GetXaxis()->SetBinLabel(8 ,"8  = 4  & 6");
 161   fhNEvents->GetXaxis()->SetBinLabel(9 ,"9  = 5  & 6");
 162   fhNEvents->GetXaxis()->SetBinLabel(10,"10 = 2  & not pileup");
 163   fhNEvents->GetXaxis()->SetBinLabel(11,"11 = 2  & good vertex");
 164   fhNEvents->GetXaxis()->SetBinLabel(12,"12 = 3  & 11");
 165   fhNEvents->GetXaxis()->SetBinLabel(13,"13 = 4  & 11");
 166   fhNEvents->GetXaxis()->SetBinLabel(14,"14 = 6  & 11");
 167   fhNEvents->GetXaxis()->SetBinLabel(15,"15 = 9  & 11");
 168   fhNEvents->GetXaxis()->SetBinLabel(16,"16 = 10 & 11");
 169   fhNEvents->GetXaxis()->SetBinLabel(17,"17 = 6  & 10");
 170   fhNEvents->GetXaxis()->SetBinLabel(17,"17 = 1  & |Z|<50");
 171   fhNEvents->GetXaxis()->SetBinLabel(18,"18 = Reject EMCAL 1");
 172   fhNEvents->GetXaxis()->SetBinLabel(19,"19 = 18 & 2");
 173   fhNEvents->GetXaxis()->SetBinLabel(20,"20 = Reject EMCAL 2");
 174   fhNEvents->GetXaxis()->SetBinLabel(21,"20 = 20 & 2");
 175
 176   fOutputContainer->Add(fhNEvents);
 177
 178   fOutputContainer->SetOwner(kTRUE);
 179
 180   PostData(1,fOutputContainer);
 181
 182 }
 183
 184 //_______________________________________________
 185 void AliAnalysisTaskCounter::UserExec(Option_t *)
 186 {
 187   // Main loop
 188   // Called for each event
 189
 190   //printf("___ Event __ %d __\n",(Int_t)Entry());
 191
 192   fhNEvents->Fill(0.5);
 193
 194   AliVEvent * event = InputEvent();
 195   if (!event)
 196   {
 197     printf("AliAnalysisTaskCounter::UserExec() - ERROR: event not available \n");
 198     return;
 199   }
 200
 201   AliESDEvent * esdevent = dynamic_cast<AliESDEvent*> (event);
 202   AliAODEvent * aodevent = dynamic_cast<AliAODEvent*> (event);
 203
 204   TString triggerclasses = "";
 205   if(esdevent) triggerclasses = esdevent->GetFiredTriggerClasses();
 206   if(aodevent) triggerclasses = aodevent->GetFiredTriggerClasses();
 207
 208   if (triggerclasses.Contains("FAST") && !triggerclasses.Contains("ALL") && !fAcceptFastCluster)
 209   {
 210     //printf("Do not count events from fast cluster, trigger name %s\n",triggerclasses.Data());
 211     return;
 212   }
 213
 214   fhNEvents->Fill(1.5);
 215
 216   //Initialize bools
 217   Bool_t bSelectVZ    = kFALSE;
 218   Bool_t bV0AND       = kFALSE;
 219   Bool_t bPileup      = kFALSE;
 220   Bool_t bGoodV       = kFALSE;
 221   Bool_t bSelectTrack = kFALSE;
 222   Int_t  trackMult    = 0;
 223
 224   //---------------------------------
 225   //Get the primary vertex, cut on Z
 226   //---------------------------------
 227   Double_t v[3];
 228   event->GetPrimaryVertex()->GetXYZ(v) ;
 229   fhXVertex->Fill(v[0]);
 230   fhYVertex->Fill(v[1]);
 231   fhZVertex->Fill(v[2]);
 232
 233   if(TMath::Abs(v[2]) < fZVertexCut)
 234   {
 235     bSelectVZ=kTRUE;
 236     fhNEvents->Fill(2.5);
 237   }
 238   //else printf("Vertex out %f \n",v[2]);
 239
 240
 241   //--------------------------------------------------
 242   //Tweak for calorimeter only productions
 243   //--------------------------------------------------
 244   if(fCaloFilterPatch && !esdevent)
 245   {
 246     if(event->GetNumberOfCaloClusters() > 0)
 247     {
 248       AliVCluster * calo = event->GetCaloCluster(0);
 249       if(calo->GetNLabels() == 4){
 250         Int_t * selection = calo->GetLabels();
 251         bPileup   = selection[0];
 252         bGoodV    = selection[1];
 253         bV0AND    = selection[2];
 254         trackMult = selection[3];
 255         //if(selection[0] || selection[1] || selection[2])
 256         //printf(" pu %d, gv %d, v0 %d, track mult %d\n ", selection[0], selection[1], selection[2], selection[3]);
 257         if(trackMult > 0 )
 258           bSelectTrack = kFALSE;
 259       }
 260       else
 261       {
 262         //First filtered AODs, track multiplicity stored there.
 263         trackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
 264       }
 265     }
 266     else
 267     {   //at least one cluster
 268
 269         //First filtered AODs, track multiplicity stored there.
 270         trackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
 271     }
 272   }
 273   else
 274   {
 275     //--------------------------------------------------
 276     //Count tracks, cut on number of tracks in eta < 0.8
 277     //--------------------------------------------------
 278     Int_t nTracks   = event->GetNumberOfTracks() ;
 279     for (Int_t itrack =  0; itrack <  nTracks; itrack++)
 280     {////////////// track loop
 281       AliVTrack * track = (AliVTrack*)event->GetTrack(itrack) ; // retrieve track from esd
 282
 283       //Only for ESDs
 284       if(esdevent && !fESDtrackCuts->AcceptTrack((AliESDtrack*)track)) continue;
 285
 286       //Do not count tracks out of acceptance cut
 287       if(TMath::Abs(track->Eta())< fTrackMultEtaCut) trackMult++;
 288     }
 289   }
 290
 291   //printf("AliAnalysisTaskCounter::UserExec() - Track Mult %d \n",trackMult);
 292
 293   //--------------------------------------------------
 294   // At least one track
 295   //--------------------------------------------------
 296   if (trackMult > 0)
 297   {
 298     bSelectTrack = kTRUE;
 299                   fhNEvents->Fill(3.5);
 300     if(bSelectVZ) fhNEvents->Fill(4.5);
 301   }
 302
 303   //---------------------------------
 304   // V0AND
 305   //---------------------------------
 306   if(esdevent && !fCaloFilterPatch) bV0AND = fTriggerAnalysis->IsOfflineTriggerFired(esdevent, AliTriggerAnalysis::kV0AND);
 307   //else if(aodevent  && !fCaloFilterPatch) bV0AND = //FIXME FOR AODs
 308
 309   if(bV0AND)
 310   {
 311                                     fhNEvents->Fill(5.5);
 312     if (bSelectVZ)                  fhNEvents->Fill(6.5);
 313     if (bSelectTrack)               fhNEvents->Fill(7.5);
 314     if (bSelectVZ &&  bSelectTrack) fhNEvents->Fill(8.5);
 315   }
 316
 317   //---------------------------------
 318   // Pileup
 319   //---------------------------------
 320   if(!fCaloFilterPatch)
 321     bPileup = event->IsPileupFromSPD(3, 0.8, 3., 2., 5.); //Default values, if not it does not compile
 322   //bPileup = event->IsPileupFromSPD();
 323
 324   if (!bPileup)
 325   {
 326                 fhNEvents->Fill(9.5);
 327     if(bV0AND)  fhNEvents->Fill(16.5);
 328   }
 329
 330   //---------------------------------
 331   // Good vertex
 332   //---------------------------------
 333   if(!fCaloFilterPatch) bGoodV = CheckForPrimaryVertex();
 334
 335   //Remove events with  vertex (0,0,0), bad vertex reconstruction
 336   if(TMath::Abs(v[0]) < 1.e-6 &&
 337      TMath::Abs(v[1]) < 1.e-6 &&
 338      TMath::Abs(v[2]) < 1.e-6) bGoodV = kFALSE;
 339
 340   if(bGoodV)
 341   {
 342     fhXGoodVertex->Fill(v[0]);
 343     fhYGoodVertex->Fill(v[1]);
 344     fhZGoodVertex->Fill(v[2]);
 345
 346                      fhNEvents->Fill(10.5);
 347     if(bSelectVZ)    fhNEvents->Fill(11.5);
 348     if(bSelectTrack) fhNEvents->Fill(12.5);
 349     if(bV0AND)       fhNEvents->Fill(13.5);
 350     if(bSelectVZ && bSelectTrack && bV0AND)
 351                      fhNEvents->Fill(14.5);
 352     if(!bPileup)     fhNEvents->Fill(15.5);
 353
 354     if(TMath::Abs(v[2]) < 10.)
 355     {
 356       if(InputEvent()->GetCentrality())
 357         fhCentrality->Fill(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"));
 358
 359       if(InputEvent()->GetEventplane())
 360       {
 361         Float_t ep = InputEvent()->GetEventplane()->GetEventplane("V0", InputEvent());
 362
 363         ep+=TMath::Pi()/2.; // put same range as for <Q> method, [0,pi]
 364
 365         fhEventPlaneAngle->Fill(ep);
 366       }
 367     }
 368
 369   }
 370
 371   //printf("AliAnalysisTaskCounter::UserExec() : z vertex %d, good vertex %d, v0and %d, pile up %d, track mult %d\n ", bSelectVZ, bGoodV, bV0AND, bPileup, trackMult);
 372
 373   // Events that could be rejected in EMCAL
 374   // LHC11a, SM4 and some SM3 events cut with this
 375   Bool_t bEMCALRejected = kFALSE;
 376   for (Int_t i = 0; i < InputEvent()->GetNumberOfCaloClusters(); i++)
 377   {
 378     AliVCluster *clus = InputEvent()->GetCaloCluster(i);
 379     if(clus->IsEMCAL()){
 380       if ((clus->E() > 500 && clus->GetNCells() > 200 ) || clus->GetNCells() > 200)
 381       {
 382
 383         //printf("Counter: Reject event with cluster: E %f, ncells %d\n",clus->E(),clus->GetNCells());
 384
 385                          fhNEvents->Fill(17.5);
 386         if(bSelectVZ)    fhNEvents->Fill(18.5);
 387         bEMCALRejected = kTRUE;
 388         break;
 389       }
 390     }
 391   }
 392
 393   //LHC11a, 3 last runs, cut with this
 394   if(!bEMCALRejected)
 395   {
 396     // Count number of cells in SM3 with energy larger than 0.1, cut on this number
 397     Int_t ncellsSM3 = 0;
 398     Int_t ncellsSM4 = 0;
 399     for(Int_t icell = 0; icell < event->GetEMCALCells()->GetNumberOfCells(); icell++)
 400     {
 401       if(event->GetEMCALCells()->GetAmplitude(icell) > 0.1 && event->GetEMCALCells()->GetCellNumber(icell)/(24*48)==3) ncellsSM3++;
 402       if(event->GetEMCALCells()->GetAmplitude(icell) > 0.1 && event->GetEMCALCells()->GetCellNumber(icell)/(24*48)==4) ncellsSM4++;
 403     }
 404
 405     Int_t ncellcut = 21;
 406     if(triggerclasses.Contains("EMC")) ncellcut = 35;
 407
 408     if( ncellsSM3 >= ncellcut || ncellsSM4 >= 100 )
 409     {
 410       //printf("Counter: reject event with ncells in SM3: ncells %d\n",ncells);
 411
 412                        fhNEvents->Fill(19.5);
 413       if(bSelectVZ)    fhNEvents->Fill(20.5);
 414     }
 415
 416   }
 417
 418   PostData(1,fOutputContainer);
 419
 420 }
 421
 422 //____________________________________________________
 423 Bool_t AliAnalysisTaskCounter::CheckForPrimaryVertex()
 424 {
 425   //Check if the vertex was well reconstructed, copy from V0Reader of conversion group
 426   //It only works for ESDs
 427
 428   AliESDEvent * event = dynamic_cast<AliESDEvent*> (InputEvent());
 429   if(!event) return 1;
 430
 431   if(event->GetPrimaryVertexTracks()->GetNContributors() > 0)
 432   {
 433     return 1;
 434   }
 435
 436   if(event->GetPrimaryVertexTracks()->GetNContributors() < 1)
 437   {
 438     // SPD vertex
 439     if(event->GetPrimaryVertexSPD()->GetNContributors() > 0)
 440     {
 441       //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
 442       return 1;
 443
 444     }
 445     if(event->GetPrimaryVertexSPD()->GetNContributors() < 1)
 446     {
 447       //      cout<<"bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
 448       return 0;
 449     }
 450   }
 451
 452   return 0;
 453   //return fInputEvent->GetPrimaryVertex()->GetNContributors()>0;
 454 }
 455
 456
 457
 458 //_____________________________________________________
 459 void AliAnalysisTaskCounter::FinishTaskOutput()
 460 {
 461   // Put in the output some event summary histograms
 462
 463   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 464   AliInputEventHandler *inputH = dynamic_cast<AliInputEventHandler*>(am->GetInputEventHandler());
 465   if (!inputH) return;
 466   TH2F *histStat = dynamic_cast<TH2F*>(inputH->GetStatistics());
 467   TH2F *histBin0 = dynamic_cast<TH2F*>(inputH->GetStatistics("BIN0"));
 468
 469   if(histStat)
 470     fOutputContainer->Add(histStat);
 471   else if(DebugLevel() > 1)
 472     printf("AliAnalysisTaskCounter::FinishTaskOutput() - Stat histogram not available check, \n if ESDs, that AliPhysicsSelection was on, \n if AODs, if EventStat_temp.root exists \n");
 473
 474   if(histBin0)
 475     fOutputContainer->Add(histBin0);
 476
 477 }