PWGLF/SPECTRA/PiKaPr/TestAOD/AliAnalysisTaskV2AllChAOD.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //-----------------------------------------------------------------
  17 //         AliAnalysisTaskV2AllChAOD class
  18 //-----------------------------------------------------------------
  19
  20 #include "TChain.h"
  21 #include "TTree.h"
  22 #include "TLegend.h"
  23 #include "TH1F.h"
  24 #include "TH1D.h"
  25 #include "TH2F.h"
  26 #include "THnSparse.h"
  27 #include "TProfile.h"
  28 #include "TCanvas.h"
  29 #include "AliAnalysisTask.h"
  30 #include "AliAODTrack.h"
  31 #include "AliAODMCParticle.h"
  32 #include "AliVParticle.h"
  33 #include "AliAODEvent.h"
  34 #include "AliAODInputHandler.h"
  35 #include "AliAnalysisTaskV2AllChAOD.h"
  36 #include "AliAnalysisTaskESDfilter.h"
  37 #include "AliAnalysisDataContainer.h"
  38 #include "AliSpectraAODTrackCuts.h"
  39 #include "AliSpectraAODEventCuts.h"
  40 #include "AliPIDCombined.h"
  41 #include "AliCentrality.h"
  42 #include "TProof.h"
  43 #include "AliVEvent.h"
  44 #include "AliStack.h"
  45 #include <TMCProcess.h>
  46 #include <TRandom.h>
  47
  48 #include <iostream>
  49
  50 using namespace std;
  51
  52 ClassImp(AliAnalysisTaskV2AllChAOD)
  53
  54 //________________________________________________________________________
  55 AliAnalysisTaskV2AllChAOD::AliAnalysisTaskV2AllChAOD(const char *name) : AliAnalysisTaskSE(name),
  56 fAOD(0x0),
  57 fTrackCuts(0x0),
  58 fEventCuts(0x0),
  59 fIsMC(0),
  60 fCharge(0),
  61 fVZEROside(0),
  62 fOutput(0x0),
  63 fOutput_lq(0x0),
  64 fOutput_sq(0x0),
  65 fnCentBins(20),
  66 fnQvecBins(100),
  67 fQvecUpperLim(100),
  68 fCutLargeQperc(90.),
  69 fCutSmallQperc(10.),
  70 fEtaGapMin(-0.5),
  71 fEtaGapMax(0.5),
  72 fTrkBit(128),
  73 fEtaCut(0.8),
  74 fMinPt(0),
  75 fMaxPt(20.0),
  76 fMinTPCNcls(70),
  77 fFillTHn(kFALSE),
  78 fCentrality(0),
  79 fQvector(0),
  80 fQvector_lq(0),
  81 fQvector_sq(0),
  82 fResSP(0),
  83 fResSP_vs_Cent(0),
  84 f2partCumQA_vs_Cent(0),
  85 f2partCumQB_vs_Cent(0),
  86 fEta_vs_Phi_bef(0),
  87 fEta_vs_PhiA(0),
  88 fEta_vs_PhiB(0),
  89 fResSP_lq(0),
  90 fResSP_vs_Cent_lq(0),
  91 f2partCumQA_vs_Cent_lq(0),
  92 f2partCumQB_vs_Cent_lq(0),
  93 fResSP_sq(0),
  94 fResSP_vs_Cent_sq(0),
  95 f2partCumQA_vs_Cent_sq(0),
  96 f2partCumQB_vs_Cent_sq(0),
  97 fResSP_inclusive(0),
  98 fv2SPGap1A_inclusive_mb(0),
  99 fv2SPGap1B_inclusive_mb(0),
 100 fv2SPGap1A_inclusive_lq(0),
 101 fv2SPGap1B_inclusive_lq(0),
 102 fv2SPGap1A_inclusive_sq(0),
 103 fv2SPGap1B_inclusive_sq(0),
 104 fResSPmc_inclusive(0),
 105 fv2SPGap1Amc_inclusive_mb(0),
 106 fv2SPGap1Bmc_inclusive_mb(0),
 107 fv2SPGap1Amc_inclusive_lq(0),
 108 fv2SPGap1Bmc_inclusive_lq(0),
 109 fv2SPGap1Amc_inclusive_sq(0),
 110 fv2SPGap1Bmc_inclusive_sq(0),
 111 fResGap1w(0),
 112 fV2IntGap1w(0),
 113 fIsRecoEff(0),
 114 fRecoEffList(0),
 115 fQvecGen(0),
 116 fQgenType(0),
 117 fnNchBins(400),
 118 fDoCentrSystCentrality(0)
 119 {
 120
 121   for (Int_t i = 0; i< 9; i++){
 122
 123     fv2SPGap1A[i] = 0;
 124     fv2SPGap1B[i] = 0;
 125
 126     fSinGap1Aq[i] = 0;
 127     fCosGap1Aq[i] = 0;
 128     fSinGap1Bq[i] = 0;
 129     fCosGap1Bq[i] = 0;
 130
 131     fSinGap1A[i] = 0;
 132     fCosGap1A[i] = 0;
 133     fSinGap1B[i] = 0;
 134     fCosGap1B[i] = 0;
 135
 136     //large q
 137     fv2SPGap1A_lq[i] = 0;
 138     fv2SPGap1B_lq[i] = 0;
 139
 140     fSinGap1Aq_lq[i] = 0;
 141     fCosGap1Aq_lq[i] = 0;
 142     fSinGap1Bq_lq[i] = 0;
 143     fCosGap1Bq_lq[i] = 0;
 144
 145     fSinGap1A_lq[i] = 0;
 146     fCosGap1A_lq[i] = 0;
 147     fSinGap1B_lq[i] = 0;
 148     fCosGap1B_lq[i] = 0;
 149
 150     //small q
 151     fv2SPGap1A_sq[i] = 0;
 152     fv2SPGap1B_sq[i] = 0;
 153
 154     fSinGap1Aq_sq[i] = 0;
 155     fCosGap1Aq_sq[i] = 0;
 156     fSinGap1Bq_sq[i] = 0;
 157     fCosGap1Bq_sq[i] = 0;
 158
 159     fSinGap1A_sq[i] = 0;
 160     fCosGap1A_sq[i] = 0;
 161     fSinGap1B_sq[i] = 0;
 162     fCosGap1B_sq[i] = 0;
 163
 164     fResSP_vs_Qvec[i] = 0;
 165     fV2IntGap1wq[i] = 0;
 166
 167   }
 168
 169   fRecoEffList=new TList();
 170   fRecoEffList->SetOwner();
 171   fRecoEffList->SetName("fRecoEffList");
 172
 173   // Default constructor
 174   DefineInput(0, TChain::Class());
 175   DefineOutput(1, TList::Class());
 176   DefineOutput(2, AliSpectraAODEventCuts::Class());
 177   DefineOutput(3, AliSpectraAODTrackCuts::Class());
 178   DefineOutput(4, TList::Class());
 179   DefineOutput(5, TList::Class());
 180 }
 181
 182 //________________________________________________________________________
 183 void AliAnalysisTaskV2AllChAOD::UserCreateOutputObjects()
 184 {
 185   // create output objects
 186   fOutput=new TList();
 187   fOutput->SetOwner();
 188   fOutput->SetName("fOutput");
 189
 190   fOutput_lq=new TList();
 191   fOutput_lq->SetOwner();
 192   fOutput_lq->SetName("fOutput_lq");
 193
 194   fOutput_sq=new TList();
 195   fOutput_sq->SetOwner();
 196   fOutput_sq->SetName("fOutput_sq");
 197
 198   if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
 199   if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
 200
 201   if( fFillTHn ){
 202     //dimensions of THnSparse for Q vector checks
 203     const Int_t nvarev=6;
 204     //                                             cent         q-rec_perc        qvec_v0a       q-rec_v0c        qvec-gen_tpc          Nch
 205     Int_t    binsHistRealEv[nvarev] = {     fnCentBins,              100,        fnQvecBins,     fnQvecBins,       fnQvecBins,     fnNchBins};
 206     Double_t xminHistRealEv[nvarev] = {             0.,               0.,                0.,             0.,               0.,            0.};
 207     Double_t xmaxHistRealEv[nvarev] = {           100.,             100.,     fQvecUpperLim,  fQvecUpperLim,    fQvecUpperLim,         2000.};
 208
 209     THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
 210     NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
 211     NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
 212
 213     NSparseHistEv->GetAxis(1)->SetTitle("q-vec rec percentile");
 214     NSparseHistEv->GetAxis(1)->SetName("Qrec_perc");
 215
 216     NSparseHistEv->GetAxis(2)->SetTitle("q-vec V0A");
 217     NSparseHistEv->GetAxis(2)->SetName("Qrec_V0A");
 218
 219     NSparseHistEv->GetAxis(3)->SetTitle("q-vec V0C");
 220     NSparseHistEv->GetAxis(3)->SetName("Qrec_V0C");
 221
 222     NSparseHistEv->GetAxis(4)->SetTitle("q-vec TPC");
 223     NSparseHistEv->GetAxis(4)->SetName("Qgen_TPC");
 224
 225     NSparseHistEv->GetAxis(5)->SetTitle("Ncharged");
 226     NSparseHistEv->GetAxis(5)->SetName("Nch");
 227     fOutput->Add(NSparseHistEv);
 228   }
 229
 230   fCentrality = new TH1D("fCentrality", "centrality distribution; centrality", 200, 0., 100);
 231   fOutput->Add(fCentrality);
 232
 233   fQvector = new TH1D("fQvector", "q-vector distribution; q-vector", fnQvecBins, 0., fQvecUpperLim);
 234   fOutput->Add(fQvector);
 235
 236   fQvector_lq = new TH1D("fQvector_lq", "q-vector distribution; q-vector", fnQvecBins, 0., fQvecUpperLim);
 237   fOutput_lq->Add(fQvector_lq);
 238
 239   fQvector_sq = new TH1D("fQvector_sq", "q-vector distribution; q-vector", fnQvecBins, 0., fQvecUpperLim);
 240   fOutput_sq->Add(fQvector_sq);
 241
 242   // binning common to all the THn
 243   //change it according to your needs + move it to global variables -> setter/getter
 244   //   Double_t ptBins[] = {0., 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.4, 2.8, 3.2, 3.6, 4.0, 4.5, 5.0, 5.5, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0, 20.0};
 245   //   const Int_t nptBins = 31;
 246   Double_t ptBins[] = {0., 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 3.0, 3.4, 3.8, 4.2, 4.6, 5.0, 5.5, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0, 20.0};
 247   const Int_t nptBins = 33;
 248
 249   fResSP = new TProfile("fResSP", "Resolution; centrality; Resolution", 9, -0.5, 8.5);
 250   fOutput->Add(fResSP);
 251
 252   fResSP_vs_Cent = new TProfile("fResSP_vs_Cent", "Resolution; centrality; Resolution", 20., 0., 100.);
 253   fOutput->Add(fResSP_vs_Cent);
 254
 255   f2partCumQA_vs_Cent = new TProfile("f2partCumQA_vs_Cent", "Resolution; centrality; Resolution", 100., 0., 100.);
 256   fOutput->Add(f2partCumQA_vs_Cent);
 257
 258   f2partCumQB_vs_Cent = new TProfile("f2partCumQB_vs_Cent", "Resolution; centrality; Resolution", 100., 0., 100.);
 259   fOutput->Add(f2partCumQB_vs_Cent);
 260
 261   fEta_vs_Phi_bef = new TH2D("fEta_vs_Phi_bef","eta vs phi distribution before eta gap;#eta;#phi",200.,-1.,1.,175.,0.,7.);
 262   fOutput->Add(fEta_vs_Phi_bef);
 263
 264   fEta_vs_PhiA = new TH2D("fEta_vs_PhiA","eta vs phi distribution;#eta;#phi",200.,-1.,1.,175.,0.,7.);
 265   fOutput->Add(fEta_vs_PhiA);
 266
 267   fEta_vs_PhiB = new TH2D("fEta_vs_PhiB","eta vs phi distribution;#eta;#phi",200.,-1.,1.,175.,0.,7.);
 268   fOutput->Add(fEta_vs_PhiB);
 269
 270   // MC closure test
 271   fResSP_inclusive = new TProfile("fResSP_inclusive", "Resolution; ese; Resolution", 3, 0., 3.);
 272   fOutput->Add(fResSP_inclusive);
 273
 274   fv2SPGap1A_inclusive_mb = new TProfile("fv2SPGap1A_inclusive_mb", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
 275   fOutput->Add(fv2SPGap1A_inclusive_mb);
 276
 277   fv2SPGap1B_inclusive_mb = new TProfile("fv2SPGap1B_inclusive_mb", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
 278   fOutput->Add(fv2SPGap1B_inclusive_mb);
 279
 280
 281   //large q
 282   fResSP_lq = new TProfile("fResSP_lq", "Resolution; centrality; Resolution", 9, -0.5, 8.5);
 283   fOutput_lq->Add(fResSP_lq);
 284
 285   fResSP_vs_Cent_lq = new TProfile("fResSP_vs_Cent_lq", "Resolution; centrality; Resolution", 20., 0., 100.);
 286   fOutput_lq->Add(fResSP_vs_Cent_lq);
 287
 288   f2partCumQA_vs_Cent_lq = new TProfile("f2partCumQA_vs_Cent_lq", "Resolution; centrality; Resolution", 100., 0., 100.);
 289   fOutput_lq->Add(f2partCumQA_vs_Cent_lq);
 290
 291   f2partCumQB_vs_Cent_lq = new TProfile("f2partCumQB_vs_Cent_lq", "Resolution; centrality; Resolution", 100., 0., 100.);
 292   fOutput_lq->Add(f2partCumQB_vs_Cent_lq);
 293
 294   // MC closure test
 295   fv2SPGap1A_inclusive_lq = new TProfile("fv2SPGap1A_inclusive_lq", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
 296   fOutput_lq->Add(fv2SPGap1A_inclusive_lq);
 297
 298   fv2SPGap1B_inclusive_lq = new TProfile("fv2SPGap1B_inclusive_lq", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
 299   fOutput_lq->Add(fv2SPGap1B_inclusive_lq);
 300
 301   //small q resolution
 302   fResSP_sq = new TProfile("fResSP_sq", "Resolution; centrality; Resolution", 9, -0.5, 8.5);
 303   fOutput_sq->Add(fResSP_sq);
 304
 305   fResSP_vs_Cent_sq = new TProfile("fResSP_vs_Cent_sq", "Resolution; centrality; Resolution", 20., 0., 100.);
 306   fOutput_sq->Add(fResSP_vs_Cent_sq);
 307
 308   f2partCumQA_vs_Cent_sq = new TProfile("f2partCumQA_vs_Cent_sq", "Resolution; centrality; Resolution", 100., 0., 100.);
 309   fOutput_sq->Add(f2partCumQA_vs_Cent_sq);
 310
 311   f2partCumQB_vs_Cent_sq = new TProfile("f2partCumQB_vs_Cent_sq", "Resolution; centrality; Resolution", 100., 0., 100.);
 312   fOutput_sq->Add(f2partCumQB_vs_Cent_sq);
 313
 314   // MC closure test
 315   fv2SPGap1A_inclusive_sq = new TProfile("fv2SPGap1A_inclusive_sq", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
 316   fOutput_sq->Add(fv2SPGap1A_inclusive_sq);
 317
 318   fv2SPGap1B_inclusive_sq = new TProfile("fv2SPGap1B_inclusive_sq", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
 319   fOutput_sq->Add(fv2SPGap1B_inclusive_sq);
 320
 321   for (Int_t iC = 0; iC < 9; iC++){
 322
 323     fv2SPGap1A[iC] = new TProfile(Form("fv2SPGap1A_%d", iC), "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
 324     fOutput->Add(fv2SPGap1A[iC]);
 325
 326     fv2SPGap1B[iC] = new TProfile(Form("fv2SPGap1B_%d", iC), "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
 327     fOutput->Add(fv2SPGap1B[iC]);
 328
 329     fSinGap1Aq[iC] = new TProfile(Form("fSinGap1Aq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
 330     fOutput->Add(fSinGap1Aq[iC]);
 331
 332     fCosGap1Aq[iC] = new TProfile(Form("fCosGap1Aq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
 333     fOutput->Add(fCosGap1Aq[iC]);
 334
 335     fSinGap1Bq[iC] = new TProfile(Form("fSinGap1Bq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
 336     fOutput->Add(fSinGap1Bq[iC]);
 337
 338     fCosGap1Bq[iC] = new TProfile(Form("fCosGap1Bq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
 339     fOutput->Add(fCosGap1Bq[iC]);
 340
 341     fSinGap1A[iC] = new TProfile(Form("fSinGap1A_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
 342     fOutput->Add(fSinGap1A[iC]);
 343
 344     fCosGap1A[iC] = new TProfile(Form("fCosGap1A_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
 345     fOutput->Add(fCosGap1A[iC]);
 346
 347     fSinGap1B[iC] = new TProfile(Form("fSinGap1B_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
 348     fOutput->Add(fSinGap1B[iC]);
 349
 350     fCosGap1B[iC] = new TProfile(Form("fCosGap1B_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
 351     fOutput->Add(fCosGap1B[iC]);
 352
 353     //large q
 354     fv2SPGap1A_lq[iC] = new TProfile(Form("fv2SPGap1A_lq_%d", iC), "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
 355     fOutput_lq->Add(fv2SPGap1A_lq[iC]);
 356
 357     fv2SPGap1B_lq[iC] = new TProfile(Form("fv2SPGap1B_lq_%d", iC), "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
 358     fOutput_lq->Add(fv2SPGap1B_lq[iC]);
 359
 360     fSinGap1Aq_lq[iC] = new TProfile(Form("fSinGap1Aq_lq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
 361     fOutput_lq->Add(fSinGap1Aq_lq[iC]);
 362
 363     fCosGap1Aq_lq[iC] = new TProfile(Form("fCosGap1Aq_lq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
 364     fOutput_lq->Add(fCosGap1Aq_lq[iC]);
 365
 366     fSinGap1Bq_lq[iC] = new TProfile(Form("fSinGap1Bq_lq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
 367     fOutput_lq->Add(fSinGap1Bq_lq[iC]);
 368
 369     fCosGap1Bq_lq[iC] = new TProfile(Form("fCosGap1Bq_lq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
 370     fOutput_lq->Add(fCosGap1Bq_lq[iC]);
 371
 372     fSinGap1A_lq[iC] = new TProfile(Form("fSinGap1A_lq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
 373     fOutput_lq->Add(fSinGap1A_lq[iC]);
 374
 375     fCosGap1A_lq[iC] = new TProfile(Form("fCosGap1A_lq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
 376     fOutput_lq->Add(fCosGap1A_lq[iC]);
 377
 378     fSinGap1B_lq[iC] = new TProfile(Form("fSinGap1B_lq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
 379     fOutput_lq->Add(fSinGap1B_lq[iC]);
 380
 381     fCosGap1B_lq[iC] = new TProfile(Form("fCosGap1B_lq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
 382     fOutput_lq->Add(fCosGap1B_lq[iC]);
 383
 384     //small q
 385     fv2SPGap1A_sq[iC] = new TProfile(Form("fv2SPGap1A_sq_%d", iC), "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
 386     fOutput_sq->Add(fv2SPGap1A_sq[iC]);
 387
 388     fv2SPGap1B_sq[iC] = new TProfile(Form("fv2SPGap1B_sq_%d", iC), "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
 389     fOutput_sq->Add(fv2SPGap1B_sq[iC]);
 390
 391     fSinGap1Aq_sq[iC] = new TProfile(Form("fSinGap1Aq_sq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
 392     fOutput_sq->Add(fSinGap1Aq_sq[iC]);
 393
 394     fCosGap1Aq_sq[iC] = new TProfile(Form("fCosGap1Aq_sq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
 395     fOutput_sq->Add(fCosGap1Aq_sq[iC]);
 396
 397     fSinGap1Bq_sq[iC] = new TProfile(Form("fSinGap1Bq_sq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
 398     fOutput_sq->Add(fSinGap1Bq_sq[iC]);
 399
 400     fCosGap1Bq_sq[iC] = new TProfile(Form("fCosGap1Bq_sq_%d", iC), "p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
 401     fOutput_sq->Add(fCosGap1Bq_sq[iC]);
 402
 403     fSinGap1A_sq[iC] = new TProfile(Form("fSinGap1A_sq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
 404     fOutput_sq->Add(fSinGap1A_sq[iC]);
 405
 406     fCosGap1A_sq[iC] = new TProfile(Form("fCosGap1A_sq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
 407     fOutput_sq->Add(fCosGap1A_sq[iC]);
 408
 409     fSinGap1B_sq[iC] = new TProfile(Form("fSinGap1B_sq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
 410     fOutput_sq->Add(fSinGap1B_sq[iC]);
 411
 412     fCosGap1B_sq[iC] = new TProfile(Form("fCosGap1B_sq_%d", iC), "p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
 413     fOutput_sq->Add(fCosGap1B_sq[iC]);
 414
 415     //v2 vs qvec...
 416     fResSP_vs_Qvec[iC] = new TProfile(Form("fResSP_vs_Qvec_%d", iC), "Resolution; Qvec (V0A); Resolution", 20., 0., 100.);
 417     fResSP_vs_Qvec[iC]->Sumw2();
 418     fOutput->Add(fResSP_vs_Qvec[iC]);
 419
 420     fV2IntGap1wq[iC] = new TProfile(Form("fV2IntGap1wq_%d", iC), "integrated v_{2} vs q-vector; Qvec (V0A); v_{2}", 20., 0., 100.);
 421     fV2IntGap1wq[iC]->Sumw2();
 422     fOutput->Add(fV2IntGap1wq[iC]);
 423   };
 424
 425   fResGap1w = new TProfile("fResGap1w", "Resolution; centrality; Resolution", 9, -0.5, 8.5);
 426   fResGap1w->Sumw2();
 427   fOutput->Add(fResGap1w);
 428
 429   fV2IntGap1w = new TProfile("fV2IntGap1w", "; centrality; v_{2}", 9, -0.5, 8.5);
 430   fV2IntGap1w->Sumw2();
 431   fOutput->Add(fV2IntGap1w);
 432
 433   if(fIsMC){
 434     fResSPmc_inclusive = new TProfile("fResSPmc_inclusive", "Resolution; ese; Resolution", 3, 0., 3.);
 435     fOutput->Add(fResSPmc_inclusive);
 436
 437     fv2SPGap1Amc_inclusive_mb = new TProfile("fv2SPGap1Amc_inclusive_mb", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
 438     fOutput->Add(fv2SPGap1Amc_inclusive_mb);
 439
 440     fv2SPGap1Bmc_inclusive_mb = new TProfile("fv2SPGap1Bmc_inclusive_mb", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
 441     fOutput->Add(fv2SPGap1Bmc_inclusive_mb);
 442
 443     //large-q
 444     fv2SPGap1Amc_inclusive_lq = new TProfile("fv2SPGap1Amc_inclusive_lq", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
 445     fOutput_lq->Add(fv2SPGap1Amc_inclusive_lq);
 446
 447     fv2SPGap1Bmc_inclusive_lq = new TProfile("fv2SPGap1Bmc_inclusive_lq", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
 448     fOutput_lq->Add(fv2SPGap1Bmc_inclusive_lq);
 449
 450     //small-q
 451     fv2SPGap1Amc_inclusive_sq = new TProfile("fv2SPGap1Amc_inclusive_sq", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
 452     fOutput_sq->Add(fv2SPGap1Amc_inclusive_sq);
 453
 454     fv2SPGap1Bmc_inclusive_sq = new TProfile("fv2SPGap1Bmc_inclusive_sq", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
 455     fOutput_sq->Add(fv2SPGap1Bmc_inclusive_sq);
 456   }
 457
 458   PostData(1, fOutput  );
 459   PostData(2, fEventCuts);
 460   PostData(3, fTrackCuts);
 461   PostData(4, fOutput_lq  );
 462   PostData(5, fOutput_sq  );
 463 }
 464
 465 //________________________________________________________________________
 466
 467 void AliAnalysisTaskV2AllChAOD::UserExec(Option_t *)
 468 {
 469   //Printf("An event");
 470   // main event loop
 471   fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
 472   if (!fAOD) {
 473     AliWarning("ERROR: AliAODEvent not available \n");
 474     return;
 475   }
 476
 477   if (strcmp(fAOD->ClassName(), "AliAODEvent"))
 478   {
 479     AliFatal("Not processing AODs");
 480   }
 481
 482   if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection
 483
 484   //Get q-vector percentile.
 485   Double_t Qvec=0.;
 486   if(fIsMC && fQvecGen) Qvec = fEventCuts->GetQvecPercentileMC(fVZEROside, fQgenType);
 487   else Qvec = fEventCuts->GetQvecPercentile(fVZEROside);
 488
 489   fQvector->Fill(Qvec);
 490   if (Qvec > fCutLargeQperc && Qvec < 100.) fQvector_lq->Fill(Qvec);
 491   if (Qvec > 0. && Qvec < fCutSmallQperc) fQvector_sq->Fill(Qvec);
 492
 493   Double_t Cent=(fDoCentrSystCentrality)?1.01*fEventCuts->GetCent():fEventCuts->GetCent();
 494   fCentrality->Fill(Cent);
 495
 496   Int_t centV0 = -1;
 497   if ((Cent > 0) && (Cent <= 5.0))
 498     centV0 = 0;
 499   else if ((Cent > 5.0) && (Cent <= 10.0))
 500     centV0 = 1;
 501   else if ((Cent > 10.0) && (Cent <= 20.0))
 502     centV0 = 2;
 503   else if ((Cent > 20.0) && (Cent <= 30.0))
 504     centV0 = 3;
 505   else if ((Cent > 30.0) && (Cent <= 40.0))
 506     centV0 = 4;
 507   else if ((Cent > 40.0) && (Cent <= 50.0))
 508     centV0 = 5;
 509   else if ((Cent > 50.0) && (Cent <= 60.0))
 510     centV0 = 6;
 511   else if ((Cent > 60.0) && (Cent <= 70.0))
 512     centV0 = 7;
 513   else if ((Cent > 70.0) && (Cent <= 80.0))
 514     centV0 = 8;
 515
 516   if(centV0==-1)return; // FIXME if the centrality is not defined or >80%... return!!!
 517
 518   if(fIsMC) MCclosure(Qvec); // fill mc histograms for montecarlo closure
 519
 520   Double_t QxGap1A = 0., QyGap1A = 0.;
 521   Double_t QxGap1B = 0., QyGap1B = 0.;
 522   Double_t multGap1A = 0, multGap1B = 0;
 523
 524   for (Int_t loop = 0; loop < 2; loop++){
 525
 526     //main loop on tracks
 527     for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
 528       AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
 529       if(!track) AliFatal("Not a standard AOD");
 530       if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge
 531       if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)
 532
 533       fEta_vs_Phi_bef->Fill( track->Eta(), track->Phi() );
 534
 535       //       if (fIsRecoEff){
 536       //
 537       //        // 2) reject randomly tracks at high pT until the reconstruction efficiency becomes flat (add the following before the loop == 0 part): (mail by Alexandru)
 538       //
 539       //        Double_t recoEff = GetRecoEff(track->Pt(), centV0);
 540       //         if (recoEff < 0){
 541       //          cout<<"No reconstruction efficiency!"<<endl;
 542       //           continue;
 543       //        }
 544       //
 545       //        Double_t rndPt = gRandom->Rndm();
 546       // //         cout<<"rndPt: "<<rndPt<<endl;
 547
 548       //        Double_t minRecPt = GetRecoEff(0.200001, centV0);
 549       // //         cout<<"minRecPt: "<<minRecPt<<endl;
 550
 551       //        if (rndPt > minRecPt/recoEff){
 552       // //       cout<<"Track rejected: "<<iTracks<<"  from "<<fAOD->GetNumberOfTracks()<<endl;
 553       //          continue;
 554       //        }
 555       //
 556       //       } // end fIsRecoEff
 557
 558       if (loop == 0) {
 559
 560         if (track->Eta() > fEtaGapMax){
 561
 562           if(fIsRecoEff){
 563             Double_t receff = GetRecoEff(track->Pt(), centV0);
 564             QxGap1A += (TMath::Cos(2.*track->Phi()))/receff;
 565             QyGap1A += (TMath::Sin(2.*track->Phi()))/receff;
 566             multGap1A+=1./receff;
 567           } else  {
 568             QxGap1A += TMath::Cos(2.*track->Phi());
 569             QyGap1A += TMath::Sin(2.*track->Phi());
 570             multGap1A++;
 571           }
 572
 573           fSinGap1Aq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
 574           fCosGap1Aq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
 575
 576           fEta_vs_PhiA->Fill( track->Eta(), track->Phi() );
 577
 578           if (Qvec > fCutLargeQperc && Qvec < 100.){
 579             fSinGap1Aq_lq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
 580             fCosGap1Aq_lq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
 581           }
 582
 583           if (Qvec > 0. && Qvec < fCutSmallQperc){
 584             fSinGap1Aq_sq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
 585             fCosGap1Aq_sq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
 586           }
 587
 588         }
 589
 590         if (track->Eta() < fEtaGapMin){
 591
 592           if(fIsRecoEff){
 593             Double_t receff = GetRecoEff(track->Pt(), centV0);
 594             QxGap1B += (TMath::Cos(2.*track->Phi()))/receff;
 595             QyGap1B += (TMath::Sin(2.*track->Phi()))/receff;
 596             multGap1B+=1./receff;
 597           } else {
 598             QxGap1B += TMath::Cos(2.*track->Phi());
 599             QyGap1B += TMath::Sin(2.*track->Phi());
 600             multGap1B++;
 601           }
 602
 603           fCosGap1Bq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
 604           fSinGap1Bq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
 605
 606           fEta_vs_PhiB->Fill( track->Eta(), track->Phi() );
 607
 608           if (Qvec > fCutLargeQperc && Qvec < 100.){
 609             fSinGap1Bq_lq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
 610             fCosGap1Bq_lq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
 611           }
 612
 613           if (Qvec > 0. && Qvec < fCutSmallQperc){
 614             fSinGap1Bq_sq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
 615             fCosGap1Bq_sq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
 616           }
 617
 618         }
 619
 620       } else {
 621
 622         //eval v2 scalar product
 623         if (track->Eta() < fEtaGapMin && multGap1A > 0){
 624
 625           Double_t v2SPGap1A = 0.;
 626           if(fIsRecoEff){
 627             Double_t receff = GetRecoEff(track->Pt(), centV0);
 628             Double_t uxGap1A = (TMath::Cos(2.*track->Phi()))/receff;
 629             Double_t uyGap1A = (TMath::Sin(2.*track->Phi()))/receff;
 630             //      multGap1A = multGap1A/receff;
 631             v2SPGap1A = (uxGap1A*QxGap1A + uyGap1A*QyGap1A)/(Double_t)multGap1A;
 632           } else v2SPGap1A = (TMath::Cos(2.*track->Phi())*QxGap1A + TMath::Sin(2.*track->Phi())*QyGap1A)/(Double_t)multGap1A;
 633
 634           //           Double_t v2SPGap1A = (TMath::Cos(2.*track->Phi())*QxGap1A + TMath::Sin(2.*track->Phi())*QyGap1A)/(Double_t)multGap1A;
 635
 636           fv2SPGap1A[centV0]->Fill(track->Pt(), v2SPGap1A);
 637
 638           fv2SPGap1A_inclusive_mb->Fill(track->Pt(), v2SPGap1A); //mb v2 for mc closure
 639
 640           fSinGap1A[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
 641           fCosGap1A[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
 642
 643           if (Qvec > fCutLargeQperc && Qvec < 100.){
 644             fv2SPGap1A_lq[centV0]->Fill(track->Pt(), v2SPGap1A);
 645             fSinGap1A_lq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
 646             fCosGap1A_lq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
 647
 648             fv2SPGap1A_inclusive_lq->Fill(track->Pt(), v2SPGap1A); //lq v2 for mc closure
 649           }
 650
 651           if (Qvec > 0. && Qvec < fCutSmallQperc){
 652             fv2SPGap1A_sq[centV0]->Fill(track->Pt(), v2SPGap1A);
 653             fSinGap1A_sq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
 654             fCosGap1A_sq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
 655
 656             fv2SPGap1A_inclusive_sq->Fill(track->Pt(), v2SPGap1A); //sq v2 for mc closure
 657           }
 658
 659         }
 660
 661         if (track->Eta() > fEtaGapMax && multGap1B > 0){
 662
 663           Double_t v2SPGap1B = 0;
 664           if(fIsRecoEff){
 665             Double_t receff = GetRecoEff(track->Pt(), centV0);
 666             Double_t uxGap1B = (TMath::Cos(2.*track->Phi()))/receff;
 667             Double_t uyGap1B = (TMath::Sin(2.*track->Phi()))/receff;
 668 //            multGap1B = multGap1B/receff;
 669             v2SPGap1B = (uxGap1B*QxGap1B + uyGap1B*QyGap1B)/(Double_t)multGap1B;
 670           }
 671           else v2SPGap1B = (TMath::Cos(2.*track->Phi())*QxGap1B + TMath::Sin(2.*track->Phi())*QyGap1B)/(Double_t)multGap1B;
 672
 673           //      Double_t v2SPGap1B = (TMath::Cos(2.*track->Phi())*QxGap1B + TMath::Sin(2.*track->Phi())*QyGap1B)/(Double_t)multGap1B;
 674
 675           fv2SPGap1B[centV0]->Fill(track->Pt(), v2SPGap1B);
 676
 677           fv2SPGap1B_inclusive_mb->Fill(track->Pt(), v2SPGap1B); //mb v2 for mc closure
 678
 679           fCosGap1B[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
 680           fSinGap1B[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
 681
 682           if (Qvec > fCutLargeQperc && Qvec < 100.){
 683             fv2SPGap1B_lq[centV0]->Fill(track->Pt(), v2SPGap1B);
 684             fSinGap1B_lq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
 685             fCosGap1B_lq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
 686
 687             fv2SPGap1B_inclusive_lq->Fill(track->Pt(), v2SPGap1B); //lq v2 for mc closure
 688           }
 689
 690           if (Qvec > 0. && Qvec < fCutSmallQperc){
 691             fv2SPGap1B_sq[centV0]->Fill(track->Pt(), v2SPGap1B);
 692             fSinGap1B_sq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
 693             fCosGap1B_sq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
 694
 695             fv2SPGap1B_inclusive_sq->Fill(track->Pt(), v2SPGap1B); //sq v2 for mc closure
 696           }
 697
 698         }
 699       }// end else
 700     } // end loop on tracks
 701   } // end loop
 702
 703
 704   if (multGap1A > 0 && multGap1B > 0){
 705     Double_t res = (QxGap1A*QxGap1B + QyGap1A*QyGap1B)/(Double_t)multGap1A/(Double_t)multGap1B;
 706     fResSP->Fill((Double_t)centV0, res);
 707
 708     fResSP_inclusive->Fill(0., res); //mb v2 for mc closure
 709     fResSP_vs_Cent->Fill(Cent, res);
 710
 711     Double_t f2partCumQA = -999.;
 712     if(multGap1A>1)
 713       f2partCumQA = ( ( (QxGap1A*QxGap1A + QyGap1A*QyGap1A) - (Double_t)multGap1A ) / ((Double_t)multGap1A*((Double_t)multGap1A-1)) );
 714     if(f2partCumQA>0)f2partCumQA_vs_Cent->Fill((Double_t)Cent,f2partCumQA);
 715
 716     Double_t f2partCumQB = -999.;
 717     if(multGap1B>1)
 718       f2partCumQB = ( ( (QxGap1B*QxGap1B + QyGap1B*QyGap1B) - (Double_t)multGap1B ) / ((Double_t)multGap1B*((Double_t)multGap1B-1)) );
 719     if(f2partCumQB>0)f2partCumQB_vs_Cent->Fill((Double_t)Cent,f2partCumQB);
 720
 721     if (Qvec > fCutLargeQperc && Qvec < 100.){
 722       fResSP_lq->Fill((Double_t)centV0, res);
 723       fResSP_vs_Cent_lq->Fill(Cent, res);
 724       if(f2partCumQA>0)f2partCumQA_vs_Cent_lq->Fill((Double_t)Cent,f2partCumQA);
 725       if(f2partCumQB>0)f2partCumQB_vs_Cent_lq->Fill((Double_t)Cent,f2partCumQB);
 726
 727       fResSP_inclusive->Fill(1., res); //lq v2 for mc closure
 728     }
 729
 730     if (Qvec > 0. && Qvec < fCutSmallQperc){
 731       fResSP_sq->Fill((Double_t)centV0, res);
 732       fResSP_vs_Cent_sq->Fill(Cent, res);
 733       if(f2partCumQA>0)f2partCumQA_vs_Cent_sq->Fill((Double_t)Cent,f2partCumQA);
 734       if(f2partCumQB>0)f2partCumQB_vs_Cent_sq->Fill((Double_t)Cent,f2partCumQB);
 735
 736       fResSP_inclusive->Fill(2., res); //sq v2 for mc closure
 737     }
 738   }// end multiplicity if
 739
 740   //v2 vs qvec
 741   if ((multGap1A > 0) && (multGap1B > 0)){
 742
 743     fResGap1w->Fill(Double_t(centV0), (QxGap1A*QxGap1B + QyGap1A*QyGap1B)/multGap1A/multGap1B, (Double_t)(multGap1A*multGap1B));
 744
 745     Double_t nGap1 = multGap1A*multGap1B;
 746     Double_t qqGap1 = (QxGap1A*QxGap1B + QyGap1A*QyGap1B)/nGap1;
 747
 748     fV2IntGap1w->Fill(Double_t(centV0), qqGap1, nGap1);
 749
 750     fResSP_vs_Qvec[centV0]->Fill(Double_t(Qvec), (QxGap1A*QxGap1B + QyGap1A*QyGap1B)/multGap1A/multGap1B, (Double_t)(multGap1A*multGap1B));
 751     fV2IntGap1wq[centV0]->Fill(Double_t(Qvec), qqGap1, nGap1);
 752
 753   }
 754
 755   if( fFillTHn ){
 756
 757     Double_t varEv[6];
 758     varEv[0]=Cent;
 759
 760     varEv[1]=(Double_t)Qvec; // qvec_rec_perc
 761
 762     varEv[2]=(Double_t)fEventCuts->GetqV0A();
 763
 764     varEv[3]=(Double_t)fEventCuts->GetqV0C();
 765
 766     varEv[4]=(Double_t)fEventCuts->GetqTPC();
 767
 768     varEv[5]=(Double_t)fEventCuts->GetNch(); // Nch
 769
 770     ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
 771
 772   }
 773
 774
 775   PostData(1, fOutput  );
 776   PostData(2, fEventCuts);
 777   PostData(3, fTrackCuts);
 778   PostData(4, fOutput_lq  );
 779   PostData(5, fOutput_sq  );
 780 }
 781
 782 //_________________________________________________________________
 783 Bool_t  AliAnalysisTaskV2AllChAOD::GetDCA(const AliAODTrack* trk, Double_t * p){
 784
 785   //AliAODTrack::DCA(): for newest AOD fTrack->DCA() always gives -999. This should fix.
 786   //FIXME should update EventCuts?
 787   //FIXME add track->GetXYZ(p) method
 788
 789   double xyz[3],cov[3];
 790
 791   if (!trk->GetXYZ(xyz)) { // dca is not stored
 792     AliExternalTrackParam etp;
 793     etp.CopyFromVTrack(trk);
 794     AliVEvent* ev = (AliVEvent*)trk->GetEvent();
 795     if (!ev) {/*printf("Event is not connected to the track\n");*/ return kFALSE;}
 796     if (!etp.PropagateToDCA(ev->GetPrimaryVertex(), ev->GetMagneticField(),999,xyz,cov)) return kFALSE; // failed, track is too far from vertex
 797   }
 798   p[0] = xyz[0];
 799   p[1] = xyz[1];
 800   return kTRUE;
 801
 802 }
 803
 804 //_________________________________________________________________
 805 void  AliAnalysisTaskV2AllChAOD::MCclosure(Double_t qvec){
 806   // First do MC to fill up the MC particle array
 807
 808   TClonesArray *arrayMC = 0;
 809   if (fIsMC)
 810   {
 811     arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 812     if (!arrayMC) {
 813       AliFatal("Error: MC particles branch not found!\n");
 814     }
 815
 816     Double_t QxGap1Amc = 0., QyGap1Amc = 0.;
 817     Double_t QxGap1Bmc = 0., QyGap1Bmc = 0.;
 818     Int_t multGap1Amc = 0, multGap1Bmc = 0;
 819
 820     for (Int_t loop = 0; loop < 2; loop++){
 821
 822       Int_t nMC = arrayMC->GetEntries();
 823
 824       for (Int_t iMC = 0; iMC < nMC; iMC++)
 825       {
 826         AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
 827         if(!partMC->Charge()) continue;//Skip neutrals
 828         if(fCharge != 0 && partMC->Charge()*fCharge < 0.) continue;//if fCharge != 0 only select fCharge
 829
 830         if (!(partMC->IsPhysicalPrimary()))
 831           continue;
 832
 833         if(partMC->Eta()<fTrackCuts->GetEtaMin() || partMC->Eta()>fTrackCuts->GetEtaMax()) continue;
 834
 835         //Printf("a particle");
 836
 837
 838         if (loop == 0) {
 839
 840           if (partMC->Eta() > fEtaGapMax){
 841             QxGap1Amc += TMath::Cos(2.*partMC->Phi());
 842             QyGap1Amc += TMath::Sin(2.*partMC->Phi());
 843             multGap1Amc++;
 844           }
 845
 846           if (partMC->Eta() < fEtaGapMin){
 847             QxGap1Bmc += TMath::Cos(2.*partMC->Phi());
 848             QyGap1Bmc += TMath::Sin(2.*partMC->Phi());
 849             multGap1Bmc++;
 850           }
 851
 852         } else {
 853
 854           //eval v2 scalar product
 855           if (partMC->Eta() < fEtaGapMin && multGap1Amc > 0){
 856             Double_t v2SPGap1Amc = (TMath::Cos(2.*partMC->Phi())*QxGap1Amc + TMath::Sin(2.*partMC->Phi())*QyGap1Amc)/(Double_t)multGap1Amc;
 857             fv2SPGap1Amc_inclusive_mb->Fill(partMC->Pt(), v2SPGap1Amc);
 858
 859             if (qvec > fCutLargeQperc && qvec < 100.){
 860               fv2SPGap1Amc_inclusive_lq->Fill(partMC->Pt(), v2SPGap1Amc);
 861             }
 862
 863             if (qvec > 0. && qvec < fCutSmallQperc){
 864               fv2SPGap1Amc_inclusive_sq->Fill(partMC->Pt(), v2SPGap1Amc);
 865             }
 866
 867           }
 868
 869           if (partMC->Eta() > fEtaGapMax && multGap1Bmc > 0){
 870             Double_t v2SPGap1Bmc = (TMath::Cos(2.*partMC->Phi())*QxGap1Bmc + TMath::Sin(2.*partMC->Phi())*QyGap1Bmc)/(Double_t)multGap1Bmc;
 871             fv2SPGap1Bmc_inclusive_mb->Fill(partMC->Pt(), v2SPGap1Bmc);
 872
 873             if (qvec > fCutLargeQperc && qvec < 100.){
 874               fv2SPGap1Bmc_inclusive_lq->Fill(partMC->Pt(), v2SPGap1Bmc);
 875             }
 876
 877             if (qvec > 0. && qvec < fCutSmallQperc){
 878               fv2SPGap1Bmc_inclusive_sq->Fill(partMC->Pt(), v2SPGap1Bmc);
 879             }
 880
 881           }
 882
 883
 884         }// end else
 885       } // end loop on partMCs
 886     } // end loop
 887
 888     if (multGap1Amc > 0 && multGap1Bmc > 0){
 889       Double_t resmc = (QxGap1Amc*QxGap1Bmc + QyGap1Amc*QyGap1Bmc)/(Double_t)multGap1Amc/(Double_t)multGap1Bmc;
 890       fResSPmc_inclusive->Fill(0.,resmc);
 891
 892       if (qvec > fCutLargeQperc && qvec < 100.){
 893         fResSPmc_inclusive->Fill(1.,resmc);
 894       }
 895
 896       if (qvec > 0. && qvec < fCutSmallQperc){
 897         fResSPmc_inclusive->Fill(2.,resmc);
 898       }
 899
 900     }
 901
 902   }// end if MC
 903 }
 904
 905 //_________________________________________________________________
 906 Double_t AliAnalysisTaskV2AllChAOD::GetRecoEff(Double_t pt, Int_t iC){
 907
 908   if(iC>8) return 1.;
 909
 910   if(pt<0.2 || pt>100.) return 1.;
 911
 912   // // //   //spectra ese binning
 913   // // //   //const Double_t ptBins[] = {0.20,0.30,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.4,2.8,3.2,3.6,4.0,5.0,6.0,7.0,8.0,9.0,10.,12.,15.,20.,25.,30.,35.,40.,50.,75.,100.};
 914   // // //   //const Int_t nptBins=34;
 915
 916   const Double_t fEpsilon=0.000001;
 917   //   cout<<"list: "<<endl;
 918   //   fRecoEffList->ls();
 919   TH1D *h = (TH1D*)fRecoEffList->At(iC);
 920
 921   Int_t bin = h->FindBin(pt);
 922
 923   Double_t lowlim = h->GetBinLowEdge(bin);
 924   Double_t uplim = h->GetBinLowEdge(bin) + h->GetBinWidth(bin);
 925
 926   Double_t eff = 1.;
 927
 928   if( pt>lowlim && pt<uplim ) eff = h->GetBinContent(bin);
 929   if( pt == lowlim ) eff = h->GetBinContent( h->FindBin(pt+fEpsilon) );
 930   if( pt == uplim ) eff = h->GetBinContent( h->FindBin(pt-fEpsilon) );
 931
 932   return eff;
 933
 934 }
 935 //_________________________________________________________________
 936 void   AliAnalysisTaskV2AllChAOD::Terminate(Option_t *)
 937 {
 938   // Terminate
 939 }