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