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)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       }
 406       if(cand&&esdEventSelected){
 407         fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexIn);
 408         fh2ESDTriggerCount->Fill(it,kSelected);
 409         AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
 410       }
 411     }
 412   }
 413
 414   if(aod){
 415     const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
 416     //      Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print();
 417     TString vtxTitle(vtxAOD->GetTitle());
 418     for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
 419       Bool_t aodTrig = kFALSE;
 420       aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
 421       Bool_t cand = (aod->GetHeader()->GetOfflineTrigger())&(AliVEvent::kMB);
 422       if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
 423       if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
 424         if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredVertex);
 425         Float_t zvtx = vtxAOD->GetZ();
 426         if(aodTrig&&aodEventSelected){
 427           fh2TriggerVtx->Fill(it,zvtx);
 428           fh2TriggerCount->Fill(it,kTriggeredVertexIn);
 429         }
 430       }
 431     }
 432   }
 433
 434   if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
 435
 436
 437   Double_t ptHard = 0;
 438   Double_t nTrials = 1; // Trials for MC trigger
 439
 440   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 441   AliMCEvent* mcEvent = MCEvent();
 442   //    AliStack *pStack = 0;
 443   if(mcEvent){
 444     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 445     if(pythiaGenHeader){
 446       nTrials = pythiaGenHeader->Trials();
 447       ptHard  = pythiaGenHeader->GetPtHard();
 448       int iProcessType = pythiaGenHeader->ProcessType();
 449       // 11 f+f -> f+f
 450       // 12 f+barf -> f+barf
 451       // 13 f+barf -> g+g
 452       // 28 f+g -> f+g
 453       // 53 g+g -> f+barf
 454       // 68 g+g -> g+g
 455       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 456       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 457       fh1PtHard->Fill(ptHard);
 458       fh1PtHardTrials->Fill(ptHard,nTrials);
 459
 460     }// if pythia gen header
 461   }
 462
 463   // trigger selection
 464
 465
 466   PostData(1, fHistList);
 467 }
 468
 469 Bool_t AliAnalysisTaskJetServices::IsEventSelectedESD(AliESDEvent* esd){
 470   if(!esd)return kFALSE;
 471   const AliESDVertex *vtx = esd->GetPrimaryVertex();
 472   // We can check the type of the vertex by its name and title
 473   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
 474   // vtx       name                  title
 475   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
 476   // TPC       TPCVertex             VertexerTracksNoConstraint
 477   // SPD       SPDVertex             vertexer: 3D
 478   // SPD       SPDVertex             vertexer: Z
 479
 480
 481
 482   if(vtx->GetNContributors()<3)return kFALSE;
 483
 484   // do not want tpc only primary vertex
 485   TString vtxName(vtx->GetName());
 486   if(vtxName.Contains("TPCVertex"))return kFALSE;
 487   Float_t zvtx = vtx->GetZ();
 488   Float_t yvtx = vtx->GetY();
 489   Float_t xvtx = vtx->GetX();
 490
 491   // here we may fill the histos for run dependence....
 492
 493   xvtx -= fVtxXMean;
 494   yvtx -= fVtxYMean;
 495   zvtx -= fVtxZMean;
 496
 497   Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 498
 499   Bool_t eventSel = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
 500   return eventSel;
 501 }
 502
 503 Bool_t AliAnalysisTaskJetServices::IsEventPileUpESD(AliESDEvent* esd){
 504   if(!esd)return kFALSE;
 505   return esd->IsPileupFromSPD();
 506 }
 507
 508 Bool_t AliAnalysisTaskJetServices::IsEventCosmicESD(AliESDEvent* esd){
 509   if(!esd)return kFALSE;
 510   // add track cuts for which we look for cosmics...
 511
 512   Bool_t isCosmic = kFALSE;
 513   Int_t nTracks = esd->GetNumberOfTracks();
 514   Int_t nCosmicCandidates = 0;
 515
 516   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
 517     AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
 518     if (!track1)  continue;
 519     UInt_t status1 = track1->GetStatus();
 520     //If track is ITS stand alone track, skip the track
 521     if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue;
 522     if(track1->Pt()<fPtMinCosmic) continue;
 523     //Start 2nd track loop to look for correlations
 524     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
 525       AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
 526       if(!track2) continue;
 527       UInt_t status2 = track2->GetStatus();
 528       //If track is ITS stand alone track, skip the track
 529       if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue;
 530       if(track2->Pt()<fPtMinCosmic) continue;
 531       //Check if back-to-back
 532       Double_t mom1[3],mom2[3];
 533       track1->GetPxPyPz(mom1);
 534       track2->GetPxPyPz(mom2);
 535       TVector3 momv1(mom1[0],mom1[1],mom1[2]);
 536       TVector3 momv2(mom2[0],mom2[1],mom2[2]);
 537       Float_t theta = (float)(momv1.Phi()-momv2.Phi());
 538       if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
 539
 540       Float_t deltaPhi = track1->Phi()-track2->Phi();
 541       if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
 542
 543       Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
 544       if(rIsol<fRIsolMinCosmic) continue;
 545
 546       if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
 547         nCosmicCandidates+=1;
 548         isCosmic = kTRUE;
 549       }
 550
 551     }
 552   }
 553
 554   fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
 555
 556   return isCosmic;
 557 }
 558
 559
 560 Bool_t AliAnalysisTaskJetServices::IsEventSelectedAOD(AliAODEvent* aod){
 561   if(!aod)return kFALSE;
 562   const AliAODVertex *vtx = aod->GetPrimaryVertex();
 563   // We can check the type of the vertex by its name and title
 564   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
 565   // vtx       name                  title
 566   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
 567   // TPC       TPCVertex             VertexerTracksNoConstraint
 568   // SPD       SPDVertex             vertexer: 3D
 569   // SPD       SPDVertex             vertexer: Z
 570
 571
 572
 573   if(vtx->GetNContributors()<3)return kFALSE;
 574
 575   // do not want tpc only primary vertex
 576   TString vtxName(vtx->GetName());
 577   if(vtxName.Contains("TPCVertex"))return kFALSE;
 578
 579   Float_t zvtx = vtx->GetZ();
 580   Float_t yvtx = vtx->GetY();
 581   Float_t xvtx = vtx->GetX();
 582
 583   // here we may fill the histos for run dependence....
 584
 585   xvtx -= fVtxXMean;
 586   yvtx -= fVtxYMean;
 587   zvtx -= fVtxZMean;
 588
 589   Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 590   Bool_t eventSel = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
 591   return eventSel;
 592 }
 593
 594
 595
 596
 597 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
 598 {
 599   // Terminate analysis
 600   //
 601 }