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