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