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