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