PWG4/JetTasks/AliAnalysisTaskJetServices.cxx

   1 // **************************************
   2 // Task used for the correction of determiantion of reconstructed jet spectra
   3 // Compares input (gen) and output (rec) jets
   4 // *******************************************
   5
   6
   7 /**************************************************************************
   8  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   9  *                                                                        *
  10  * Author: The ALICE Off-line Project.                                    *
  11  * Contributors are mentioned in the code where appropriate.              *
  12  *                                                                        *
  13  * Permission to use, copy, modify and distribute this software and its   *
  14  * documentation strictly for non-commercial purposes is hereby granted   *
  15  * without fee, provided that the above copyright notice appears in all   *
  16  * copies and that both the copyright notice and this permission notice   *
  17  * appear in the supporting documentation. The authors make no claims     *
  18  * about the suitability of this software for any purpose. It is          *
  19  * provided "as is" without express or implied warranty.                  *
  20  **************************************************************************/
  21
  22
  23 #include <TROOT.h>
  24 #include <TRandom.h>
  25 #include <TString.h>
  26 #include <TSystem.h>
  27 #include <TInterpreter.h>
  28 #include <TChain.h>
  29 #include <TFile.h>
  30 #include <TKey.h>
  31 #include <TH1F.h>
  32 #include <TH2F.h>
  33 #include <TH3F.h>
  34 #include <TProfile.h>
  35 #include <TList.h>
  36 #include <TLorentzVector.h>
  37 #include <TClonesArray.h>
  38 #include "TDatabasePDG.h"
  39
  40 #include "AliAnalysisTaskJetServices.h"
  41 #include "AliAnalysisDataContainer.h"
  42 #include "AliAnalysisDataSlot.h"
  43 #include "AliAnalysisManager.h"
  44 #include "AliJetFinder.h"
  45 #include "AliJetHeader.h"
  46 #include "AliJetReader.h"
  47 #include "AliJetReaderHeader.h"
  48 #include "AliUA1JetHeaderV1.h"
  49 #include "AliESDEvent.h"
  50 #include "AliAODEvent.h"
  51 #include "AliAODHandler.h"
  52 #include "AliAODTrack.h"
  53 #include "AliAODJet.h"
  54 #include "AliAODMCParticle.h"
  55 #include "AliMCEventHandler.h"
  56 #include "AliMCEvent.h"
  57 #include "AliStack.h"
  58 #include "AliGenPythiaEventHeader.h"
  59 #include "AliJetKineReaderHeader.h"
  60 #include "AliGenCocktailEventHeader.h"
  61 #include "AliInputEventHandler.h"
  62 #include "AliPhysicsSelection.h"
  63
  64
  65 #include "AliAnalysisHelperJetTasks.h"
  66
  67 ClassImp(AliAnalysisTaskJetServices)
  68
  69 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(): AliAnalysisTaskSE(),
  70   fUseAODInput(kFALSE),
  71   fUsePhysicsSelection(kFALSE),
  72   fRealData(kFALSE),
  73   fSelectionInfoESD(0),
  74   fAvgTrials(1),
  75   fVtxXMean(0),
  76   fVtxYMean(0),
  77   fVtxZMean(0),
  78   fVtxRCut(1.),
  79   fVtxZCut(8.),
  80   fPtMinCosmic(5.),
  81   fRIsolMinCosmic(3.),
  82   fMaxCosmicAngle(0.01),
  83   fh1Xsec(0x0),
  84   fh1Trials(0x0),
  85   fh1PtHard(0x0),
  86   fh1PtHardTrials(0x0),
  87   fh1SelectionInfoESD(0x0),
  88   fh2TriggerCount(0x0),
  89   fh2ESDTriggerCount(0x0),
  90   fh2TriggerVtx(0x0),
  91   fh2ESDTriggerVtx(0x0),
  92   fh2ESDTriggerRun(0x0),
  93   fh2VtxXY(0x0),
  94   fh1NCosmicsPerEvent(0x0),
  95   fHistList(0x0)
  96 {
  97   fRunRange[0] = fRunRange[1] = 0;
  98 }
  99
 100 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
 101   AliAnalysisTaskSE(name),
 102   fUseAODInput(kFALSE),
 103   fUsePhysicsSelection(kFALSE),
 104   fRealData(kFALSE),
 105   fSelectionInfoESD(0),
 106   fAvgTrials(1),
 107   fVtxXMean(0),
 108   fVtxYMean(0),
 109   fVtxZMean(0),
 110   fVtxRCut(1.),
 111   fVtxZCut(8.),
 112   fPtMinCosmic(5.),
 113   fRIsolMinCosmic(3.),
 114   fMaxCosmicAngle(0.01),
 115   fh1Xsec(0x0),
 116   fh1Trials(0x0),
 117   fh1PtHard(0x0),
 118   fh1PtHardTrials(0x0),
 119   fh1SelectionInfoESD(0x0),
 120   fh2TriggerCount(0x0),
 121   fh2ESDTriggerCount(0x0),
 122   fh2TriggerVtx(0x0),
 123   fh2ESDTriggerVtx(0x0),
 124   fh2ESDTriggerRun(0x0),
 125   fh2VtxXY(0x0),
 126   fh1NCosmicsPerEvent(0x0),
 127   fHistList(0x0)
 128 {
 129   fRunRange[0] = fRunRange[1] = 0;
 130   DefineOutput(1,TList::Class());
 131 }
 132
 133
 134
 135 Bool_t AliAnalysisTaskJetServices::Notify()
 136 {
 137   //
 138   // Implemented Notify() to read the cross sections
 139   // and number of trials from pyxsec.root
 140   //
 141
 142   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 143   Float_t xsection = 0;
 144   Float_t ftrials  = 1;
 145
 146   fAvgTrials = 1;
 147   if(tree){
 148     TFile *curfile = tree->GetCurrentFile();
 149     if (!curfile) {
 150       Error("Notify","No current file");
 151       return kFALSE;
 152     }
 153     if(!fh1Xsec||!fh1Trials){
 154       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 155       return kFALSE;
 156     }
 157     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 158     fh1Xsec->Fill("<#sigma>",xsection);
 159     // construct a poor man average trials
 160     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 161     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 162   }
 163   return kTRUE;
 164 }
 165
 166 void AliAnalysisTaskJetServices::UserCreateOutputObjects()
 167 {
 168
 169   //
 170   // Create the output container
 171   //
 172
 173
 174   if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
 175
 176   OpenFile(1);
 177   if(!fHistList)fHistList = new TList();
 178
 179   Bool_t oldStatus = TH1::AddDirectoryStatus();
 180   TH1::AddDirectory(kFALSE);
 181   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 182   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 183   fHistList->Add(fh1Xsec);
 184
 185   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 186   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 187   fHistList->Add(fh1Trials);
 188
 189   const Int_t nBinPt = 100;
 190   Double_t binLimitsPt[nBinPt+1];
 191   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 192     if(iPt == 0){
 193       binLimitsPt[iPt] = 0.0;
 194     }
 195     else {// 1.0
 196       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
 197     }
 198   }
 199
 200   fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
 201   fHistList->Add(fh2TriggerCount);
 202
 203   fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
 204   fHistList->Add(fh2ESDTriggerCount);
 205
 206   fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
 207   fHistList->Add(fh2TriggerVtx);
 208
 209   fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
 210   fHistList->Add(fh2ESDTriggerVtx);
 211
 212
 213   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 214   fHistList->Add(fh1PtHard);
 215   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 216   fHistList->Add(fh1PtHardTrials);
 217   fh1SelectionInfoESD = new TH1F("fh1SelectionInfoESD","Bit Masks that satisfy the selection info",
 218                                  AliAnalysisHelperJetTasks::kTotalSelections,0.5,AliAnalysisHelperJetTasks::kTotalSelections+0.5);
 219
 220   fHistList->Add(fh1SelectionInfoESD);
 221   // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range
 222   // 3 triggers BB BE/EB EE
 223
 224   fh2ESDTriggerRun = new TH2F("fh2ESDTriggerRun","Trigger vs run number:run;trigger",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5,10,-0.5,9.5);
 225   fHistList->Add(fh2ESDTriggerRun);
 226
 227   fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
 228   fHistList->Add(fh2VtxXY);
 229   // =========== Switch on Sumw2 for all histos ===========
 230   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 231     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 232     if (h1){
 233       h1->Sumw2();
 234       continue;
 235     }
 236     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 237     if(hn)hn->Sumw2();
 238   }
 239
 240   fh1NCosmicsPerEvent = new TH1F("fh1NCosmicsPerEvent","Number of cosmic candidates per event",10,0.,10.);
 241   fHistList->Add(fh1NCosmicsPerEvent),
 242
 243
 244   TH1::AddDirectory(oldStatus);
 245 }
 246
 247 void AliAnalysisTaskJetServices::Init()
 248 {
 249   //
 250   // Initialization
 251   //
 252   if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
 253
 254 }
 255
 256 void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
 257 {
 258
 259   //
 260   // Execute analysis for current event
 261   //
 262
 263   AliAODEvent *aod = 0;
 264   AliESDEvent *esd = 0;
 265
 266   AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
 267   fSelectionInfoESD = 0; // reset
 268   AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
 269
 270
 271   if(fUseAODInput){
 272     aod = dynamic_cast<AliAODEvent*>(InputEvent());
 273     if(!aod){
 274       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
 275       return;
 276     }
 277     // fethc the header
 278   }
 279   else{
 280     //  assume that the AOD is in the general output...
 281     aod  = AODEvent();
 282     if(!aod){
 283       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 284       return;
 285     }
 286     esd = dynamic_cast<AliESDEvent*>(InputEvent());
 287   }
 288
 289   fSelectionInfoESD |= AliAnalysisHelperJetTasks::kNone;
 290
 291   if(esd){
 292     Float_t run = (Float_t)esd->GetRunNumber();
 293     const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
 294     Float_t zvtx = -999;
 295     Float_t xvtx = -999;
 296     Float_t yvtx = -999;
 297
 298     if(vtxESD->GetNContributors()>2){
 299       zvtx = vtxESD->GetZ();
 300       yvtx = vtxESD->GetY();
 301       xvtx = vtxESD->GetX();
 302     }
 303
 304     Int_t iTrig = -1;
 305     if(esd->GetFiredTriggerClasses().Contains("CINT1B")
 306        ||esd->GetFiredTriggerClasses().Contains("CSMBB")
 307        ||esd->GetFiredTriggerClasses().Contains("MB1")
 308        ||esd->GetFiredTriggerClasses().Contains("CINT6B")){
 309       iTrig = 0;
 310       fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kBunchBunch;
 311     }
 312     else if(esd->GetFiredTriggerClasses().Contains("CINT1A")
 313             ||esd->GetFiredTriggerClasses().Contains("CSMBA")
 314             ||esd->GetFiredTriggerClasses().Contains("CINT6A")
 315             ||esd->GetFiredTriggerClasses().Contains("CINT1C")
 316             ||esd->GetFiredTriggerClasses().Contains("CSMBC")
 317             ||esd->GetFiredTriggerClasses().Contains("CINT6C")){
 318       // empty bunch or bunch empty
 319       iTrig = 1;
 320       fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kBunchEmpty;
 321     }
 322     else if(esd->GetFiredTriggerClasses().Contains("CINT1-E")
 323        ||esd->GetFiredTriggerClasses().Contains("CINT6-E")){
 324       iTrig = 2;
 325       fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kEmptyEmpty;
 326     }
 327
 328
 329     if(iTrig>=0){
 330       iTrig *= 3;
 331       fh2ESDTriggerRun->Fill(run,iTrig+1);
 332       if(vtxESD->GetNContributors()>2){
 333         fh2ESDTriggerRun->Fill(run,iTrig+2);
 334         fh2VtxXY->Fill(xvtx,yvtx);
 335       }
 336       xvtx -= fVtxXMean;
 337       yvtx -= fVtxYMean;
 338       zvtx -= fVtxZMean;
 339       Float_t r2 = xvtx *xvtx + yvtx *yvtx;
 340       if(TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut))fh2ESDTriggerRun->Fill(run,iTrig+3);
 341     }
 342     else{
 343       fh2ESDTriggerRun->Fill(run,0);
 344     }
 345
 346
 347   }
 348
 349
 350   // Apply additional constraints
 351   Bool_t esdEventSelected = IsEventSelectedESD(esd);
 352   Bool_t esdEventPileUp = IsEventPileUpESD(esd);
 353   Bool_t esdEventCosmic = IsEventCosmicESD(esd);
 354
 355   Bool_t aodEventSelected = IsEventSelectedAOD(aod);
 356
 357   Bool_t physicsSelection = ((fInputHandler->IsEventSelected())&AliVEvent::kMB);
 358
 359   if(esdEventSelected) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kVertexIn;
 360   if(esdEventPileUp)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsPileUp;
 361   if(esdEventCosmic)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsCosmic;
 362   if(physicsSelection) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kPhysicsSelection;
 363
 364
 365   // here we have all selection information, fill histogram
 366   for(unsigned int i = 1;i<(UInt_t)fh1SelectionInfoESD->GetNbinsX();i++){
 367     if((i&fSelectionInfoESD)==i)fh1SelectionInfoESD->Fill(i);
 368   }
 369   AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
 370
 371   if(esd&&aod&&fDebug){
 372     if(esdEventSelected&&!aodEventSelected){
 373       Printf("%s:%d Different Selection for ESD and AOD",(char*)__FILE__,__LINE__);
 374       const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
 375       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
 376       Printf("ESD Vtx %s %s %d",vtxESD->GetName(),vtxESD->GetTitle(),vtxESD->GetNContributors());
 377       vtxESD->Print();
 378       Printf("AOD Vtx %s %s %d",vtxAOD->GetName(),vtxAOD->GetTitle(),vtxAOD->GetNContributors());
 379       vtxAOD->Print();
 380     }
 381   }
 382
 383   // loop over all possible triggers for esd
 384
 385   if(esd){
 386     const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
 387     //      Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print();
 388     TString vtxName(vtxESD->GetName());
 389     for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
 390       Bool_t esdTrig = kFALSE;
 391       esdTrig = AliAnalysisHelperJetTasks::IsTriggerFired(esd,(AliAnalysisHelperJetTasks::Trigger)it);
 392       if(esdTrig)fh2ESDTriggerCount->Fill(it,kAllTriggered);
 393       Bool_t cand = physicsSelection;
 394       if(cand){
 395         fh2ESDTriggerCount->Fill(it,kSelectedALICE);
 396       }
 397       if(!fUsePhysicsSelection)cand =  AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
 398       if(vtxESD->GetNContributors()>2&&!vtxName.Contains("TPCVertex")){
 399         if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredVertex);
 400         Float_t zvtx = vtxESD->GetZ();
 401         if(esdEventSelected&&esdTrig){
 402           fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn);
 403           fh2ESDTriggerVtx->Fill(it,zvtx);
 404         }
 405         if(cand)fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexValid);
 406       }
 407       if(cand&&esdEventSelected){
 408         fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexIn);
 409         fh2ESDTriggerCount->Fill(it,kSelected);
 410         AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
 411       }
 412     }
 413   }
 414
 415   if(aod){
 416     const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
 417     //      Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print();
 418     TString vtxTitle(vtxAOD->GetTitle());
 419     for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
 420       Bool_t aodTrig = kFALSE;
 421       aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
 422       Bool_t cand = ((aod->GetHeader()->GetOfflineTrigger())&(AliVEvent::kMB))==AliVEvent::kMB;
 423       if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
 424       if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
 425         if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredVertex);
 426         Float_t zvtx = vtxAOD->GetZ();
 427         if(cand&&aodEventSelected){
 428           fh2TriggerCount->Fill(it,kSelected);
 429         }
 430         if(aodTrig&&aodEventSelected){
 431           fh2TriggerVtx->Fill(it,zvtx);
 432           fh2TriggerCount->Fill(it,kTriggeredVertexIn);
 433         }
 434         if(cand){
 435           fh2TriggerCount->Fill(it,kSelectedALICEVertexValid);
 436           if(aodEventSelected)fh2TriggerCount->Fill(it,kSelectedALICEVertexIn);
 437         }
 438
 439       }
 440     }
 441   }
 442
 443   if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
 444
 445
 446   Double_t ptHard = 0;
 447   Double_t nTrials = 1; // Trials for MC trigger
 448
 449   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 450   AliMCEvent* mcEvent = MCEvent();
 451   //    AliStack *pStack = 0;
 452   if(mcEvent){
 453     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 454     if(pythiaGenHeader){
 455       nTrials = pythiaGenHeader->Trials();
 456       ptHard  = pythiaGenHeader->GetPtHard();
 457       int iProcessType = pythiaGenHeader->ProcessType();
 458       // 11 f+f -> f+f
 459       // 12 f+barf -> f+barf
 460       // 13 f+barf -> g+g
 461       // 28 f+g -> f+g
 462       // 53 g+g -> f+barf
 463       // 68 g+g -> g+g
 464       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 465       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 466       fh1PtHard->Fill(ptHard);
 467       fh1PtHardTrials->Fill(ptHard,nTrials);
 468
 469     }// if pythia gen header
 470   }
 471
 472   // trigger selection
 473
 474
 475   PostData(1, fHistList);
 476 }
 477
 478 Bool_t AliAnalysisTaskJetServices::IsEventSelectedESD(AliESDEvent* esd){
 479   if(!esd)return kFALSE;
 480   const AliESDVertex *vtx = esd->GetPrimaryVertex();
 481   // We can check the type of the vertex by its name and title
 482   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
 483   // vtx       name                  title
 484   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
 485   // TPC       TPCVertex             VertexerTracksNoConstraint
 486   // SPD       SPDVertex             vertexer: 3D
 487   // SPD       SPDVertex             vertexer: Z
 488
 489
 490
 491   if(vtx->GetNContributors()<3)return kFALSE;
 492
 493   // do not want tpc only primary vertex
 494   TString vtxName(vtx->GetName());
 495   if(vtxName.Contains("TPCVertex"))return kFALSE;
 496   Float_t zvtx = vtx->GetZ();
 497   Float_t yvtx = vtx->GetY();
 498   Float_t xvtx = vtx->GetX();
 499
 500   // here we may fill the histos for run dependence....
 501
 502   xvtx -= fVtxXMean;
 503   yvtx -= fVtxYMean;
 504   zvtx -= fVtxZMean;
 505
 506   Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 507
 508   Bool_t eventSel = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
 509   return eventSel;
 510 }
 511
 512 Bool_t AliAnalysisTaskJetServices::IsEventPileUpESD(AliESDEvent* esd){
 513   if(!esd)return kFALSE;
 514   return esd->IsPileupFromSPD();
 515 }
 516
 517 Bool_t AliAnalysisTaskJetServices::IsEventCosmicESD(AliESDEvent* esd){
 518   if(!esd)return kFALSE;
 519   // add track cuts for which we look for cosmics...
 520
 521   Bool_t isCosmic = kFALSE;
 522   Int_t nTracks = esd->GetNumberOfTracks();
 523   Int_t nCosmicCandidates = 0;
 524
 525   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
 526     AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
 527     if (!track1)  continue;
 528     UInt_t status1 = track1->GetStatus();
 529     //If track is ITS stand alone track, skip the track
 530     if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue;
 531     if(track1->Pt()<fPtMinCosmic) continue;
 532     //Start 2nd track loop to look for correlations
 533     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
 534       AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
 535       if(!track2) continue;
 536       UInt_t status2 = track2->GetStatus();
 537       //If track is ITS stand alone track, skip the track
 538       if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue;
 539       if(track2->Pt()<fPtMinCosmic) continue;
 540       //Check if back-to-back
 541       Double_t mom1[3],mom2[3];
 542       track1->GetPxPyPz(mom1);
 543       track2->GetPxPyPz(mom2);
 544       TVector3 momv1(mom1[0],mom1[1],mom1[2]);
 545       TVector3 momv2(mom2[0],mom2[1],mom2[2]);
 546       Float_t theta = (float)(momv1.Phi()-momv2.Phi());
 547       if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
 548
 549       Float_t deltaPhi = track1->Phi()-track2->Phi();
 550       if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
 551
 552       Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
 553       if(rIsol<fRIsolMinCosmic) continue;
 554
 555       if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
 556         nCosmicCandidates+=1;
 557         isCosmic = kTRUE;
 558       }
 559
 560     }
 561   }
 562
 563   fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
 564
 565   return isCosmic;
 566 }
 567
 568
 569 Bool_t AliAnalysisTaskJetServices::IsEventSelectedAOD(AliAODEvent* aod){
 570   if(!aod)return kFALSE;
 571   const AliAODVertex *vtx = aod->GetPrimaryVertex();
 572   // We can check the type of the vertex by its name and title
 573   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
 574   // vtx       name                  title
 575   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
 576   // TPC       TPCVertex             VertexerTracksNoConstraint
 577   // SPD       SPDVertex             vertexer: 3D
 578   // SPD       SPDVertex             vertexer: Z
 579
 580
 581
 582   if(vtx->GetNContributors()<3)return kFALSE;
 583
 584   // do not want tpc only primary vertex
 585   TString vtxName(vtx->GetName());
 586   if(vtxName.Contains("TPCVertex"))return kFALSE;
 587
 588   Float_t zvtx = vtx->GetZ();
 589   Float_t yvtx = vtx->GetY();
 590   Float_t xvtx = vtx->GetX();
 591
 592   // here we may fill the histos for run dependence....
 593
 594   xvtx -= fVtxXMean;
 595   yvtx -= fVtxYMean;
 596   zvtx -= fVtxZMean;
 597
 598   Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 599   Bool_t eventSel = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
 600   return eventSel;
 601 }
 602
 603
 604
 605
 606 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
 607 {
 608   // Terminate analysis
 609   //
 610 }