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