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