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