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