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   fEventCutInfoESD(0),
  75   fAvgTrials(1),
  76   fVtxXMean(0),
  77   fVtxYMean(0),
  78   fVtxZMean(0),
  79   fVtxRCut(1.),
  80   fVtxZCut(8.),
  81   fPtMinCosmic(5.),
  82   fRIsolMinCosmic(3.),
  83   fMaxCosmicAngle(0.01),
  84   fh1Xsec(0x0),
  85   fh1Trials(0x0),
  86   fh1PtHard(0x0),
  87   fh1PtHardTrials(0x0),
  88   fh1SelectionInfoESD(0x0),
  89   fh1EventCutInfoESD(0),
  90   fh2TriggerCount(0x0),
  91   fh2ESDTriggerCount(0x0),
  92   fh2TriggerVtx(0x0),
  93   fh2ESDTriggerVtx(0x0),
  94   fh2ESDTriggerRun(0x0),
  95   fh2VtxXY(0x0),
  96   fh1NCosmicsPerEvent(0x0),
  97   fHistList(0x0)
  98 {
  99   fRunRange[0] = fRunRange[1] = 0;
 100 }
 101
 102 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
 103   AliAnalysisTaskSE(name),
 104   fUseAODInput(kFALSE),
 105   fUsePhysicsSelection(kFALSE),
 106   fRealData(kFALSE),
 107   fSelectionInfoESD(0),
 108   fEventCutInfoESD(0),
 109   fAvgTrials(1),
 110   fVtxXMean(0),
 111   fVtxYMean(0),
 112   fVtxZMean(0),
 113   fVtxRCut(1.),
 114   fVtxZCut(8.),
 115   fPtMinCosmic(5.),
 116   fRIsolMinCosmic(3.),
 117   fMaxCosmicAngle(0.01),
 118   fh1Xsec(0x0),
 119   fh1Trials(0x0),
 120   fh1PtHard(0x0),
 121   fh1PtHardTrials(0x0),
 122   fh1SelectionInfoESD(0x0),
 123   fh1EventCutInfoESD(0),
 124   fh2TriggerCount(0x0),
 125   fh2ESDTriggerCount(0x0),
 126   fh2TriggerVtx(0x0),
 127   fh2ESDTriggerVtx(0x0),
 128   fh2ESDTriggerRun(0x0),
 129   fh2VtxXY(0x0),
 130   fh1NCosmicsPerEvent(0x0),
 131   fHistList(0x0)
 132 {
 133   fRunRange[0] = fRunRange[1] = 0;
 134   DefineOutput(1,TList::Class());
 135 }
 136
 137
 138
 139 Bool_t AliAnalysisTaskJetServices::Notify()
 140 {
 141   //
 142   // Implemented Notify() to read the cross sections
 143   // and number of trials from pyxsec.root
 144   //
 145
 146   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 147   Float_t xsection = 0;
 148   Float_t ftrials  = 1;
 149
 150   fAvgTrials = 1;
 151   if(tree){
 152     TFile *curfile = tree->GetCurrentFile();
 153     if (!curfile) {
 154       Error("Notify","No current file");
 155       return kFALSE;
 156     }
 157     if(!fh1Xsec||!fh1Trials){
 158       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 159       return kFALSE;
 160     }
 161     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 162     fh1Xsec->Fill("<#sigma>",xsection);
 163     // construct a poor man average trials
 164     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 165     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 166   }
 167   return kTRUE;
 168 }
 169
 170 void AliAnalysisTaskJetServices::UserCreateOutputObjects()
 171 {
 172
 173   //
 174   // Create the output container
 175   //
 176
 177
 178   if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
 179
 180   OpenFile(1);
 181   if(!fHistList)fHistList = new TList();
 182
 183   Bool_t oldStatus = TH1::AddDirectoryStatus();
 184   TH1::AddDirectory(kFALSE);
 185   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 186   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 187   fHistList->Add(fh1Xsec);
 188
 189   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 190   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 191   fHistList->Add(fh1Trials);
 192
 193   const Int_t nBinPt = 100;
 194   Double_t binLimitsPt[nBinPt+1];
 195   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 196     if(iPt == 0){
 197       binLimitsPt[iPt] = 0.0;
 198     }
 199     else {// 1.0
 200       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
 201     }
 202   }
 203
 204   fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
 205   fHistList->Add(fh2TriggerCount);
 206
 207   fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
 208   fHistList->Add(fh2ESDTriggerCount);
 209
 210   fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
 211   fHistList->Add(fh2TriggerVtx);
 212
 213   fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
 214   fHistList->Add(fh2ESDTriggerVtx);
 215
 216
 217   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 218   fHistList->Add(fh1PtHard);
 219   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 220   fHistList->Add(fh1PtHardTrials);
 221   fh1SelectionInfoESD = new TH1F("fh1SelectionInfoESD","Bit Masks that satisfy the selection info",
 222                                  AliAnalysisHelperJetTasks::kTotalSelections,0.5,AliAnalysisHelperJetTasks::kTotalSelections+0.5);
 223   fHistList->Add(fh1SelectionInfoESD);
 224
 225   fh1EventCutInfoESD = new TH1F("fh1EventCutInfoESD","Bit Masks that for the events after each step of cuts",
 226                                 kTotalEventCuts,0.5,kTotalEventCuts+0.5);
 227   fHistList->Add(fh1EventCutInfoESD);
 228
 229   // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range
 230   // 3 triggers BB BE/EB EE
 231
 232   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);
 233   fHistList->Add(fh2ESDTriggerRun);
 234
 235   fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
 236   fHistList->Add(fh2VtxXY);
 237   // =========== Switch on Sumw2 for all histos ===========
 238   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 239     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 240     if (h1){
 241       h1->Sumw2();
 242       continue;
 243     }
 244     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 245     if(hn)hn->Sumw2();
 246   }
 247
 248   fh1NCosmicsPerEvent = new TH1F("fh1NCosmicsPerEvent","Number of cosmic candidates per event",10,0.,10.);
 249   fHistList->Add(fh1NCosmicsPerEvent),
 250
 251
 252   TH1::AddDirectory(oldStatus);
 253 }
 254
 255 void AliAnalysisTaskJetServices::Init()
 256 {
 257   //
 258   // Initialization
 259   //
 260   if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
 261
 262 }
 263
 264 void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
 265 {
 266
 267   //
 268   // Execute analysis for current event
 269   //
 270
 271   AliAODEvent *aod = 0;
 272   AliESDEvent *esd = 0;
 273
 274   AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
 275   fSelectionInfoESD = 0; // reset
 276   fEventCutInfoESD = 0; // reset
 277   AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
 278
 279
 280   if(fUseAODInput){
 281     aod = dynamic_cast<AliAODEvent*>(InputEvent());
 282     if(!aod){
 283       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
 284       return;
 285     }
 286     // fethc the header
 287   }
 288   else{
 289     //  assume that the AOD is in the general output...
 290     aod  = AODEvent();
 291     if(!aod){
 292       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 293       return;
 294     }
 295     esd = dynamic_cast<AliESDEvent*>(InputEvent());
 296   }
 297
 298   fSelectionInfoESD |= AliAnalysisHelperJetTasks::kNone;
 299   fEventCutInfoESD |= AliAnalysisHelperJetTasks::kNone;
 300
 301   Bool_t esdVtxValid = false;
 302   Bool_t esdVtxIn = false;
 303   Bool_t aodVtxValid = false;
 304   Bool_t aodVtxIn = false;
 305
 306   if(esd){
 307     Float_t run = (Float_t)esd->GetRunNumber();
 308     const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
 309     esdVtxValid =
 310     esdVtxIn = IsVertexIn(vtxESD);
 311     Float_t zvtx = -999;
 312     Float_t xvtx = -999;
 313     Float_t yvtx = -999;
 314
 315     if(esdVtxValid){
 316       zvtx = vtxESD->GetZ();
 317       yvtx = vtxESD->GetY();
 318       xvtx = vtxESD->GetX();
 319     }
 320
 321
 322
 323
 324     // CKB this can be cleaned up a bit...
 325
 326     Int_t iTrig = -1;
 327     if(esd->GetFiredTriggerClasses().Contains("CINT1B")
 328        ||esd->GetFiredTriggerClasses().Contains("CSMBB")
 329        ||esd->GetFiredTriggerClasses().Contains("MB1")
 330        ||esd->GetFiredTriggerClasses().Contains("CINT6B")){
 331       iTrig = 0;
 332       fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kBunchBunch;
 333     }
 334     else if(esd->GetFiredTriggerClasses().Contains("CINT1A")
 335             ||esd->GetFiredTriggerClasses().Contains("CSMBA")
 336             ||esd->GetFiredTriggerClasses().Contains("CINT6A")
 337             ||esd->GetFiredTriggerClasses().Contains("CINT1C")
 338             ||esd->GetFiredTriggerClasses().Contains("CSMBC")
 339             ||esd->GetFiredTriggerClasses().Contains("CINT6C")){
 340       // empty bunch or bunch empty
 341       iTrig = 1;
 342       fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kBunchEmpty;
 343     }
 344     else if(esd->GetFiredTriggerClasses().Contains("CINT1-E")
 345        ||esd->GetFiredTriggerClasses().Contains("CINT6-E")){
 346       iTrig = 2;
 347       fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kEmptyEmpty;
 348     }
 349
 350
 351     if(iTrig>=0){
 352       iTrig *= 3;
 353       fh2ESDTriggerRun->Fill(run,iTrig+1);
 354       if(vtxESD->GetNContributors()>2){
 355         fh2ESDTriggerRun->Fill(run,iTrig+2);
 356         fh2VtxXY->Fill(xvtx,yvtx);
 357       }
 358       xvtx -= fVtxXMean;
 359       yvtx -= fVtxYMean;
 360       zvtx -= fVtxZMean;
 361       Float_t r2 = xvtx *xvtx + yvtx *yvtx;
 362       if(TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut))fh2ESDTriggerRun->Fill(run,iTrig+3);
 363     }
 364     else{
 365       fh2ESDTriggerRun->Fill(run,0);
 366     }
 367     // BKC
 368   }
 369
 370
 371   // Apply additional constraints
 372   Bool_t esdEventSelected = IsEventSelected(esd);
 373   Bool_t esdEventPileUp = IsEventPileUp(esd);
 374   Bool_t esdEventCosmic = IsEventCosmic(esd);
 375
 376   Bool_t aodEventSelected = IsEventSelected(aod);
 377
 378   Bool_t physicsSelection = ((fInputHandler->IsEventSelected())&AliVEvent::kMB);
 379   fEventCutInfoESD |= kPhysicsSelectionCut; // other alreay set via IsEventSelected
 380   fh1EventCutInfoESD->Fill(fEventCutInfoESD);
 381
 382
 383
 384   if(esdEventSelected) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kVertexIn;
 385   if(esdEventPileUp)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsPileUp;
 386   if(esdEventCosmic)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsCosmic;
 387   if(physicsSelection) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kPhysicsSelection;
 388
 389
 390   // here we have all selection information, fill histogram
 391   for(unsigned int i = 1;i<(UInt_t)fh1SelectionInfoESD->GetNbinsX();i++){
 392     if((i&fSelectionInfoESD)==i)fh1SelectionInfoESD->Fill(i);
 393   }
 394   AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD);
 395
 396   if(esd&&aod&&fDebug){
 397     if(esdEventSelected&&!aodEventSelected){
 398       Printf("%s:%d Different Selection for ESD and AOD",(char*)__FILE__,__LINE__);
 399       const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
 400       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
 401       Printf("ESD Vtx %s %s %d",vtxESD->GetName(),vtxESD->GetTitle(),vtxESD->GetNContributors());
 402       vtxESD->Print();
 403       Printf("AOD Vtx %s %s %d",vtxAOD->GetName(),vtxAOD->GetTitle(),vtxAOD->GetNContributors());
 404       vtxAOD->Print();
 405     }
 406   }
 407
 408   // loop over all possible triggers for esd
 409
 410   if(esd){
 411     const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
 412     //      Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print();
 413     TString vtxName(vtxESD->GetName());
 414     for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
 415       Bool_t esdTrig = kFALSE;
 416       esdTrig = AliAnalysisHelperJetTasks::IsTriggerFired(esd,(AliAnalysisHelperJetTasks::Trigger)it);
 417       if(esdTrig)fh2ESDTriggerCount->Fill(it,kAllTriggered);
 418       Bool_t cand = physicsSelection;
 419       if(cand){
 420         fh2ESDTriggerCount->Fill(it,kSelectedALICE);
 421       }
 422       if(!fUsePhysicsSelection)cand =  AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
 423       if(vtxESD->GetNContributors()>2&&!vtxName.Contains("TPCVertex")){
 424         if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredVertex);
 425         Float_t zvtx = vtxESD->GetZ();
 426         if(esdEventSelected&&esdTrig){
 427           fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn);
 428           fh2ESDTriggerVtx->Fill(it,zvtx);
 429         }
 430         if(cand)fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexValid);
 431       }
 432       if(cand&&esdEventSelected){
 433         fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexIn);
 434         fh2ESDTriggerCount->Fill(it,kSelected);
 435         AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
 436       }
 437     }
 438   }
 439
 440   if(aod){
 441     const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
 442     aodVtxValid = IsVertexValid(vtxAOD);
 443     aodVtxIn = IsVertexIn(vtxAOD);
 444
 445     for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
 446       Bool_t aodTrig = kFALSE;
 447       aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
 448       Bool_t cand = aod->GetHeader()->GetOfflineTrigger()&AliVEvent::kMB;
 449       if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
 450       if(aodVtxValid){
 451         if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredVertex);
 452         Float_t zvtx = vtxAOD->GetZ();
 453         if(cand&&aodEventSelected){
 454           fh2TriggerCount->Fill(it,kSelected);
 455         }
 456         if(aodTrig&&aodEventSelected){
 457           fh2TriggerVtx->Fill(it,zvtx);
 458           fh2TriggerCount->Fill(it,kTriggeredVertexIn);
 459         }
 460         if(cand){
 461           fh2TriggerCount->Fill(it,kSelectedALICEVertexValid);
 462           if(aodEventSelected)fh2TriggerCount->Fill(it,kSelectedALICEVertexIn);
 463         }
 464       }
 465     }
 466   }
 467
 468   if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
 469
 470
 471   Double_t ptHard = 0;
 472   Double_t nTrials = 1; // Trials for MC trigger
 473
 474   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 475   AliMCEvent* mcEvent = MCEvent();
 476   //    AliStack *pStack = 0;
 477   if(mcEvent){
 478     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 479     if(pythiaGenHeader){
 480       nTrials = pythiaGenHeader->Trials();
 481       ptHard  = pythiaGenHeader->GetPtHard();
 482       int iProcessType = pythiaGenHeader->ProcessType();
 483       // 11 f+f -> f+f
 484       // 12 f+barf -> f+barf
 485       // 13 f+barf -> g+g
 486       // 28 f+g -> f+g
 487       // 53 g+g -> f+barf
 488       // 68 g+g -> g+g
 489       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 490       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 491       fh1PtHard->Fill(ptHard);
 492       fh1PtHardTrials->Fill(ptHard,nTrials);
 493
 494     }// if pythia gen header
 495   }
 496
 497   // trigger selection
 498
 499
 500   PostData(1, fHistList);
 501 }
 502
 503 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
 504   if(!esd)return kFALSE;
 505   const AliESDVertex *vtx = esd->GetPrimaryVertex();
 506   return IsVertexIn(vtx); // vertex in calls vertex valid
 507 }
 508
 509
 510 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliAODEvent* aod) const {
 511   if(!aod)return kFALSE;
 512   const AliAODVertex *vtx = aod->GetPrimaryVertex();
 513   return IsVertexIn(vtx); // VertexIn calls VertexValid
 514 }
 515
 516 Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliESDVertex* vtx) {
 517
 518   // We can check the type of the vertex by its name and title
 519   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
 520   // vtx       name                  title
 521   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
 522   // Tracks    PrimaryVertex         VertexerTracksWithConstraint
 523   // TPC       TPCVertex             VertexerTracksNoConstraint
 524   // TPC       TPCVertex             VertexerTracksWithConstraint
 525   // SPD       SPDVertex             vertexer: 3D
 526   // SPD       SPDVertex             vertexer: Z
 527
 528   Int_t nCont = vtx->GetNContributors();
 529
 530   if(nCont>=1){
 531     fEventCutInfoESD |= kContributorsCut1;
 532     if(nCont>=2){
 533       fEventCutInfoESD |= kContributorsCut2;
 534       if(nCont>=3){
 535         fEventCutInfoESD |= kContributorsCut3;
 536       }
 537     }
 538   }
 539
 540   if(nCont<3)return kFALSE;
 541
 542   // do not want tpc only primary vertex
 543   TString vtxName(vtx->GetName());
 544   if(vtxName.Contains("TPCVertex")){
 545     fEventCutInfoESD |= kVertexTPC;
 546     return kFALSE;
 547   }
 548   if(vtxName.Contains("SPDVertex"))fEventCutInfoESD |= kVertexSPD;
 549   if(vtxName.Contains("PrimaryVertex"))fEventCutInfoESD |= kVertexGlobal;
 550
 551
 552   TString vtxTitle(vtx->GetTitle());
 553   if(vtxTitle.Contains("vertexer: Z")){
 554     if(vtx->GetDispersion()>0.02)return kFALSE;
 555   }
 556   fEventCutInfoESD |= kSPDDispersionCut;
 557   return kTRUE;
 558 }
 559
 560
 561 Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliAODVertex* vtx) const {
 562
 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   if(vtx->GetNContributors()<3)return kFALSE;
 572   // do not want tpc only primary vertex
 573   TString vtxName(vtx->GetName());
 574   if(vtxName.Contains("TPCVertex"))return kFALSE;
 575
 576   // no dispersion yet...
 577   /*
 578   TString vtxTitle(vtx->GetTitle());
 579   if(vtxTitle.Contains("vertexer: Z")){
 580     if(vtx->GetDispersion()>0.02)return kFALSE;
 581   }
 582   */
 583   return kTRUE;
 584 }
 585
 586
 587 Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx) {
 588
 589   if(!IsVertexValid(vtx))return kFALSE;
 590
 591   Float_t zvtx = vtx->GetZ();
 592   Float_t yvtx = vtx->GetY();
 593   Float_t xvtx = vtx->GetX();
 594
 595   xvtx -= fVtxXMean;
 596   yvtx -= fVtxYMean;
 597   zvtx -= fVtxZMean;
 598
 599
 600
 601   if(TMath::Abs(zvtx)>fVtxZCut){
 602     return kFALSE;
 603   }
 604   fEventCutInfoESD |= kVertexZCut;
 605   Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 606   if(r2>(fVtxRCut*fVtxRCut)){
 607     return kFALSE;
 608   }
 609   fEventCutInfoESD |= kVertexRCut;
 610   return kTRUE;
 611 }
 612
 613
 614 Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliAODVertex* vtx) const {
 615
 616   if(!IsVertexValid(vtx))return kFALSE;
 617
 618   Float_t zvtx = vtx->GetZ();
 619   Float_t yvtx = vtx->GetY();
 620   Float_t xvtx = vtx->GetX();
 621
 622   xvtx -= fVtxXMean;
 623   yvtx -= fVtxYMean;
 624   zvtx -= fVtxZMean;
 625
 626   Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 627
 628   Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
 629   return vertexIn;
 630 }
 631
 632 Bool_t AliAnalysisTaskJetServices::IsEventPileUp(const AliESDEvent* esd) const{
 633   if(!esd)return kFALSE;
 634   return esd->IsPileupFromSPD();
 635 }
 636
 637 Bool_t AliAnalysisTaskJetServices::IsEventCosmic(const AliESDEvent* esd) const {
 638   if(!esd)return kFALSE;
 639   // add track cuts for which we look for cosmics...
 640
 641   Bool_t isCosmic = kFALSE;
 642   Int_t nTracks = esd->GetNumberOfTracks();
 643   Int_t nCosmicCandidates = 0;
 644
 645   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
 646     AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
 647     if (!track1)  continue;
 648     UInt_t status1 = track1->GetStatus();
 649     //If track is ITS stand alone track, skip the track
 650     if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue;
 651     if(track1->Pt()<fPtMinCosmic) continue;
 652     //Start 2nd track loop to look for correlations
 653     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
 654       AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
 655       if(!track2) continue;
 656       UInt_t status2 = track2->GetStatus();
 657       //If track is ITS stand alone track, skip the track
 658       if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue;
 659       if(track2->Pt()<fPtMinCosmic) continue;
 660       //Check if back-to-back
 661       Double_t mom1[3],mom2[3];
 662       track1->GetPxPyPz(mom1);
 663       track2->GetPxPyPz(mom2);
 664       TVector3 momv1(mom1[0],mom1[1],mom1[2]);
 665       TVector3 momv2(mom2[0],mom2[1],mom2[2]);
 666       Float_t theta = (float)(momv1.Phi()-momv2.Phi());
 667       if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
 668
 669       Float_t deltaPhi = track1->Phi()-track2->Phi();
 670       if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
 671
 672       Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
 673       if(rIsol<fRIsolMinCosmic) continue;
 674
 675       if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
 676         nCosmicCandidates+=1;
 677         isCosmic = kTRUE;
 678       }
 679
 680     }
 681   }
 682
 683   fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
 684
 685   return isCosmic;
 686 }
 687
 688
 689
 690
 691 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
 692 {
 693   // Terminate analysis
 694   //
 695 }