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