PWGJE/AliAnalysisTaskJetServices.cxx

   1
   2 // **************************************
   3 // Task used for the correction of determiantion of reconstructed jet spectra
   4 // Compares input (gen) and output (rec) jets
   5 // *******************************************
   6
   7
   8 /**************************************************************************
   9  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  10  *                                                                        *
  11  * Author: The ALICE Off-line Project.                                    *
  12  * Contributors are mentioned in the code where appropriate.              *
  13  *                                                                        *
  14  * Permission to use, copy, modify and distribute this software and its   *
  15  * documentation strictly for non-commercial purposes is hereby granted   *
  16  * without fee, provided that the above copyright notice appears in all   *
  17  * copies and that both the copyright notice and this permission notice   *
  18  * appear in the supporting documentation. The authors make no claims     *
  19  * about the suitability of this software for any purpose. It is          *
  20  * provided "as is" without express or implied warranty.                  *
  21  **************************************************************************/
  22
  23
  24 #include <TROOT.h>
  25 #include <TRandom.h>
  26 #include <TString.h>
  27 #include <TSystem.h>
  28 #include <TInterpreter.h>
  29 #include <TChain.h>
  30 #include <TFile.h>
  31 #include <TKey.h>
  32 #include <TH1F.h>
  33 #include <TH2F.h>
  34 #include <TH3F.h>
  35 #include <TProfile.h>
  36 #include <TList.h>
  37 #include <TLorentzVector.h>
  38 #include <TRandom3.h>
  39 #include <TClonesArray.h>
  40 #include "TDatabasePDG.h"
  41
  42 #include "AliAnalysisTaskJetServices.h"
  43 #include "AliCentrality.h"
  44 #include "AliAnalysisDataContainer.h"
  45 #include "AliAnalysisDataSlot.h"
  46 #include "AliAnalysisManager.h"
  47 #include "AliJetFinder.h"
  48 #include "AliJetHeader.h"
  49 #include "AliJetReader.h"
  50 #include "AliJetReaderHeader.h"
  51 #include "AliUA1JetHeaderV1.h"
  52 #include "AliESDEvent.h"
  53 #include "AliAODEvent.h"
  54 #include "AliAODHandler.h"
  55 #include "AliAODTrack.h"
  56 #include "AliAODVZERO.h"
  57 #include "AliAODJet.h"
  58 #include "AliAODMCParticle.h"
  59 #include "AliMCEventHandler.h"
  60 #include "AliMCEvent.h"
  61 #include "AliStack.h"
  62 #include "AliGenPythiaEventHeader.h"
  63 #include "AliJetKineReaderHeader.h"
  64 #include "AliGenCocktailEventHeader.h"
  65 #include "AliInputEventHandler.h"
  66 #include "AliPhysicsSelection.h"
  67 #include "AliTriggerAnalysis.h"
  68
  69 #include "AliAnalysisHelperJetTasks.h"
  70
  71 ClassImp(AliAnalysisTaskJetServices)
  72
  73 AliAODVZERO*  AliAnalysisTaskJetServices::fgAODVZERO = NULL;
  74 AliAODHeader*  AliAnalysisTaskJetServices::fgAODHeader = NULL;
  75 TClonesArray*   AliAnalysisTaskJetServices::fgAODVertices = NULL;
  76
  77 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices():
  78   AliAnalysisTaskSE(),
  79   fUseAODInput(kFALSE),
  80   fUsePhysicsSelection(kFALSE),
  81   fMC(kFALSE),
  82   fFilterAODCollisions(kFALSE),
  83   fPhysicsSelectionFlag(AliVEvent::kMB),
  84   fSelectionInfoESD(0),
  85   fEventCutInfoESD(0),
  86   fFilterMask(0),
  87   fRPMethod(0),
  88   fCollisionType(kPbPb),
  89   fNTrigger(0),
  90   fAvgTrials(1),
  91   fVtxXMean(0),
  92   fVtxYMean(0),
  93   fVtxZMean(0),
  94   fVtxRCut(1.),
  95   fVtxZCut(8.),
  96   fVtxDeltaZCut(0.1),
  97   fPtMinCosmic(5.),
  98   fRIsolMinCosmic(3.),
  99   fMaxCosmicAngle(0.01),
 100   fCentrality(101),
 101   fTrackRecEtaWindow(0.9),
 102   fMinTrackPt(0.15),
 103   fRPAngle(0),
 104   fPsiVZEROA(0),
 105   fPsiVZEROC(0),
 106   fTriggerBit(0x0),
 107   fRandomizer(0),
 108   fNonStdFile(""),
 109   fTriggerName(0x0),
 110   fh1Xsec(0x0),
 111   fh1Trials(0x0),
 112   fh1AvgTrials(0x0),
 113   fh1PtHard(0x0),
 114   fh1PtHardTrials(0x0),
 115   fh1SelectionInfoESD(0x0),
 116   fh1EventCutInfoESD(0),
 117   fh1CentralityESD(0),
 118   fh1Centrality(0),
 119   fh1RP(0),
 120   fh2ReducedTrigger(0),
 121   fh2CentralityTriggerESD(0),
 122   fh2CentralityTrigger(0),
 123   fh2TriggerCount(0x0),
 124   fh2ESDTriggerCount(0x0),
 125   fh2TriggerVtx(0x0),
 126   fh2ESDTriggerVtx(0x0),
 127   fh2ESDTriggerRun(0x0),
 128   fh2VtxXY(0x0),
 129   fh1NCosmicsPerEvent(0x0),
 130   fp1RPXA(0x0),
 131   fp1RPYA(0x0),
 132   fp1RPXC(0x0),
 133   fp1RPYC(0x0),
 134   fp1CalibRPXA(0x0),
 135   fp1CalibRPYA(0x0),
 136   fp1CalibRPXC(0x0),
 137   fp1CalibRPYC(0x0),
 138   fh2RPAC(0x0),
 139   fh2RPAT(0x0),
 140   fh2RPCT(0x0),
 141   fh2XYA(0x0),
 142   fh2XYC(0x0),
 143   fh2RPCentrality(0x0),
 144   fh2RPACentrality(0x0),
 145   fh2RPCCentrality(0x0),
 146   fTriggerAnalysis(0x0),
 147   fHistList(0x0)
 148 {
 149   fRunRange[0] = fRunRange[1] = 0;
 150 }
 151
 152 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
 153   AliAnalysisTaskSE(name),
 154   fUseAODInput(kFALSE),
 155   fUsePhysicsSelection(kFALSE),
 156   fMC(kFALSE),
 157   fFilterAODCollisions(kFALSE),
 158   fPhysicsSelectionFlag(AliVEvent::kMB),
 159   fSelectionInfoESD(0),
 160   fEventCutInfoESD(0),
 161   fFilterMask(0),
 162   fRPMethod(0),
 163   fCollisionType(kPbPb),
 164   fNTrigger(0),
 165   fAvgTrials(1),
 166   fVtxXMean(0),
 167   fVtxYMean(0),
 168   fVtxZMean(0),
 169   fVtxRCut(1.),
 170   fVtxZCut(8.),
 171   fVtxDeltaZCut(0.1),
 172   fPtMinCosmic(5.),
 173   fRIsolMinCosmic(3.),
 174   fMaxCosmicAngle(0.01),
 175   fCentrality(101),
 176   fTrackRecEtaWindow(0.9),
 177   fMinTrackPt(0.15),
 178   fRPAngle(0),
 179   fPsiVZEROA(0),
 180   fPsiVZEROC(0),
 181   fTriggerBit(0x0),
 182   fRandomizer(0),
 183   fNonStdFile(""),
 184   fTriggerName(0x0),
 185   fh1Xsec(0x0),
 186   fh1Trials(0x0),
 187   fh1AvgTrials(0x0),
 188   fh1PtHard(0x0),
 189   fh1PtHardTrials(0x0),
 190   fh1SelectionInfoESD(0x0),
 191   fh1EventCutInfoESD(0),
 192   fh1CentralityESD(0),
 193   fh1Centrality(0),
 194   fh1RP(0),
 195   fh2ReducedTrigger(0),
 196   fh2CentralityTriggerESD(0),
 197   fh2CentralityTrigger(0),
 198   fh2TriggerCount(0x0),
 199   fh2ESDTriggerCount(0x0),
 200   fh2TriggerVtx(0x0),
 201   fh2ESDTriggerVtx(0x0),
 202   fh2ESDTriggerRun(0x0),
 203   fh2VtxXY(0x0),
 204   fh1NCosmicsPerEvent(0x0),
 205   fp1RPXA(0x0),
 206   fp1RPYA(0x0),
 207   fp1RPXC(0x0),
 208   fp1RPYC(0x0),
 209   fp1CalibRPXA(0x0),
 210   fp1CalibRPYA(0x0),
 211   fp1CalibRPXC(0x0),
 212   fp1CalibRPYC(0x0),
 213   fh2RPAC(0x0),
 214   fh2RPAT(0x0),
 215   fh2RPCT(0x0),
 216   fh2XYA(0x0),
 217   fh2XYC(0x0),
 218   fh2RPCentrality(0x0),
 219   fh2RPACentrality(0x0),
 220   fh2RPCCentrality(0x0),
 221   fTriggerAnalysis(0x0),
 222   fHistList(0x0)
 223 {
 224   fRunRange[0] = fRunRange[1] = 0;
 225   DefineOutput(1,TList::Class());
 226 }
 227
 228
 229
 230 Bool_t AliAnalysisTaskJetServices::Notify()
 231 {
 232   //
 233   // Implemented Notify() to read the cross sections
 234   // and number of trials from pyxsec.root
 235   //
 236
 237   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 238   Float_t xsection = 0;
 239   Float_t ftrials  = 1;
 240
 241   fAvgTrials = 1;
 242   if(tree){
 243     TFile *curfile = tree->GetCurrentFile();
 244     if (!curfile) {
 245       Error("Notify","No current file");
 246       return kFALSE;
 247     }
 248     if(!fh1Xsec||!fh1Trials){
 249       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 250       return kFALSE;
 251     }
 252     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 253     fh1Xsec->Fill("<#sigma>",xsection);
 254     // construct a poor man average trials
 255     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 256     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 257     fh1Trials->Fill("#sum{avg ntrials}",ftrials);
 258   }
 259   return kTRUE;
 260 }
 261
 262 void AliAnalysisTaskJetServices::UserCreateOutputObjects()
 263 {
 264
 265   //
 266   // Create the output container
 267   //
 268
 269   fRandomizer = new TRandom3(0);
 270   if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
 271
 272   OpenFile(1);
 273   if(!fHistList)fHistList = new TList();
 274   fHistList->SetOwner();
 275   PostData(1, fHistList); // post data in any case once
 276
 277   Bool_t oldStatus = TH1::AddDirectoryStatus();
 278   TH1::AddDirectory(kFALSE);
 279   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 280   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 281   fHistList->Add(fh1Xsec);
 282
 283   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 284   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 285   fHistList->Add(fh1Trials);
 286
 287   fh1AvgTrials = new TH1F("fh1AvgTrials","avg trials root file",1,0,1);
 288   fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
 289   fHistList->Add(fh1AvgTrials);
 290
 291   const Int_t nBinPt = 125;
 292   Double_t binLimitsPt[nBinPt+1];
 293   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 294     if(iPt == 0){
 295       binLimitsPt[iPt] = 0.0;
 296     }
 297     else {// 1.0
 298       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.;
 299     }
 300   }
 301
 302
 303   fh1CentralityESD = new TH1F("fh1CentralityESD","cent",103,-1,102);
 304   fHistList->Add(fh1CentralityESD);
 305
 306   fh1Centrality = new TH1F("fh1Centrality","cent",103,-1,102);
 307   fHistList->Add(fh1Centrality);
 308
 309   fh1RP = new TH1F("fh1RP","RP;#Psi",440, -1.*TMath::Pi(), 2.*TMath::Pi());
 310   fHistList->Add(fh1RP);
 311
 312   fh2ReducedTrigger = new TH2F("fh2ReducedTrigger","red trigger;red trigger;centrality",1<<fNTrigger,-0.5,(1<<fNTrigger)-0.5,100,-0.5,99.5);
 313   fHistList->Add(fh2ReducedTrigger);
 314
 315   fh2CentralityTriggerESD = new TH2F("fh2CentralityTriggerESD",";cent;trigger no",103,-1,102,fNTrigger,-0.5,fNTrigger-0.5);
 316   fHistList->Add(fh2CentralityTriggerESD);
 317
 318   fh2CentralityTrigger = new TH2F("fh2CentralityTrigger",";cent;trigger no",103,-1,102,fNTrigger,-0.5,fNTrigger-0.5);
 319   fHistList->Add(fh2CentralityTrigger);
 320
 321   for(int i = 0;i<fNTrigger;++i){
 322     fh2CentralityTriggerESD->GetYaxis()->SetBinLabel(i+1,fTriggerName[i].Data());
 323     fh2CentralityTrigger->GetYaxis()->SetBinLabel(i+1,fTriggerName[i].Data());
 324   }
 325
 326
 327   fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",6,-0.5,5.5,kConstraints,-0.5,kConstraints-0.5);
 328   fHistList->Add(fh2TriggerCount);
 329
 330   fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",6,-0.5,5.5,kConstraints,-0.5,kConstraints-0.5);
 331   fHistList->Add(fh2ESDTriggerCount);
 332   const Int_t nBins = 6*kConstraints;
 333   fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Constraint No. * (trig no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.);
 334   fHistList->Add(fh2TriggerVtx);
 335
 336   fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Constraint No.* (trg no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.);
 337   fHistList->Add(fh2ESDTriggerVtx);
 338
 339
 340   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 341   fHistList->Add(fh1PtHard);
 342   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 343   fHistList->Add(fh1PtHardTrials);
 344   fh1SelectionInfoESD = new TH1F("fh1SelectionInfoESD","Bit Masks that satisfy the selection info",
 345                                  AliAnalysisHelperJetTasks::kTotalSelections,0.5,AliAnalysisHelperJetTasks::kTotalSelections+0.5);
 346   fHistList->Add(fh1SelectionInfoESD);
 347
 348   fh1EventCutInfoESD = new TH1F("fh1EventCutInfoESD","Bit Masks that for the events after each step of cuts",
 349                                 kTotalEventCuts,0.5,kTotalEventCuts+0.5);
 350   fHistList->Add(fh1EventCutInfoESD);
 351
 352   // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range
 353   // 3 triggers BB BE/EB EE
 354
 355   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);
 356   fHistList->Add(fh2ESDTriggerRun);
 357
 358   fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
 359   fHistList->Add(fh2VtxXY);
 360   // =========== Switch on Sumw2 for all histos ===========
 361   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 362     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 363     if (h1){
 364       h1->Sumw2();
 365       continue;
 366     }
 367     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 368     if(hn)hn->Sumw2();
 369   }
 370
 371   fh1NCosmicsPerEvent = new TH1F("fh1NCosmicsPerEvent","Number of cosmic candidates per event",10,0.,10.);
 372   fHistList->Add(fh1NCosmicsPerEvent),
 373
 374
 375   fh2RPCentrality = new TH2F("fh2RPCentrality" ,"Reaction Plane Angle from tracks" , 20, 0.,100., 180, 0, TMath::Pi());
 376   fHistList->Add(fh2RPCentrality);
 377
 378
 379   fh2RPACentrality = new TH2F("fh2RPACentrality" ,"Reaction Plane Angle from vzero A" , 20, 0.,100., 180, 0, TMath::Pi());
 380   fHistList->Add(fh2RPACentrality);
 381
 382   fh2RPCCentrality = new TH2F("fh2RPCCentrality" ,"Reaction Plane Angle from vzero C" , 20, 0.,100., 180, 0, TMath::Pi());
 383   fHistList->Add(fh2RPCCentrality);
 384
 385   fh2RPAC = new TH2F("fh2RPAC" ,"Reaction Plane Angle vzero a vs c" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
 386   fHistList->Add( fh2RPAC);
 387
 388   fh2RPAT = new TH2F("fh2RPAT" ,"Reaction Plane Angle vzero a vs tracks" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
 389   fHistList->Add( fh2RPAT);
 390
 391   fh2RPCT = new TH2F("fh2RPCT" ,"Reaction Plane Angle vzero c vs tracks" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
 392   fHistList->Add( fh2RPCT);
 393
 394   fh2XYA = new TH2F("fh2XYA" ,"XY vzeroa ;X;Y;" ,100,-0.3,0.3,100,-0.3,0.3);
 395   fHistList->Add(fh2XYA);
 396   fh2XYC = new TH2F("fh2XYC" ,"XY vzeroc ;X;Y;" ,100,-0.3,0.3,100,-0.3,0.3);
 397   fHistList->Add(fh2XYC);
 398
 399   // profiles for mean
 400   fp1RPXA = new TProfile("fp1RPXA","mean vzeroa x vs run number;run;x",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5);
 401   fHistList->Add(fp1RPXA);
 402
 403   fp1RPYA = new TProfile("fp1RPYA","mean vzeroa y vs run number;run;y",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5);
 404   fHistList->Add(fp1RPYA);
 405
 406
 407   fp1RPXC = new TProfile("fp1RPXC","mean vzeroc x vs run number;run;x",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5);
 408   fHistList->Add(fp1RPXC);
 409
 410   fp1RPYC = new TProfile("fp1RPYC","mean vzeroa y vs run number;run;y",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5);
 411   fHistList->Add(fp1RPYC);
 412
 413
 414   TH1::AddDirectory(oldStatus);
 415
 416   // Add an AOD branch for replication
 417   if(fNonStdFile.Length()){
 418      if (fDebug > 1) AliInfo("Replicating header");
 419      fgAODHeader = new AliAODHeader;
 420      AddAODBranch("AliAODHeader",&fgAODHeader,fNonStdFile.Data());
 421      if (fDebug > 1) AliInfo("Replicating vzeros");
 422      fgAODVZERO = new AliAODVZERO;
 423      AddAODBranch("AliAODVZERO",&fgAODVZERO,fNonStdFile.Data());
 424      if (fDebug > 1) AliInfo("Replicating primary vertices");
 425      fgAODVertices = new TClonesArray("AliAODVertex",3);
 426      fgAODVertices->SetName("vertices");
 427      AddAODBranch("TClonesArray",&fgAODVertices,fNonStdFile.Data());
 428   }
 429 }
 430
 431 void AliAnalysisTaskJetServices::Init()
 432 {
 433   //
 434   // Initialization
 435   //
 436   if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
 437
 438 }
 439
 440 void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
 441 {
 442
 443   //
 444   // Execute analysis for current event
 445   //
 446
 447   AliAODEvent *aod = 0;
 448   AliESDEvent *esd = 0;
 449
 450
 451
 452   AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
 453   AliAnalysisHelperJetTasks::EventClass(kTRUE,0);
 454   AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,0); // set slection to false
 455   fSelectionInfoESD = 0; // reset
 456   fEventCutInfoESD = 0; // reset
 457   AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
 458
 459
 460   static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 461
 462
 463
 464   if(fUseAODInput){
 465     aod = dynamic_cast<AliAODEvent*>(InputEvent());
 466     if(!aod){
 467       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
 468       return;
 469     }
 470     // fethc the header
 471   }
 472   else{
 473     //  assume that the AOD is in the general output...
 474     aod  = AODEvent();
 475     if(!aod){
 476       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 477       return;
 478     }
 479     esd = dynamic_cast<AliESDEvent*>(InputEvent());
 480   }
 481   if(aod&&fDebug>2){
 482     aod->Print();
 483     Printf("Vertices %d",aod->GetNumberOfVertices());
 484     Printf("tracks %d",aod->GetNumberOfTracks());
 485     Printf("jets %d",aod->GetNJets());
 486   }
 487   fSelectionInfoESD |= kNoEventCut;
 488   fEventCutInfoESD |= kNoEventCut;
 489
 490   Bool_t esdVtxValid = false;
 491   Bool_t esdVtxIn = false;
 492   Bool_t aodVtxValid = false;
 493   Bool_t aodVtxIn = false;
 494
 495   // Apply additional constraints
 496   Bool_t esdEventSelected = IsEventSelected(esd);
 497   Bool_t esdEventPileUp = IsEventPileUp(esd);
 498   Bool_t esdEventCosmic = IsEventCosmic(esd);
 499
 500   Bool_t aodEventSelected = IsEventSelected(aod);
 501
 502   Bool_t physicsSelection = ((fInputHandler->IsEventSelected())&fPhysicsSelectionFlag);
 503
 504
 505   fEventCutInfoESD |= kPhysicsSelectionCut; // other alreay set via IsEventSelected
 506   fh1EventCutInfoESD->Fill (fEventCutInfoESD);
 507
 508   if(esdEventSelected) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kVertexIn;
 509   if(esdEventPileUp)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsPileUp;
 510   if(esdEventCosmic)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsCosmic;
 511   if(physicsSelection) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kPhysicsSelection;
 512
 513
 514
 515
 516   // here we have all selection information, fill histogram
 517   for(unsigned int i = 1;i<(UInt_t)fh1SelectionInfoESD->GetNbinsX();i++){
 518     if((i&fSelectionInfoESD)==i)fh1SelectionInfoESD->Fill(i);
 519   }
 520   AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD);
 521
 522   if(esd&&aod&&fDebug){
 523     if(esdEventSelected&&!aodEventSelected){
 524       Printf("%s:%d Different Selection for ESD and AOD",(char*)__FILE__,__LINE__);
 525       const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
 526       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
 527       Printf("ESD Vtx %s %s %d",vtxESD->GetName(),vtxESD->GetTitle(),vtxESD->GetNContributors());
 528       vtxESD->Print();
 529       Printf("AOD Vtx %s %s %d",vtxAOD->GetName(),vtxAOD->GetTitle(),vtxAOD->GetNContributors());
 530       vtxAOD->Print();
 531     }
 532   }
 533
 534   // loop over all possible triggers for esd
 535
 536   Float_t cent = 0;
 537   if(fCollisionType==kPbPb){
 538     if(aod)cent = aod->GetHeader()->GetCentrality();
 539     if(fDebug)Printf("%s:%d %3.3f",(char*)__FILE__,__LINE__,cent);
 540     if(cent<0)cent = 101;
 541   }
 542   fCentrality = cent;
 543   fRPAngle = 0;
 544   fPsiVZEROA = 0;
 545   fPsiVZEROC = 0;
 546
 547
 548   if(esd){
 549     const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
 550     const AliESDVertex *vtxESDSPD = esd->GetPrimaryVertexSPD();
 551     esdVtxValid = IsVertexValid(vtxESD);
 552     esdVtxIn = IsVertexIn(vtxESD,vtxESDSPD);
 553     if(aodH&&physicsSelection&&fFilterAODCollisions&&aod){
 554       if(fDebug)Printf("%s:%d Centrality %3.3f vtxin %d",(char*)__FILE__,__LINE__,cent,esdVtxIn);
 555       if(cent<=80&&esdVtxIn){
 556         aodH->SetFillAOD(kTRUE);
 557         aodH->SetFillExtension(kTRUE);
 558       }
 559     }
 560
 561
 562     Float_t zvtx = vtxESD->GetZ();
 563     Int_t  iCl = GetEventClass(esd);
 564     AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
 565     Bool_t cand = physicsSelection;
 566
 567     if(fDebug)Printf("%s:%d %d %d %d Icl %d",(char*)__FILE__,__LINE__,esdVtxValid,esdVtxIn,cand,iCl);
 568     fh2ESDTriggerCount->Fill(0.,kAllTriggered);
 569     fh2ESDTriggerCount->Fill(iCl,kAllTriggered);
 570     if(cand){
 571       fh2ESDTriggerCount->Fill(0.,kSelectedALICE);
 572       fh2ESDTriggerCount->Fill(iCl,kSelectedALICE);
 573       fh2ESDTriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
 574     }
 575     //    if(!fUsePhysicsSelection)cand =  AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
 576     if(esdVtxValid){
 577       fh2ESDTriggerCount->Fill(0.,kTriggeredVertex);
 578       fh2ESDTriggerCount->Fill(iCl,kTriggeredVertex);
 579       fh2ESDTriggerVtx->Fill(iCl,zvtx);
 580       if(esdVtxIn){
 581         fh2ESDTriggerCount->Fill(0.,kTriggeredVertexIn);
 582         fh2ESDTriggerCount->Fill(iCl,kTriggeredVertexIn);
 583         fh2ESDTriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
 584       }
 585       if(cand){
 586         fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexValid);
 587         fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexValid);
 588         fh2ESDTriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
 589       }
 590     }
 591
 592     if(cand&&esdVtxIn&&iCl<5){
 593       fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexIn);
 594       fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexIn);
 595       fh2ESDTriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
 596       fh2ESDTriggerVtx->Fill(kSelected*(iCl+1),zvtx);
 597       fh2ESDTriggerCount->Fill(iCl,kSelected);
 598       fh2ESDTriggerCount->Fill(0.,kSelected);
 599       AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
 600       if(esd->GetCentrality()){
 601         Float_t tmpCent = 100;
 602         tmpCent = esd->GetCentrality()->GetCentralityPercentile("V0M");
 603         if(tmpCent<0)tmpCent = 101;
 604         fh1CentralityESD->Fill(tmpCent);
 605         UInt_t ir = 0;
 606         for(int it = 0;it<fNTrigger;++it){
 607           if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
 608             fh2CentralityTriggerESD->Fill(tmpCent,it);
 609             ir |= (1<<it);
 610           }
 611         }
 612         fh2ReducedTrigger->Fill(ir,cent);
 613       }
 614     }
 615   }
 616
 617
 618
 619   if(aod){
 620     const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
 621     const AliAODVertex *vtxAODSPD = aod->GetPrimaryVertexSPD();
 622     aodVtxValid = IsVertexValid(vtxAOD);
 623     aodVtxIn = IsVertexIn(vtxAOD,vtxAODSPD);
 624     Float_t zvtx = vtxAOD->GetZ();
 625     Int_t  iCl = GetEventClass(aod);
 626     AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
 627     Bool_t cand = aod->GetHeader()->GetOfflineTrigger()&fPhysicsSelectionFlag;
 628     if(fDebug)Printf("%s:%d AOD selection %d %d",(char*)__FILE__,__LINE__,cand,aod->GetHeader()->GetOfflineTrigger());
 629     fh2TriggerCount->Fill(0.,kAllTriggered);
 630     fh2TriggerCount->Fill(iCl,kAllTriggered);
 631     if(cand){
 632       fh2TriggerCount->Fill(0.,kSelectedALICE);
 633       fh2TriggerCount->Fill(iCl,kSelectedALICE);
 634       fh2TriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
 635     }
 636     if(aodVtxValid){
 637       fh2TriggerCount->Fill(0.,kTriggeredVertex);
 638       fh2TriggerCount->Fill(iCl,kTriggeredVertex);
 639       fh2TriggerVtx->Fill(iCl,zvtx);
 640       if(aodVtxIn){
 641         fh2TriggerCount->Fill(0.,kTriggeredVertexIn);
 642         fh2TriggerCount->Fill(iCl,kTriggeredVertexIn);
 643         fh2TriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
 644       }
 645       if(cand){
 646         fh2TriggerCount->Fill(0.,kSelectedALICEVertexValid);
 647         fh2TriggerCount->Fill(iCl,kSelectedALICEVertexValid);
 648         fh2TriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
 649       }
 650     }
 651
 652     if(cand&&aodVtxIn&&iCl<5){
 653       fh2TriggerCount->Fill(0.,kSelectedALICEVertexIn);
 654       fh2TriggerCount->Fill(iCl,kSelectedALICEVertexIn);
 655       fh2TriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
 656       fh2TriggerVtx->Fill(kSelected*(iCl+1),zvtx);
 657       fh2TriggerCount->Fill(iCl,kSelected);
 658       fh2TriggerCount->Fill(0.,kSelected);
 659       fh1Centrality->Fill(cent);
 660       for(int it = 0;it<fNTrigger;++it){
 661         if(fInputHandler->IsEventSelected()&fTriggerBit[it])fh2CentralityTrigger->Fill(cent,it);
 662       }
 663       AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
 664       if(aodH&&cand&&fFilterAODCollisions&&!esd){
 665         if(fCentrality<=80&&aodVtxIn){
 666           aodH->SetFillAOD(kTRUE);
 667           aodH->SetFillExtension(kTRUE);
 668         }
 669       }
 670
 671       TList recTracks;
 672       GetListOfTracks(&recTracks);
 673       CalculateReactionPlaneAngleVZERO(aod);
 674       fRPAngle = aod->GetHeader()->GetEventplane();
 675       fh1RP->Fill(fRPAngle);
 676       fh2RPCentrality->Fill(fCentrality,fRPAngle);
 677       fh2RPACentrality->Fill(fCentrality,fPsiVZEROA);
 678       fh2RPCCentrality->Fill(fCentrality,fPsiVZEROC);
 679       fh2RPAC->Fill(fPsiVZEROA,fPsiVZEROC);
 680       fh2RPAT->Fill(fPsiVZEROA,fRPAngle);
 681       fh2RPCT->Fill(fPsiVZEROC,fRPAngle);
 682       if(fRPMethod==kRPTracks)AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fRPAngle); // set slection to false
 683       else if(fRPMethod==kRPVZEROA)AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fPsiVZEROA); // set slection to false
 684       else if(fRPMethod==kRPVZEROC)AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fPsiVZEROA); // set slection to false
 685
 686       if(fUseAODInput&&fCentrality<=80){
 687         if(fFilterAODCollisions&&aod){
 688           aodH->SetFillAOD(kTRUE);
 689           aodH->SetFillExtension(kTRUE);
 690         }
 691       }
 692     }
 693   }
 694
 695   if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
 696
 697
 698   Double_t ptHard = 0;
 699   Double_t nTrials = 1; // Trials for MC trigger
 700
 701   fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
 702   AliMCEvent* mcEvent = MCEvent();
 703   //    AliStack *pStack = 0;
 704   if(mcEvent){
 705     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 706     if(pythiaGenHeader){
 707       nTrials = pythiaGenHeader->Trials();
 708       ptHard  = pythiaGenHeader->GetPtHard();
 709       int iProcessType = pythiaGenHeader->ProcessType();
 710       // 11 f+f -> f+f
 711       // 12 f+barf -> f+barf
 712       // 13 f+barf -> g+g
 713       // 28 f+g -> f+g
 714       // 53 g+g -> f+barf
 715       // 68 g+g -> g+g
 716       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 717       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 718       fh1PtHard->Fill(ptHard);
 719       fh1PtHardTrials->Fill(ptHard,nTrials);
 720
 721     }// if pythia gen header
 722   }
 723
 724   // trigger selection
 725
 726   // replication of
 727   if(fNonStdFile.Length()&&aod){
 728     if (fgAODHeader){
 729       *fgAODHeader =  *(dynamic_cast<AliAODHeader*>(aod->GetHeader()));
 730       Double_t q[kRPMethods] = {fRPAngle,fPsiVZEROA,fPsiVZEROC};
 731       fgAODHeader->SetQTheta(q,kRPMethods);
 732     }
 733     if (fgAODVZERO){
 734       *fgAODVZERO =  *(dynamic_cast<AliAODVZERO*>(aod->GetVZEROData()));
 735     }
 736     if(fgAODVertices){
 737       fgAODVertices->Delete();
 738       TClonesArray &vertices = *fgAODVertices;
 739       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
 740       // we only use some basic information,
 741
 742
 743       Double_t pos[3];
 744       Double_t covVtx[6];
 745
 746       vtxAOD->GetXYZ(pos); // position
 747       vtxAOD->GetCovMatrix(covVtx); //covariance matrix
 748       Int_t jVertices = 0;
 749       AliAODVertex * vtx = new(vertices[jVertices++])
 750         AliAODVertex(pos, covVtx, vtxAOD->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
 751       vtx->SetName(vtxAOD->GetName());
 752       vtx->SetTitle(vtxAOD->GetTitle());
 753       TString vtitle = vtxAOD->GetTitle();
 754       vtx->SetNContributors(vtxAOD->GetNContributors());
 755
 756       // Add SPD "main" vertex
 757       const AliAODVertex *vtxS = aod->GetPrimaryVertexSPD();
 758       vtxS->GetXYZ(pos); // position
 759       vtxS->GetCovMatrix(covVtx); //covariance matrix
 760       AliAODVertex * mVSPD = new(vertices[jVertices++])
 761         AliAODVertex(pos, covVtx, vtxS->GetChi2perNDF(), NULL, -1, AliAODVertex::kMainSPD);
 762       mVSPD->SetName(vtxS->GetName());
 763       mVSPD->SetTitle(vtxS->GetTitle());
 764       mVSPD->SetNContributors(vtxS->GetNContributors());
 765
 766       // Add tpc only vertex
 767       if(esd){
 768         const AliESDVertex *vtxT =  esd->GetPrimaryVertexTPC();
 769         vtxT->GetXYZ(pos); // position
 770         vtxT->GetCovMatrix(covVtx); //covariance matrix
 771         AliAODVertex * mVTPC = new(vertices[jVertices++])
 772           AliAODVertex(pos, covVtx, vtxT->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainTPC);
 773         mVTPC->SetName(vtxT->GetName());
 774         mVTPC->SetTitle(vtxT->GetTitle());
 775         mVTPC->SetNContributors(vtxT->GetNContributors());
 776       }
 777     }
 778   }
 779
 780   PostData(1, fHistList);
 781 }
 782
 783 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
 784   if(!esd)return kFALSE;
 785   const AliESDVertex *vtx = esd->GetPrimaryVertex();
 786   const AliESDVertex *vtxSPD = esd->GetPrimaryVertex();
 787   return IsVertexIn(vtx,vtxSPD); // vertex in calls vertex valid
 788 }
 789
 790 AliAnalysisTaskJetServices::~AliAnalysisTaskJetServices(){
 791   if(fgAODVertices){
 792     fgAODVertices->Delete();
 793     delete fgAODVertices;
 794   }
 795   delete fRandomizer;
 796   if(fgAODHeader)delete fgAODHeader;
 797   if(fgAODVZERO)delete fgAODVZERO;
 798   delete fp1CalibRPXA;
 799   delete fp1CalibRPYA;
 800   delete fp1CalibRPXC;
 801   delete fp1CalibRPYC;
 802   delete [] fTriggerBit;
 803   delete [] fTriggerName;
 804
 805 }
 806
 807
 808 void AliAnalysisTaskJetServices::SetV0Centroids(TProfile *xa,TProfile *ya,TProfile *xc, TProfile *yc){
 809
 810   if(xa){
 811     if(fp1CalibRPXA)delete fp1CalibRPXA;
 812     fp1CalibRPXA =  (TProfile*)xa->Clone(Form("%sCalib",xa->GetName()));
 813   }
 814   else{
 815     Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
 816   }
 817   if(ya){
 818     if(fp1CalibRPYA)delete fp1CalibRPYA;
 819     fp1CalibRPYA =  (TProfile*)ya->Clone(Form("%sCalib",ya->GetName()));
 820   }
 821   else{
 822     Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
 823   }
 824   if(xc){
 825     if(fp1CalibRPXC)delete fp1CalibRPXC;
 826     fp1CalibRPXC =  (TProfile*)xc->Clone(Form("%sCalib",xc->GetName()));
 827   }
 828   else{
 829     Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
 830   }
 831   if(ya){
 832     if(fp1CalibRPYC)delete fp1CalibRPYC;
 833     fp1CalibRPYC =  (TProfile*)yc->Clone(Form("%sCalib",yc->GetName()));
 834   }
 835   else{
 836     Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
 837   }
 838 }
 839
 840 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliAODEvent* aod) const {
 841   if(!aod)return kFALSE;
 842   const AliAODVertex *vtx = aod->GetPrimaryVertex();
 843   const AliAODVertex *vtxSPD = aod->GetPrimaryVertexSPD();
 844   return IsVertexIn(vtx,vtxSPD); // VertexIn calls VertexValid
 845 }
 846
 847 Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliESDVertex* vtx) {
 848
 849   // We can check the type of the vertex by its name and title
 850   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
 851   // vtx       name                  title
 852   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
 853   // Tracks    PrimaryVertex         VertexerTracksWithConstraint
 854   // TPC       TPCVertex             VertexerTracksNoConstraint
 855   // TPC       TPCVertex             VertexerTracksWithConstraint
 856   // SPD       SPDVertex             vertexer: 3D
 857   // SPD       SPDVertex             vertexer: Z
 858
 859   if(!vtx){
 860     Printf("%s:%d No ESD vertex found",(char*)__FILE__,__LINE__);
 861     return kFALSE;
 862   }
 863   Int_t nCont = vtx->GetNContributors();
 864   if(nCont>=1){
 865     fEventCutInfoESD |= kContributorsCut1;
 866     if(nCont>=2){
 867       fEventCutInfoESD |= kContributorsCut2;
 868       if(nCont>=3){
 869         fEventCutInfoESD |= kContributorsCut3;
 870       }
 871     }
 872   }
 873
 874   if(nCont<3)return kFALSE;
 875
 876   // do not want tpc only primary vertex
 877   TString vtxName(vtx->GetName());
 878   if(vtxName.Contains("TPCVertex")){
 879     fEventCutInfoESD |= kVertexTPC;
 880     return kFALSE;
 881   }
 882   if(vtxName.Contains("SPDVertex"))fEventCutInfoESD |= kVertexSPD;
 883   if(vtxName.Contains("PrimaryVertex"))fEventCutInfoESD |= kVertexGlobal;
 884
 885
 886   TString vtxTitle(vtx->GetTitle());
 887   if(vtxTitle.Contains("vertexer: Z")){
 888     if(vtx->GetDispersion()>0.02)return kFALSE;
 889   }
 890   fEventCutInfoESD |= kSPDDispersionCut;
 891   return kTRUE;
 892 }
 893
 894
 895 Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliAODVertex* vtx) const {
 896
 897   // We can check the type of the vertex by its name and title
 898   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
 899   // vtx       name                  title
 900   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
 901   // TPC       TPCVertex             VertexerTracksNoConstraint
 902   // SPD       SPDVertex             vertexer: 3D
 903   // SPD       SPDVertex             vertexer: Z
 904
 905   if(!vtx){
 906     Printf("%s:%d No AOD vertex found",(char*)__FILE__,__LINE__);
 907     return kFALSE;
 908   }
 909
 910   TString vtxName(vtx->GetName());
 911   if(fDebug){
 912     Printf(" n contrib %d",vtx->GetNContributors());
 913     Printf("vtxName: %s",vtxName.Data());
 914     vtx->Print();
 915   }
 916
 917   //  if(vtx->GetNContributors()<3)return kFALSE;
 918   // do not want tpc only primary vertex
 919
 920   if(vtxName.Contains("TPCVertex"))return kFALSE;
 921
 922   // no dispersion yet...
 923   /*
 924   TString vtxTitle(vtx->GetTitle());
 925   if(vtxTitle.Contains("vertexer: Z")){
 926     if(vtx->GetDispersion()>0.02)return kFALSE;
 927   }
 928   */
 929   return kTRUE;
 930 }
 931
 932
 933 Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx,const AliESDVertex* vtxSPD) {
 934
 935   if(!IsVertexValid(vtx))return kFALSE;
 936
 937   Float_t zvtx = vtx->GetZ();
 938   Float_t yvtx = vtx->GetY();
 939   Float_t xvtx = vtx->GetX();
 940   Float_t deltaZ = zvtx - vtxSPD->GetZ();
 941   if(fDebug)Printf("%s:%d deltaz %3.3f ",__FILE__,__LINE__,deltaZ);
 942
 943   xvtx -= fVtxXMean;
 944   yvtx -= fVtxYMean;
 945   zvtx -= fVtxZMean;
 946
 947
 948
 949   if(TMath::Abs(zvtx)>fVtxZCut){
 950     return kFALSE;
 951   }
 952   fEventCutInfoESD |= kVertexZCut;
 953   Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 954   if(r2>(fVtxRCut*fVtxRCut)){
 955     return kFALSE;
 956   }
 957   fEventCutInfoESD |= kVertexRCut;
 958
 959
 960
 961
 962   if(TMath::Abs(deltaZ)>fVtxDeltaZCut)return kFALSE;
 963   fEventCutInfoESD |= kVertexDeltaZCut;
 964   return kTRUE;
 965 }
 966
 967
 968 Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliAODVertex* vtx,const AliAODVertex* vtxSPD) const {
 969
 970   if(!IsVertexValid(vtx))return kFALSE;
 971
 972   Float_t zvtx = vtx->GetZ();
 973   Float_t yvtx = vtx->GetY();
 974   Float_t xvtx = vtx->GetX();
 975   Float_t deltaZ = zvtx - vtxSPD->GetZ();
 976   if(fDebug)Printf("%s:%d deltaz %3.3f ",__FILE__,__LINE__,deltaZ);
 977
 978   xvtx -= fVtxXMean;
 979   yvtx -= fVtxYMean;
 980   zvtx -= fVtxZMean;
 981
 982   Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 983
 984   Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
 985
 986
 987   if(TMath::Abs(deltaZ)>fVtxDeltaZCut)vertexIn = kFALSE;
 988   if(fDebug)Printf("%s:%d vertexIn %d ",__FILE__,__LINE__,vertexIn);
 989
 990
 991   return vertexIn;
 992 }
 993
 994 Bool_t AliAnalysisTaskJetServices::IsEventPileUp(const AliESDEvent* esd) const{
 995   if(!esd)return kFALSE;
 996   return esd->IsPileupFromSPD();
 997 }
 998
 999 Bool_t AliAnalysisTaskJetServices::IsEventCosmic(const AliESDEvent* esd) const {
1000   if(!esd)return kFALSE;
1001   // add track cuts for which we look for cosmics...
1002
1003   Bool_t isCosmic = kFALSE;
1004   Int_t nTracks = esd->GetNumberOfTracks();
1005   Int_t nCosmicCandidates = 0;
1006
1007   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
1008     AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
1009     if (!track1)  continue;
1010     UInt_t status1 = track1->GetStatus();
1011     //If track is ITS stand alone track, skip the track
1012     if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue;
1013     if(track1->Pt()<fPtMinCosmic) continue;
1014     //Start 2nd track loop to look for correlations
1015     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
1016       AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
1017       if(!track2) continue;
1018       UInt_t status2 = track2->GetStatus();
1019       //If track is ITS stand alone track, skip the track
1020       if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue;
1021       if(track2->Pt()<fPtMinCosmic) continue;
1022       //Check if back-to-back
1023       Double_t mom1[3],mom2[3];
1024       track1->GetPxPyPz(mom1);
1025       track2->GetPxPyPz(mom2);
1026       TVector3 momv1(mom1[0],mom1[1],mom1[2]);
1027       TVector3 momv2(mom2[0],mom2[1],mom2[2]);
1028       Float_t theta = (float)(momv1.Phi()-momv2.Phi());
1029       if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
1030
1031       Float_t deltaPhi = track1->Phi()-track2->Phi();
1032       if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
1033
1034       Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
1035       if(rIsol<fRIsolMinCosmic) continue;
1036
1037       if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
1038         nCosmicCandidates+=1;
1039         isCosmic = kTRUE;
1040       }
1041
1042     }
1043   }
1044
1045   fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
1046
1047   return isCosmic;
1048 }
1049
1050
1051 Int_t AliAnalysisTaskJetServices::GetEventClass(AliESDEvent *esd){
1052
1053   Float_t cent = 0;
1054   if(fCollisionType==kPbPb){
1055     if(esd->GetCentrality()){
1056       cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
1057     }
1058     if(cent<0)cent = 101;
1059     if(cent>80||cent<0)return 5;
1060     if(cent>50)return 4;
1061     if(cent>30)return 3;
1062     if(cent>10)return 2;
1063     return 1;
1064   }
1065   return 1;
1066 }
1067
1068
1069 Int_t AliAnalysisTaskJetServices::GetEventClass(AliAODEvent *aod){
1070
1071   if(fCollisionType==kPbPb){
1072     Float_t cent = aod->GetHeader()->GetCentrality();
1073     if(cent>80||cent<0)return 5;
1074     if(cent>50)return 4;
1075     if(cent>30)return 3;
1076     if(cent>10)return 2;
1077     return 1;
1078   }
1079   return 1;
1080
1081 }
1082
1083
1084 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
1085 {
1086   // Terminate analysis
1087   //
1088 }
1089
1090 Bool_t AliAnalysisTaskJetServices::CalculateReactionPlaneAngleVZERO(AliAODEvent *aod){
1091
1092   //  const Double_t arr_eta[11]={-3.7, -3.2, -2.7, -2.2, -1.7,0, 2.8, 3.4, 3.9, 4.5,5.1};
1093   if(!aod)return kFALSE;
1094   static bool bFirst = true;
1095   static Double_t v0phi[64] = {0,};
1096
1097   if(bFirst){
1098     int is=0;
1099     for(int iArm = 0; iArm<2; iArm++){
1100       for(int iRing = 0; iRing<4 ; iRing++){
1101         for(int iSec = 0; iSec<8 ; iSec++){
1102           v0phi[is] = 22.5 + 45. * iSec;
1103           v0phi[is] *= TMath::Pi()/180;
1104           // cout<< is <<" "<< v0phi[is]<<endl;
1105           is++;
1106         }
1107       }
1108     }
1109     bFirst = false;
1110   }
1111
1112   //
1113   const AliAODVZERO *aodVZERO = aod->GetVZEROData();
1114   Double_t numYZNA = 0,numXZNA = 0,sumZNA = 0;
1115   Double_t meanXA = 0,meanYA = 0;
1116
1117   Double_t numYZNC = 0,numXZNC = 0,sumZNC = 0;
1118   Double_t meanXC = 0,meanYC = 0;
1119
1120
1121
1122   static Int_t iOldRun = -1;
1123   static Int_t iFoundBin = -1;
1124
1125   if(aod->GetRunNumber()!=iOldRun&&(fp1CalibRPYA)){
1126     // search only or the bin in case of new runs
1127     iFoundBin = -1;
1128     Int_t ib = fp1CalibRPYA->FindBin(aod->GetRunNumber());
1129     Float_t err = fp1CalibRPYA->GetBinError(ib);
1130     if(err>0){// value can be zero...
1131       iFoundBin = ib;
1132     }
1133     else{
1134       Int_t ibLo = ib-1;
1135       Int_t ibUp = ib+1;
1136       while(iFoundBin<0&&(ibLo>0||ibUp<=fp1CalibRPYA->GetNbinsX())){
1137         err = fp1CalibRPYA->GetBinError(ibLo);
1138         if(err>0){
1139           iFoundBin = ibLo;
1140         }
1141         else{
1142           err = fp1CalibRPYA->GetBinError(ibUp);
1143           if(err>0)iFoundBin = ibUp;
1144         }
1145         ibUp++;
1146         ibLo--;
1147       }
1148     }
1149     iOldRun = aod->GetRunNumber();
1150   }
1151
1152   if(fDebug)Printf("%s:%d iFoundBin %d",(char*)__FILE__,__LINE__,iFoundBin);
1153
1154   if(iFoundBin>0&&(fp1CalibRPYA)){
1155     meanXA = fp1CalibRPXA->GetBinContent(iFoundBin);
1156     meanYA = fp1CalibRPYA->GetBinContent(iFoundBin);
1157     meanXC = fp1CalibRPXC->GetBinContent(iFoundBin);
1158     meanYC = fp1CalibRPYC->GetBinContent(iFoundBin);
1159   }
1160
1161   if(fDebug)Printf("%s:%d iFoundBin %1.3E %1.3E %1.3E %1.3E",(char*)__FILE__,__LINE__,meanXA,meanYA,meanXC,meanYC);
1162
1163   for (int i=0; i<64; i++) {
1164     Double_t mult = aodVZERO->GetMultiplicity(i);
1165     Double_t phi = v0phi[i];
1166     if (mult>0) {
1167       if (i<32) { //C-side
1168         Double_t wZNC= mult;
1169         numYZNC += sin(2.*phi)*wZNC;
1170         numXZNC += cos(2.*phi)*wZNC;
1171         sumZNC+=wZNC;
1172       }
1173       else if(i>31){ //A-side
1174         Double_t wZNA=mult;
1175         numYZNA += sin(2.*phi)*wZNA;
1176         numXZNA += cos(2.*phi)*wZNA;
1177         sumZNA+=wZNA;
1178       }
1179     }// mult>0
1180   }// 64 sectors
1181
1182
1183   Double_t XC = 0;
1184   Double_t YC = 0;
1185   Double_t XA = 0;
1186   Double_t YA = 0;
1187
1188   if(sumZNC!=0){
1189
1190     XC = numXZNC/sumZNC;
1191     YC = numYZNC/sumZNC;
1192     fPsiVZEROC = 0.5*TMath::ATan2(YC-meanYC, XC-meanXC);
1193     if(fPsiVZEROC>TMath::Pi()){fPsiVZEROC-=TMath::Pi();}
1194     if(fPsiVZEROC<0){fPsiVZEROC+=TMath::Pi();}
1195
1196
1197     fh2XYC->Fill(XC-meanXC,YC-meanYC); // control
1198     fp1RPXC->Fill(aod->GetRunNumber(),XC);
1199     fp1RPYC->Fill(aod->GetRunNumber(),YC);
1200
1201
1202   }
1203
1204   if(sumZNA!=0){
1205     XA = numXZNA/sumZNA;
1206     YA = numYZNA/sumZNA;
1207     fPsiVZEROA = 0.5*TMath::ATan2(YA-meanYA, XA-meanXA);
1208     if(fPsiVZEROA>TMath::Pi()){fPsiVZEROA-=TMath::Pi();}
1209     if(fPsiVZEROA<0){fPsiVZEROA+=TMath::Pi();}
1210
1211     fh2XYA->Fill(XA-meanXA,YA-meanYA); // control
1212     fp1RPXA->Fill(aod->GetRunNumber(),XA);
1213     fp1RPYA->Fill(aod->GetRunNumber(),YA);
1214
1215   }
1216   return kTRUE;
1217
1218 }
1219
1220 Int_t  AliAnalysisTaskJetServices::GetListOfTracks(TList *list){
1221   Int_t iCount = 0;
1222   AliAODEvent *aod = 0;
1223   if(fUseAODInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1224   else aod = AODEvent();
1225   if(!aod){
1226     return iCount;
1227   }
1228   for(int it = 0;it < aod->GetNumberOfTracks();++it){
1229     AliAODTrack *tr = aod->GetTrack(it);
1230     if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1231     if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1232     if(tr->Pt()<fMinTrackPt)continue;
1233     list->Add(tr);
1234     iCount++;
1235   }
1236   list->Sort();
1237   return iCount;
1238
1239 }
1240
1241 void AliAnalysisTaskJetServices::SetNTrigger(Int_t n){
1242   if(n>0){
1243       fNTrigger = n;
1244       delete [] fTriggerName;
1245       fTriggerName = new TString [fNTrigger];
1246       delete [] fTriggerBit;fTriggerBit = 0;
1247       fTriggerBit = new UInt_t [fNTrigger];
1248   }
1249   else{
1250     fNTrigger = 0;
1251   }
1252 }
1253
1254 void AliAnalysisTaskJetServices::SetTrigger(Int_t i,UInt_t it,const char* c){
1255   if(i<fNTrigger){
1256     Printf("%d",it);
1257     Printf("%p",c);
1258     Printf("%s",c);
1259     Printf("%p",&fTriggerName[i]);
1260
1261     fTriggerBit[i] = it;
1262     fTriggerName[i] =  c; // placement new
1263     //    new(&fTriggerName[i]) TString(c); // placement new
1264     Printf("%s",fTriggerName[i].Data());
1265
1266   }
1267 }