PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskNucleiv2.cxx

   1 /**************************************************************************
   2  * Contributors are not mentioned at all.                                 *
   3  *                                                                        *
   4  * Permission to use, copy, modify and distribute this software and its   *
   5  * documentation strictly for non-commercial purposes is hereby granted   *
   6  * without fee, provided that the above copyright noticxse appears in all   *
   7  * copies and that both the copyright notice and this permission notice   *
   8  * appear in the supporting documentation. The authors make no claims     *
   9  * about the suitability of this software for any purpose. It is          *
  10  * provided "as is" without express or implied warranty.                  *
  11  **************************************************************************/
  12 //-----------------------------------------------------------------
  13 //                 AliAnalysisTaskNucleiv2 class
  14 //-----------------------------------------------------------------
  15
  16 class TTree;
  17 class TParticle;
  18 class TVector3;
  19
  20 #include "AliAnalysisManager.h"
  21 #include <AliMCEventHandler.h>
  22 #include <AliMCEvent.h>
  23 #include <AliStack.h>
  24
  25 class AliESDVertex;
  26 class AliAODVertex;
  27 class AliESDv0;
  28 class AliAODv0;
  29 class AliCascadeVertexer;
  30
  31 #include <iostream>
  32
  33 #include "TList.h"
  34 #include "TH1.h"
  35 #include "TH2.h"
  36 #include "TH3.h"
  37 #include "TNtuple.h"
  38 #include "TGraph.h"
  39 #include "TF1.h"
  40 #include "TCanvas.h"
  41 #include "TMath.h"
  42 #include "TChain.h"
  43 #include "Riostream.h"
  44 #include "AliLog.h"
  45 #include "AliCascadeVertexer.h"
  46 #include "AliESDEvent.h"
  47 #include "AliESDtrack.h"
  48 #include "AliExternalTrackParam.h"
  49 #include "AliAODEvent.h"
  50 #include "AliInputEventHandler.h"
  51 #include "AliESDcascade.h"
  52 #include "AliAODcascade.h"
  53 #include "AliAnalysisTaskNucleiv2.h"
  54 #include "AliESDtrackCuts.h"
  55 #include "AliCentrality.h"
  56 #include "AliFlowCandidateTrack.h"
  57 #include "AliFlowTrackCuts.h"
  58 #include "AliFlowEventSimple.h"
  59 #include "AliFlowCommonConstants.h"
  60 #include "AliFlowEvent.h"
  61 #include "AliFlowTrack.h"
  62 #include "TProfile.h"
  63
  64 ClassImp(AliAnalysisTaskNucleiv2)
  65
  66 using std::cout;
  67 using std::endl;
  68
  69 //________________________________________________________________________
  70 AliAnalysisTaskNucleiv2::AliAnalysisTaskNucleiv2()
  71 : AliAnalysisTaskSE(),
  72   fAnalysisType(0),
  73   fCollidingSystems(0),
  74   fDataType(0),
  75   fFillNtuple(0),
  76   fCentralityMin(0),
  77   fCentralityMax(0),
  78   fCutsRP(0),
  79   fNullCuts(0),
  80   fFlowEvent(0),
  81   fListHist(0),
  82   fHistEventMultiplicity(0),
  83   fHistTrackMultiplicity(0),
  84   fHistTrackMultiplicityCentral(0),
  85   fHistTrackMultiplicitySemiCentral(0),
  86   fHistTrackMultiplicityMB(0),
  87   fhBB(0),
  88   fhBBDeu(0),
  89   fhPtDeu(0),
  90   fhTOF(0),
  91   fhMassTOF(0),
  92   EPVzAvsCentrality(0),
  93   EPVzCvsCentrality(0),
  94   EPTPCvsCentrality(0),
  95   EPVzvsCentrality(0),
  96   EPTPCpvsCentrality(0),
  97   EPTPCnvsCentrality(0),
  98   fSubEventDPhiv205(0),
  99   fSubEventDPhiv2new05(0),
 100   fSubEventDPhiv22040(0),
 101   fSubEventDPhiv2new2040(0),
 102   fSubEventDPhiv24060(0),
 103   fSubEventDPhiv2new4060(0),
 104   hCos2DeltaPhiVzAvsCentrality(0),
 105   hCos2DeltaPhiVzCvsCentrality(0),
 106   hCos2DeltaPhiVzMvsCentrality(0),
 107   hCos2DeltaPhiTPCfvsCentrality(0),
 108   hCos2DeltaPhiTPCpvsCentrality(0),
 109   hCos2DeltaPhiTPCnvsCentrality(0),
 110   hEvPlaneTPCvsEvPVz05(0),
 111   hEvPlaneTPCvsEvPVz075(0),
 112   hEvPlaneTPCvsEvPVz1530(0),
 113   hEvPlaneTPCvsEvPVz3050(0),
 114   hEvPlaneTPCvsEvPVz2040(0),
 115   hEvPlaneTPCvsEvPVz4060(0),
 116   hCos2DeltaPhivsPt075(0),
 117   hCos2DeltaPhiVZEROvsPt075(0),
 118   hCos2DeltaPhivsPt1530(0),
 119   hCos2DeltaPhiVZEROvsPt1530(0),
 120   hCos2DeltaPhivsPt3050(0),
 121   hCos2DeltaPhiVZEROvsPt3050(0),
 122   hCos2DeltaPhivsPt05(0),
 123   hCos2DeltaPhiVZEROvsPt05(0),
 124   hCos2DeltaPhivsPt2040(0),
 125   hCos2DeltaPhiVZEROvsPt2040(0),
 126   hCos2DeltaPhivsPt4060(0),
 127   hCos2DeltaPhiVZEROvsPt4060(0),
 128   fESDtrackCuts(0),
 129   fESDtrackCutsEP(0),
 130   fPIDResponse(0),
 131   fNtuple1(0),
 132   tCharge(0),
 133   tPhi(0),
 134   trpangleTPC(0),
 135   tPDGCode(0),
 136   tPDGCodeMum(0),
 137   tIsPrimaryTr(0),
 138   fNtuple2(0) ,
 139   tCentralityMC(0),
 140   tPDGCodeMC(0),
 141   tPDGCodeMumMC(0),
 142   tIsPrimary(0),
 143   tEtaMC(0),
 144   tPtMC(0),
 145   tYMC(0)
 146 {
 147   // Dummy Constructor
 148   fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
 149   fESDtrackCutsEP = new AliESDtrackCuts("AliESDtrackCutsEP","AliESDtrackCutsEP");
 150   //
 151   Initialize();
 152 }
 153
 154 //________________________________________________________________________
 155 AliAnalysisTaskNucleiv2::AliAnalysisTaskNucleiv2(const char *name)
 156   :  AliAnalysisTaskSE(name),
 157      fAnalysisType(0),
 158      fCollidingSystems(0),
 159      fDataType(0),
 160      fFillNtuple(0),
 161      fCentralityMin(0),
 162      fCentralityMax(0),
 163      fCutsRP(0),
 164      fNullCuts(0),
 165      fFlowEvent(0),
 166      fListHist(0),
 167      fHistEventMultiplicity(0),
 168      fHistTrackMultiplicity(0),
 169      fHistTrackMultiplicityCentral(0),
 170      fHistTrackMultiplicitySemiCentral(0),
 171      fHistTrackMultiplicityMB(0),
 172      fhBB(0),
 173      fhBBDeu(0),
 174      fhPtDeu(0),
 175      fhTOF(0),
 176      fhMassTOF(0),
 177      EPVzAvsCentrality(0),
 178      EPVzCvsCentrality(0),
 179      EPTPCvsCentrality(0),
 180      EPVzvsCentrality(0),
 181      EPTPCpvsCentrality(0),
 182      EPTPCnvsCentrality(0),
 183      fSubEventDPhiv205(0),
 184      fSubEventDPhiv2new05(0),
 185      fSubEventDPhiv22040(0),
 186      fSubEventDPhiv2new2040(0),
 187      fSubEventDPhiv24060(0),
 188      fSubEventDPhiv2new4060(0),
 189      hCos2DeltaPhiVzAvsCentrality(0),
 190      hCos2DeltaPhiVzCvsCentrality(0),
 191      hCos2DeltaPhiVzMvsCentrality(0),
 192      hCos2DeltaPhiTPCfvsCentrality(0),
 193      hCos2DeltaPhiTPCpvsCentrality(0),
 194      hCos2DeltaPhiTPCnvsCentrality(0),
 195      hEvPlaneTPCvsEvPVz05(0),
 196      hEvPlaneTPCvsEvPVz075(0),
 197      hEvPlaneTPCvsEvPVz1530(0),
 198      hEvPlaneTPCvsEvPVz3050(0),
 199      hEvPlaneTPCvsEvPVz2040(0),
 200      hEvPlaneTPCvsEvPVz4060(0),
 201      hCos2DeltaPhivsPt075(0),
 202      hCos2DeltaPhiVZEROvsPt075(0),
 203      hCos2DeltaPhivsPt1530(0),
 204      hCos2DeltaPhiVZEROvsPt1530(0),
 205      hCos2DeltaPhivsPt3050(0),
 206      hCos2DeltaPhiVZEROvsPt3050(0),
 207      hCos2DeltaPhivsPt05(0),
 208      hCos2DeltaPhiVZEROvsPt05(0),
 209      hCos2DeltaPhivsPt2040(0),
 210      hCos2DeltaPhiVZEROvsPt2040(0),
 211      hCos2DeltaPhivsPt4060(0),
 212      hCos2DeltaPhiVZEROvsPt4060(0),
 213      fESDtrackCuts(0),
 214      fESDtrackCutsEP(0),
 215      fPIDResponse(0),
 216      fNtuple1(0),
 217      tCharge(0),
 218      tPhi(0),
 219      trpangleTPC(0),
 220      tPDGCode(0),
 221      tPDGCodeMum(0),
 222      tIsPrimaryTr(0),
 223      fNtuple2(0) ,
 224      tCentralityMC(0),
 225      tPDGCodeMC(0),
 226      tPDGCodeMumMC(0),
 227      tIsPrimary(0),
 228      tEtaMC(0),
 229      tPtMC(0),
 230      tYMC(0)
 231
 232 {
 233   // Define input and output slots here
 234   // Input slot #0 works with a TChain
 235   //DefineInput(0, TChain::Class());
 236   // Output slot #0 writes into a TList container ()
 237
 238   //
 239   // create track cuts
 240   //
 241   fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
 242   //
 243   Initialize();
 244
 245   DefineInput(0, TChain::Class());
 246
 247   DefineOutput(1, TList::Class());
 248   DefineOutput(2, TTree::Class());
 249   DefineOutput(3, TTree::Class());
 250   DefineOutput(4, AliFlowEventSimple::Class());
 251 }
 252
 253 void AliAnalysisTaskNucleiv2::Initialize()
 254 {
 255   //
 256   // updating parameters in case of changes
 257   //
 258   fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,kTRUE);
 259   fESDtrackCuts->SetMaxDCAToVertexXY(3);
 260   fESDtrackCuts->SetMaxDCAToVertexZ(2);
 261   fESDtrackCuts->SetEtaRange(-0.8,0.8);
 262   //
 263   //
 264
 265   fESDtrackCutsEP = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 266   //  Printf("Initizialize\n");
 267
 268 }
 269
 270 //________________________________________________________________________
 271 Float_t AliAnalysisTaskNucleiv2::GetEventPlaneForCandidate(AliESDtrack* track0, const TVector2* q,AliEventplane *pl, const TVector2* qsub1, const TVector2* qsub2){
 272
 273   // remove autocorrelations
 274
 275   TArrayF* qx = 0x0;
 276   TArrayF* qy = 0x0;
 277   TVector2 qcopy;
 278   // if(!fEtaGap){
 279   qx = pl->GetQContributionXArray();
 280   qy = pl->GetQContributionYArray();
 281   qcopy = *q;
 282   // }else {
 283   // if(d->Eta()<0.){
 284   //   qx = pl->GetQContributionXArraysub1();
 285   //   qy = pl->GetQContributionYArraysub1();
 286   //   qcopy = *qsub1;
 287   // }else{
 288   //   qx = pl->GetQContributionXArraysub2();
 289   //   qy = pl->GetQContributionYArraysub2();
 290   //   qcopy = *qsub2;
 291   // }
 292   //}
 293
 294   TVector2 q0;
 295   if((track0->GetID()) < qx->fN){
 296     q0.Set(qx->At(track0->GetID()),qy->At(track0->GetID()));
 297   }
 298
 299   qcopy = qcopy - q0;
 300
 301   return qcopy.Phi()/2.;
 302
 303 }
 304 //_____________________________________________________________________________
 305 template <typename T> void AliAnalysisTaskNucleiv2::SetNullCuts(T* event)
 306 {
 307     //Set null cuts
 308   if (fDebug) cout << " fCutsRP " << fCutsRP << endl;
 309   fCutsRP->SetEvent(event, MCEvent());
 310   fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
 311   fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
 312   fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
 313   fNullCuts->SetEvent(event, MCEvent());
 314 }
 315 //______________________________________________________________________
 316 void AliAnalysisTaskNucleiv2::PrepareFlowEvent(Int_t iMulti, AliFlowEvent *FlowEv) const
 317 {
 318     //Prepare flow events
 319     FlowEv->ClearFast();
 320     FlowEv->Fill(fCutsRP, fNullCuts);
 321     FlowEv->SetReferenceMultiplicity(iMulti);
 322     FlowEv->DefineDeadZone(0, 0, 0, 0);
 323     //  FlowEv->TagSubeventsInEta(-0.7, 0, 0, 0.7);
 324 }
 325
 326 //________________________________________________________________________
 327 Float_t AliAnalysisTaskNucleiv2::GetPhi0Pi(Float_t phi){
 328   // Sets the phi angle in the range 0-pi
 329   Float_t result=phi;
 330   while(result<0){
 331     result=result+TMath::Pi();
 332   }
 333   while(result>TMath::Pi()){
 334     result=result-TMath::Pi();
 335   }
 336    return result;
 337 }
 338
 339 //_____________________________________________________________________________
 340 void AliAnalysisTaskNucleiv2::SetCentralityParameters(Double_t CentralityMin, Double_t CentralityMax)
 341 {
 342     // Set a centrality range ]min, max]
 343   fCentralityMin = CentralityMin;
 344   fCentralityMax = CentralityMax;
 345
 346 }
 347 //==================DEFINITION OF OUTPUT OBJECTS==============================
 348
 349 void AliAnalysisTaskNucleiv2::UserCreateOutputObjects()
 350 {
 351   //-------------------------------------------------------
 352   fNullCuts = new AliFlowTrackCuts("null_cuts");
 353
 354   AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
 355   cc->SetNbinsMult(10000);
 356   cc->SetMultMin(0);
 357   cc->SetMultMax(10000);
 358
 359   cc->SetNbinsPt(100);
 360   cc->SetPtMin(0);
 361   cc->SetPtMax(50);
 362
 363   cc->SetNbinsPhi(180);
 364   cc->SetPhiMin(0.0);
 365   cc->SetPhiMax(TMath::TwoPi());
 366
 367   cc->SetNbinsEta(30);
 368   cc->SetEtaMin(-7.0);
 369   cc->SetEtaMax(+7.0);
 370
 371   cc->SetNbinsQ(500);
 372   cc->SetQMin(0.0);
 373   cc->SetQMax(3.0);
 374
 375   //------------------------------------------------------
 376
 377   fListHist = new TList();
 378   fListHist->SetOwner();  // IMPORTANT!
 379
 380   if(! fHistEventMultiplicity ){
 381
 382     fHistEventMultiplicity   = new TH1F( "fHistEventMultiplicity" , "Nb of Events" , 12 , -0.5,11.5);
 383
 384     fHistEventMultiplicity->GetXaxis()->SetBinLabel(1,"All Events");
 385     fHistEventMultiplicity->GetXaxis()->SetBinLabel(2,"Events w/PV");
 386     fHistEventMultiplicity->GetXaxis()->SetBinLabel(3,"Events w/|Vz|<10cm");
 387     fHistEventMultiplicity->GetXaxis()->SetBinLabel(4,"Central Events");
 388     fHistEventMultiplicity->GetXaxis()->SetBinLabel(5,"Semi-Central Events");
 389     fHistEventMultiplicity->GetXaxis()->SetBinLabel(6,"MB Events");
 390     //from HF
 391     fHistEventMultiplicity->GetXaxis()->SetBinLabel(7,"nEventsAnal");
 392     fHistEventMultiplicity->GetXaxis()->SetBinLabel(8,"nEvSelected");
 393     fHistEventMultiplicity->GetXaxis()->SetBinLabel(9,"nCandidatesSelected");
 394     fHistEventMultiplicity->GetXaxis()->SetBinLabel(10,"out of pt bounds");
 395     fHistEventMultiplicity->GetXaxis()->SetBinLabel(11,"mismatch lab");
 396     fHistEventMultiplicity->GetXaxis()->SetBinLabel(12,"non valid TPC EP");
 397     fListHist->Add(fHistEventMultiplicity);
 398   }
 399
 400   if(! fHistTrackMultiplicity ){
 401     fHistTrackMultiplicity  = new TH2F( "fHistTrackMultiplicity", "Nb of Tracks MB Events |Vz| < 10", 25000,0, 25000,105,-0.5,104.5);
 402     fHistTrackMultiplicity->GetXaxis()->SetTitle("Number of tracks");
 403     fHistTrackMultiplicity->GetYaxis()->SetTitle("Percentile");
 404     fListHist->Add(fHistTrackMultiplicity);
 405   }
 406
 407   if(! fHistTrackMultiplicityCentral ){
 408     fHistTrackMultiplicityCentral  = new TH2F( "fHistTrackMultiplicityCentral", "Nb of Tracks MB Events |Vz| < 10", 25000,0, 25000,105,-0.5,104.5);
 409     fHistTrackMultiplicityCentral->GetXaxis()->SetTitle("Number of tracks");
 410     fHistTrackMultiplicityCentral->GetYaxis()->SetTitle("Percentile");
 411     fListHist->Add(fHistTrackMultiplicityCentral);
 412   }
 413   if(! fHistTrackMultiplicitySemiCentral ){
 414     fHistTrackMultiplicitySemiCentral  = new TH2F( "fHistTrackMultiplicitySemiCentral", "Nb of Tracks MB Events |Vz| < 10", 25000,0, 25000,105,-0.5,104.5);
 415     fHistTrackMultiplicitySemiCentral->GetXaxis()->SetTitle("Number of tracks");
 416     fHistTrackMultiplicitySemiCentral->GetYaxis()->SetTitle("Percentile");
 417     fListHist->Add(fHistTrackMultiplicitySemiCentral);
 418   }
 419   if(! fHistTrackMultiplicityMB ){
 420     fHistTrackMultiplicityMB  = new TH2F( "fHistTrackMultiplicityMB", "Nb of Tracks MB Events |Vz| < 10", 25000,0, 25000,105,-0.5,104.5);
 421     fHistTrackMultiplicityMB->GetXaxis()->SetTitle("Number of tracks");
 422     fHistTrackMultiplicityMB->GetYaxis()->SetTitle("Percentile");
 423     fListHist->Add(fHistTrackMultiplicityMB);
 424   }
 425
 426   if(! fhBB ){
 427     fhBB = new TH2F( "fhBB" , "BetheBlochTPC" , 240,-6,6,1000,0,1000);
 428     fListHist->Add(fhBB);
 429   }
 430
 431   if(! fhBBDeu ){
 432     fhBBDeu = new TH2F( "fhBBDeu" , "BetheBlochTPC - Deuteron" , 240,-6,6,1000,0,1000);
 433     fListHist->Add(fhBBDeu);
 434   }
 435
 436   if(!fhPtDeu  ){
 437     fhPtDeu = new TH2F( "fhPtDeu" , "pt corretto vs pt track - Deuteron" , 120,0,6,120,0,6);
 438     fListHist->Add(fhPtDeu);
 439   }
 440
 441   if(! fhTOF ){
 442     fhTOF = new TH2F( "fhTOF" , "Scatter Plot TOF" , 240,-6,6,1000,0,1.2);
 443     fListHist->Add(fhTOF);
 444   }
 445   if(! fhMassTOF){
 446     fhMassTOF=new TH1F ("fhMassTOF","Particle Mass - TOF", 300,-5 ,5);
 447     fListHist->Add(fhMassTOF);
 448   }
 449
 450   EPVzAvsCentrality = new TH2D("EPVzAvsCentrality", "EPVzAvsCentrality", 80, -2, 2,105,-0.5,105.5);
 451   fListHist->Add(EPVzAvsCentrality);
 452   EPVzCvsCentrality = new TH2D("EPVzCvsCentrality", "EPVzCvsCentrality", 80, -2, 2,105,-0.5,105.5);
 453   fListHist->Add(EPVzCvsCentrality);
 454   EPTPCvsCentrality = new TH2D("EPTPCvsCentrality", "EPTPCvsCentrality", 80, -2, 2,105,-0.5,105.5);
 455   fListHist->Add(EPTPCvsCentrality);
 456
 457
 458   EPVzvsCentrality = new TH2D("EPVzvsCentrality", "EPVzvsCentrality", 80, -2, 2,105,-0.5,105.5);
 459   fListHist->Add(EPVzvsCentrality);
 460   EPTPCpvsCentrality = new TH2D("EPTPCpvsCentrality", "EPTPCpvsCentrality", 80, -2, 2,105,-0.5,105.5);
 461   fListHist->Add(EPTPCpvsCentrality);
 462   EPTPCnvsCentrality = new TH2D("EPTPCnvsCentrality", "EPTPCnvsCentrality", 80, -2, 2,105,-0.5,105.5);
 463   fListHist->Add(EPTPCnvsCentrality);
 464
 465   //1
 466
 467   fSubEventDPhiv205 = new TProfile("fSubEventDPhiv205", "fSubEventDPhiv205", 3, 0, 3);
 468   fSubEventDPhiv205->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
 469   fSubEventDPhiv205->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
 470   fSubEventDPhiv205->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
 471   fListHist->Add(fSubEventDPhiv205);
 472
 473   fSubEventDPhiv2new05 = new TProfile("fSubEventDPhiv2new05", "fSubEventDPhiv2new05", 3, 0, 3);
 474   fSubEventDPhiv2new05->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
 475   fSubEventDPhiv2new05->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
 476   fSubEventDPhiv2new05->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
 477   fListHist->Add(fSubEventDPhiv2new05);
 478
 479   //2
 480
 481   fSubEventDPhiv22040 = new TProfile("fSubEventDPhiv22040", "fSubEventDPhiv22040", 3, 0, 3);
 482   fSubEventDPhiv22040->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
 483   fSubEventDPhiv22040->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
 484   fSubEventDPhiv22040->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
 485   fListHist->Add(fSubEventDPhiv22040);
 486
 487   fSubEventDPhiv2new2040 = new TProfile("fSubEventDPhiv2new2040", "fSubEventDPhiv2new2040", 3, 0, 3);
 488   fSubEventDPhiv2new2040->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
 489   fSubEventDPhiv2new2040->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
 490   fSubEventDPhiv2new2040->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
 491   fListHist->Add(fSubEventDPhiv2new2040);
 492
 493   //3
 494
 495   fSubEventDPhiv24060 = new TProfile("fSubEventDPhiv24060", "fSubEventDPhiv24060", 3, 0, 3);
 496   fSubEventDPhiv24060->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
 497   fSubEventDPhiv24060->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
 498   fSubEventDPhiv24060->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
 499   fListHist->Add(fSubEventDPhiv24060);
 500
 501   fSubEventDPhiv2new4060 = new TProfile("fSubEventDPhiv2new4060", "fSubEventDPhiv2new4060", 3, 0, 3);
 502   fSubEventDPhiv2new4060->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
 503   fSubEventDPhiv2new4060->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
 504   fSubEventDPhiv2new4060->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
 505   fListHist->Add(fSubEventDPhiv2new4060);
 506
 507   //------------------
 508
 509   hCos2DeltaPhiVzAvsCentrality  = new TH2F("hCos2DeltaPhiVzAvsCentrality" ,"hCos2DeltaPhiVzAvsCentrality" ,100,-1.1,1.1,105,-0.5,105.5);
 510   hCos2DeltaPhiVzCvsCentrality  = new TH2F("hCos2DeltaPhiVzCvsCentrality" ,"hCos2DeltaPhiVzCvsCentrality" ,100,-1.1,1.1,105,-0.5,105.5);
 511   hCos2DeltaPhiVzMvsCentrality  = new TH2F("hCos2DeltaPhiVzMvsCentrality" ,"hCos2DeltaPhiVzMvsCentrality" ,100,-1.1,1.1,105,-0.5,105.5);
 512   hCos2DeltaPhiTPCfvsCentrality = new TH2F("hCos2DeltaPhiTPCfvsCentrality","hCos2DeltaPhiTPCfvsCentrality",100,-1.1,1.1,105,-0.5,105.5);
 513   hCos2DeltaPhiTPCpvsCentrality = new TH2F("hCos2DeltaPhiTPCpvsCentrality","hCos2DeltaPhiTPCpvsCentrality",100,-1.1,1.1,105,-0.5,105.5);
 514   hCos2DeltaPhiTPCnvsCentrality = new TH2F("hCos2DeltaPhiTPCnvsCentrality","hCos2DeltaPhiTPCnvsCentrality",100,-1.1,1.1,105,-0.5,105.5);
 515
 516   fListHist->Add(hCos2DeltaPhiVzAvsCentrality);
 517   fListHist->Add(hCos2DeltaPhiVzCvsCentrality);
 518   fListHist->Add(hCos2DeltaPhiVzMvsCentrality);
 519   fListHist->Add(hCos2DeltaPhiTPCfvsCentrality);
 520   fListHist->Add(hCos2DeltaPhiTPCpvsCentrality);
 521   fListHist->Add(hCos2DeltaPhiTPCnvsCentrality);
 522
 523   hEvPlaneTPCvsEvPVz05   = new TH2F("hEvPlaneTPCvsEvPVz05"  ,"hEvPlaneTPCvsEvPVz05"  ,100,-2*TMath::Pi(),2*TMath::Pi(),100,-2*TMath::Pi(),2*TMath::Pi());
 524   hEvPlaneTPCvsEvPVz075  = new TH2F("hEvPlaneTPCvsEvPVz075" ,"hEvPlaneTPCvsEvPVz075" ,100,-2*TMath::Pi(),2*TMath::Pi(),100,-2*TMath::Pi(),2*TMath::Pi());
 525   hEvPlaneTPCvsEvPVz1530 = new TH2F("hEvPlaneTPCvsEvPVz1530","hEvPlaneTPCvsEvPVz1530",100,-2*TMath::Pi(),2*TMath::Pi(),100,-2*TMath::Pi(),2*TMath::Pi());
 526   hEvPlaneTPCvsEvPVz3050 = new TH2F("hEvPlaneTPCvsEvPVz3050","hEvPlaneTPCvsEvPVz3050",100,-2*TMath::Pi(),2*TMath::Pi(),100,-2*TMath::Pi(),2*TMath::Pi());
 527   hEvPlaneTPCvsEvPVz2040 = new TH2F("hEvPlaneTPCvsEvPVz2040","hEvPlaneTPCvsEvPVz2040",100,-2*TMath::Pi(),2*TMath::Pi(),100,-2*TMath::Pi(),2*TMath::Pi());
 528   hEvPlaneTPCvsEvPVz4060 = new TH2F("hEvPlaneTPCvsEvPVz4060","hEvPlaneTPCvsEvPVz4060",100,-2*TMath::Pi(),2*TMath::Pi(),100,-2*TMath::Pi(),2*TMath::Pi());
 529
 530
 531   fListHist->Add(hEvPlaneTPCvsEvPVz05);
 532   fListHist->Add(hEvPlaneTPCvsEvPVz075);
 533   fListHist->Add(hEvPlaneTPCvsEvPVz1530);
 534   fListHist->Add(hEvPlaneTPCvsEvPVz3050);
 535   fListHist->Add(hEvPlaneTPCvsEvPVz2040);
 536   fListHist->Add(hEvPlaneTPCvsEvPVz4060);
 537
 538   hCos2DeltaPhivsPt075          = new TH2F("hCos2DeltaPhivsPt075",      "hCos2DeltaPhivsPt075"       ,75, -1.1,1.1,210,-5.5,5.5);
 539   hCos2DeltaPhiVZEROvsPt075     = new TH2F("hCos2DeltaPhiVZEROvsPt075", "hCos2DeltaPhiVZEROvsPt075"  ,75, -1.1,1.1,210,-5.5,5.5);
 540   hCos2DeltaPhivsPt1530         = new TH2F("hCos2DeltaPhivsPt1530",     "hCos2DeltaPhivsPt1530"      ,75, -1.1,1.1,210,-5.5,5.5);
 541   hCos2DeltaPhiVZEROvsPt1530    = new TH2F("hCos2DeltaPhiVZEROvsPt1530","hCos2DeltaPhiVZEROvsPt1530" ,75, -1.1,1.1,210,-5.5,5.5);
 542   hCos2DeltaPhivsPt3050         = new TH2F("hCos2DeltaPhivsPt3050",     "hCos2DeltaPhivsPt3050"      ,75, -1.1,1.1,210,-5.5,5.5);
 543   hCos2DeltaPhiVZEROvsPt3050    = new TH2F("hCos2DeltaPhiVZEROvsPt3050","hCos2DeltaPhiVZEROvsPt3050" ,75, -1.1,1.1,210,-5.5,5.5);
 544   hCos2DeltaPhivsPt05           = new TH2F("hCos2DeltaPhivsPt05",       "hCos2DeltaPhivsPt05"        ,75, -1.1,1.1,210,-5.5,5.5);
 545   hCos2DeltaPhiVZEROvsPt05      = new TH2F("hCos2DeltaPhiVZEROvsPt05",  "hCos2DeltaPhiVZEROvsPt05"   ,75, -1.1,1.1,210,-5.5,5.5);
 546   hCos2DeltaPhivsPt2040         = new TH2F("hCos2DeltaPhivsPt2040",     "hCos2DeltaPhivsPt2040"      ,75, -1.1,1.1,210,-5.5,5.5);
 547   hCos2DeltaPhiVZEROvsPt2040    = new TH2F("hCos2DeltaPhiVZEROvsPt2040","hCos2DeltaPhiVZEROvsPt2040" ,75, -1.1,1.1,210,-5.5,5.5);
 548   hCos2DeltaPhivsPt4060         = new TH2F("hCos2DeltaPhivsPt4060",     "hCos2DeltaPhivsPt4060"      ,75, -1.1,1.1,210,-5.5,5.5);
 549   hCos2DeltaPhiVZEROvsPt4060    = new TH2F("hCos2DeltaPhiVZEROvsPt4060","hCos2DeltaPhiVZEROvsPt4060" ,75, -1.1,1.1,210,-5.5,5.5);
 550
 551   //--------------
 552   fListHist->Add(hCos2DeltaPhivsPt075);
 553   fListHist->Add(hCos2DeltaPhiVZEROvsPt075);
 554   fListHist->Add(hCos2DeltaPhivsPt1530);
 555   fListHist->Add(hCos2DeltaPhiVZEROvsPt1530);
 556   fListHist->Add(hCos2DeltaPhivsPt3050);
 557   fListHist->Add(hCos2DeltaPhiVZEROvsPt3050);
 558   fListHist->Add(hCos2DeltaPhivsPt05);
 559   fListHist->Add(hCos2DeltaPhiVZEROvsPt05);
 560   fListHist->Add(hCos2DeltaPhivsPt2040);
 561   fListHist->Add(hCos2DeltaPhiVZEROvsPt2040);
 562   fListHist->Add(hCos2DeltaPhivsPt4060);
 563   fListHist->Add(hCos2DeltaPhiVZEROvsPt4060);
 564
 565
 566   if(! fNtuple1 ) {
 567
 568     fNtuple1 = new TTree("fNtuple1","fNtuple1");
 569
 570     fNtuple1->Branch("tCentrality"   ,&tCentrality   ,"tCentrality[2]/D"    );
 571     fNtuple1->Branch("tPulls"        ,&tPulls        ,"tPulls[3]/D"   );
 572     fNtuple1->Branch("tMomentum"     ,&tMomentum     ,"tMomentum[4]/D"   );
 573     fNtuple1->Branch("tDCA"          ,&tDCA          ,"tDCA[2]/D"        );
 574     fNtuple1->Branch("tisTOF"        ,&tisTOF        ,"tisTOF[2]/I"      );
 575     fNtuple1->Branch("tTOFtrack"     ,&tTOFtrack     ,"tTOFtrack[3]/D"   );
 576     fNtuple1->Branch("tCharge"       ,&tCharge       ,"tCharge/I"        );
 577     fNtuple1->Branch("tPhi"          ,&tPhi          ,"tPhi/D"           );
 578     fNtuple1->Branch("trpangleTPC"   ,&trpangleTPC   ,"trpangleTPC/D"    );
 579     fNtuple1->Branch("trpangleVZERO" ,&trpangleVZERO ,"trpangleVZERO[3]/D"  );
 580
 581     if(fDataType == "SIM"){
 582       fNtuple1->Branch("tPDGCode"      ,&tPDGCode      ,"tPDGCode/D"       );
 583       fNtuple1->Branch("tPDGCodeMum"   ,&tPDGCode      ,"tPDGCodeMum/D"    );
 584       fNtuple1->Branch("tIsPrimaryTr"  ,&tIsPrimaryTr  ,"tIsPrimaryTr/D"   );
 585       fNtuple1->Branch("tIsSecondaryTr",&tIsSecondaryTr,"tIsSecondaryTr[2]/D" );
 586     }
 587
 588   }
 589
 590
 591   if(fDataType == "SIM"){
 592
 593     if(! fNtuple2 ) {
 594
 595       fNtuple2 = new TTree("fNtuple2","fNtuple2");
 596
 597       fNtuple2->Branch("tCentralityMC"   ,&tCentralityMC   ,"tCentralityMC/D"    );
 598       fNtuple2->Branch("tVertexCoordMC"  ,&tVertexCoordMC  ,"tVertexCoordMC[3]/D");
 599       fNtuple2->Branch("tMomentumMC"     ,&tMomentumMC     ,"tMomentumMC[3]/D"   );
 600       fNtuple2->Branch("tPDGCodeMC"      ,&tPDGCodeMC      ,"tPDGCodeMC/D"       );
 601       fNtuple2->Branch("tPDGCodeMumMC"   ,&tPDGCodeMC      ,"tPDGCodeMumMC/D"    );
 602       fNtuple2->Branch("tIsPrimary"      ,&tIsPrimary      ,"tIsPrimary/D"       );
 603       fNtuple2->Branch("tIsSecondary"    ,&tIsSecondary    ,"tIsSecondary[2]/D"     );
 604       fNtuple2->Branch("tEtaMC"          ,&tEtaMC          ,"tEtaMC/D"         );
 605       fNtuple2->Branch("tPtMC"           ,&tPtMC           ,"tPtMC/D"          );
 606       fNtuple2->Branch("tYMC"            ,&tYMC            ,"tYMC/D"           );
 607
 608     }
 609
 610   }
 611
 612   PostData(1,  fListHist);
 613   PostData(2,  fNtuple1);
 614   PostData(3,  fNtuple2);
 615
 616   fFlowEvent = new AliFlowEvent(10000);
 617   PostData(4, fFlowEvent);
 618
 619 }// end UserCreateOutputObjects
 620
 621
 622 //====================== USER EXEC ========================
 623
 624 void AliAnalysisTaskNucleiv2::UserExec(Option_t *)
 625 {
 626   // Main loop
 627   // Called for EACH event
 628   //  cout<<"AliAnalysisTaskNucleiv2 Starting UserExec"<<endl;
 629
 630   Info("AliAnalysisTaskNucleiv2","Starting UserExec");
 631
 632   AliVEvent *event = InputEvent();
 633   if (!event) { Printf("ERROR: Could not retrieve event"); return; }
 634
 635   AliStack *stack = 0;
 636
 637   if(fDataType == "SIM"){
 638     AliMCEvent *mcEvent = MCEvent();
 639     if (!mcEvent) {
 640       Printf("ERROR: Could not retrieve MC event");
 641       return;
 642     }
 643     Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
 644     stack = mcEvent->Stack();
 645     if( !stack ) {
 646       Printf( "Stack not available");
 647       return;
 648     }
 649   }
 650
 651   // create pointer to event
 652
 653   if(fAnalysisType == "ESD"){
 654
 655     AliESDEvent* lESDevent = dynamic_cast<AliESDEvent*>(event);
 656     if (!lESDevent) {
 657       AliError("Cannot get the ESD event");
 658       return;
 659     }
 660
 661     fHistEventMultiplicity->Fill(1);
 662     fHistEventMultiplicity->Fill(7);
 663
 664     //_____________________________________________________
 665     //   Centrality
 666
 667     AliCentrality *centrality = lESDevent->GetCentrality();
 668     Float_t percentile=centrality->GetCentralityPercentile("V0M");
 669
 670     Int_t TrackNumber = lESDevent->GetNumberOfTracks();
 671     fHistTrackMultiplicity->Fill(TrackNumber,percentile); //tracce per evento
 672
 673     //______________________________________________________
 674     // PID
 675
 676     AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 677     AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 678     fPIDResponse=inputHandler->GetPIDResponse();
 679
 680     //=================================================================
 681
 682     Float_t impactXY=-999., impactZ=-999.;
 683     Double_t pinTPC=0.,poutTPC=0.,TPCSignal=0.;
 684
 685     ULong_t  status=0;
 686     Bool_t   isTPC=kFALSE;
 687
 688
 689     // Primary vertex cut
 690
 691     const AliESDVertex *vtx = lESDevent->GetPrimaryVertexTracks();
 692
 693     if(vtx->GetNContributors()<1) {
 694
 695       // SPD vertex cut
 696       vtx = lESDevent->GetPrimaryVertexSPD();
 697
 698       if(vtx->GetNContributors()<1) {
 699         Info("AliAnalysisTaskHelium3Pi","No good vertex, skip event");
 700         return; // NO GOOD VERTEX, SKIP EVENT
 701       }
 702     }
 703
 704     fHistEventMultiplicity->Fill(2); // analyzed events with PV
 705
 706     if(TMath::Abs(vtx->GetZv())>10) return;
 707     fHistEventMultiplicity->Fill(3);
 708
 709     Bool_t isSelectedCentral     = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
 710     Bool_t isSelectedSemiCentral = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
 711     Bool_t isSelectedMB          = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 712
 713     fHistTrackMultiplicity->Fill(TrackNumber,percentile);
 714
 715     Int_t eventtype = -999;
 716
 717     //  cout<<"ET 1: "<<eventtype<<endl;
 718
 719     if(fDataType == "REAL"){
 720
 721       if(isSelectedCentral){
 722         // if(percentile<0)return;
 723         // if(percentile>=7.5)return;
 724         fHistEventMultiplicity->Fill(4);
 725         fHistTrackMultiplicityCentral->Fill(TrackNumber,percentile);
 726         eventtype =1;
 727       }
 728
 729       if(isSelectedSemiCentral){
 730         // if(percentile<15)return;
 731         // if(percentile>=50)return;
 732         fHistEventMultiplicity->Fill(5);
 733         fHistTrackMultiplicitySemiCentral->Fill(TrackNumber,percentile);
 734         eventtype =2;
 735       }
 736
 737       if(isSelectedMB){
 738         if(percentile<0)return;
 739         if(percentile>=80)return;
 740         fHistEventMultiplicity->Fill(6);
 741         fHistTrackMultiplicityMB->Fill(TrackNumber,percentile);
 742         eventtype =3;
 743       }
 744
 745       //    cout<<"ET 2: "<<eventtype<<endl;
 746
 747       if(eventtype!=1 && eventtype!=2 && eventtype!=3 )return;
 748     }
 749
 750     if(fDataType == "SIM"){
 751       cout<<"Take SIM event"<<endl;
 752       eventtype = -999;
 753       if(eventtype!=-999)return;
 754     }
 755     //----------------
 756     SetNullCuts(lESDevent);
 757     PrepareFlowEvent(TrackNumber,fFlowEvent);    //Calculate event plane Qvector and EP resolution for inclusive
 758
 759     AliEventplane *pl=lESDevent->GetEventplane();
 760
 761     if(fDataType == "REAL"){
 762
 763       if(!pl ){
 764         AliError("AliAnalysisTaskSENucleiv2::UserExec:no eventplane! v2 analysis without eventplane not possible!\n");
 765         fHistEventMultiplicity->Fill(12);
 766       }
 767     }
 768
 769
 770     //Event plane from FLOW
 771     //=============================================V0EP from Alex======================================================================
 772
 773     Double_t qxEPa = 0, qyEPa = 0;
 774     Double_t qxEPc = 0, qyEPc = 0;
 775     Double_t qxEP = 0, qyEP = 0;
 776
 777     Double_t evPlAngV0A = pl->CalculateVZEROEventPlane(lESDevent, 8, 2, qxEPa, qyEPa);
 778     Double_t evPlAngV0C = pl->CalculateVZEROEventPlane(lESDevent, 9, 2, qxEPc, qyEPc);
 779     Double_t evPlAngV0  = pl->CalculateVZEROEventPlane(lESDevent,10, 2, qxEP,  qyEP);
 780
 781     Double_t Qx2  = 0, Qy2  = 0;
 782     Double_t Qx2p = 0, Qy2p = 0;
 783     Double_t Qx2n = 0, Qy2n = 0;
 784
 785     for (Int_t iT = 0; iT < TrackNumber; iT++){
 786       AliESDtrack* track = lESDevent->GetTrack(iT);
 787
 788       if (!track)
 789         continue;
 790
 791       if ((TMath::Abs(track->Eta()) > 0.8) || (track->Pt() < 0.2) || (track->GetTPCNcls() < 70) || (track->Pt() >= 20.0))
 792         continue;
 793
 794       if(!fESDtrackCutsEP->AcceptTrack(track))
 795         continue;
 796
 797       if(track->Eta()>0 && track->Eta()<0.8){
 798
 799         Qx2p += TMath::Cos(2*track->Phi());
 800         Qy2p += TMath::Sin(2*track->Phi());
 801       }
 802       if(track->Eta()<0 && track->Eta()> -0.8){
 803
 804         Qx2n += TMath::Cos(2*track->Phi());
 805         Qy2n += TMath::Sin(2*track->Phi());
 806       }
 807
 808       Qx2 += TMath::Cos(2*track->Phi());
 809       Qy2 += TMath::Sin(2*track->Phi());
 810
 811     }
 812
 813     Double_t evPlAngTPC  = TMath::ATan2(Qy2, Qx2)/2.;
 814     Double_t evPlAngTPCn = TMath::ATan2(Qy2n, Qx2n)/2.;
 815     Double_t evPlAngTPCp = TMath::ATan2(Qy2p, Qx2p)/2.;
 816
 817     EPVzAvsCentrality->Fill(evPlAngV0A,percentile);
 818     EPVzCvsCentrality->Fill(evPlAngV0C,percentile);
 819     EPTPCvsCentrality->Fill(evPlAngTPC,percentile);
 820
 821     EPTPCnvsCentrality->Fill(evPlAngTPCn,percentile);
 822     EPTPCpvsCentrality->Fill(evPlAngTPCp,percentile);
 823     EPVzvsCentrality->Fill(evPlAngV0,percentile);
 824
 825     // For TPC resolution
 826     // Cos(2*(tpc-v0A)) vs Centrality
 827     // Cos(2*(tpc-v0C)) vs Centrality
 828     // Cos(2*(v0A-v0C)) vs Centrality
 829     // For VZEROM resolution :
 830     // Cos(2*(V0-TPCp)) vs Centrality
 831     // Cos(2*(V0-TPCn)) vs Centrality
 832     // Cos(2*(TPCp-TPCn)) vs Centrality
 833
 834     hCos2DeltaPhiVzAvsCentrality ->Fill(TMath::Cos(2.*(evPlAngV0-evPlAngTPCp))  ,percentile);
 835     hCos2DeltaPhiVzCvsCentrality ->Fill(TMath::Cos(2.*(evPlAngV0-evPlAngTPCn))  ,percentile);
 836     hCos2DeltaPhiVzMvsCentrality ->Fill(TMath::Cos(2.*(evPlAngTPCp-evPlAngTPCn)),percentile);
 837     hCos2DeltaPhiTPCfvsCentrality->Fill(TMath::Cos(2.*(evPlAngTPC-evPlAngV0A))  ,percentile);
 838     hCos2DeltaPhiTPCpvsCentrality->Fill(TMath::Cos(2.*(evPlAngTPC-evPlAngV0C))  ,percentile);
 839     hCos2DeltaPhiTPCnvsCentrality->Fill(TMath::Cos(2.*(evPlAngV0C-evPlAngV0A))  ,percentile);
 840
 841
 842     // e volendo 3 per il vzero come sotto
 843
 844     if(percentile>=0 || percentile<=5){
 845       //This is v0 resolution
 846       fSubEventDPhiv205->Fill(0.5, TMath::Cos(2.*(evPlAngV0A-evPlAngTPC))); // vzeroa - tpc
 847       fSubEventDPhiv205->Fill(1.5, TMath::Cos(2.*(evPlAngV0A-evPlAngV0C))); // vzeroa - vzeroc
 848       fSubEventDPhiv205->Fill(2.5, TMath::Cos(2.*(evPlAngV0C-evPlAngTPC))); // tpc - vzeroc
 849
 850
 851       fSubEventDPhiv2new05->Fill(0.5, TMath::Cos(2.*(evPlAngV0-evPlAngTPCp))  ); // vzero - tpcp
 852       fSubEventDPhiv2new05->Fill(1.5, TMath::Cos(2.*(evPlAngV0-evPlAngTPCn))  ); // vzero - tpcn
 853       fSubEventDPhiv2new05->Fill(2.5, TMath::Cos(2.*(evPlAngTPCp-evPlAngTPCn))); // tpcp - tpcn
 854
 855       hEvPlaneTPCvsEvPVz05  ->Fill(  evPlAngTPC,evPlAngV0);
 856     }
 857
 858     if(percentile>=20 || percentile<=40){
 859       fSubEventDPhiv22040->Fill(0.5, TMath::Cos(2.*(evPlAngV0A-evPlAngTPC))); // vzeroa - tpc
 860       fSubEventDPhiv22040->Fill(1.5, TMath::Cos(2.*(evPlAngV0A-evPlAngV0C))); // vzeroa - vzeroc
 861       fSubEventDPhiv22040->Fill(2.5, TMath::Cos(2.*(evPlAngV0C-evPlAngTPC))); // tpc - vzeroc
 862
 863
 864       fSubEventDPhiv2new2040->Fill(0.5, TMath::Cos(2.*(evPlAngV0-evPlAngTPCp))); // vzero - tpcp
 865       fSubEventDPhiv2new2040->Fill(1.5, TMath::Cos(2.*(evPlAngV0-evPlAngTPCn))); // vzero - tpcn
 866       fSubEventDPhiv2new2040->Fill(2.5, TMath::Cos(2.*(evPlAngTPCp-evPlAngTPCn))); // tpcp - tpcn
 867
 868       hEvPlaneTPCvsEvPVz2040   ->Fill(  evPlAngTPC,evPlAngV0);
 869
 870     }
 871
 872     if(percentile>=40 || percentile<=60){
 873       fSubEventDPhiv24060->Fill(0.5, TMath::Cos(2.*(evPlAngV0A-evPlAngTPC))); // vzeroa - tpc
 874       fSubEventDPhiv24060->Fill(1.5, TMath::Cos(2.*(evPlAngV0A-evPlAngV0C))); // vzeroa - vzeroc
 875       fSubEventDPhiv24060->Fill(2.5, TMath::Cos(2.*(evPlAngV0C-evPlAngTPC))); // tpc - vzeroc
 876
 877       fSubEventDPhiv2new4060->Fill(0.5, TMath::Cos(2.*(evPlAngV0-evPlAngTPCp))); // vzero - tpcp
 878       fSubEventDPhiv2new4060->Fill(1.5, TMath::Cos(2.*(evPlAngV0-evPlAngTPCn))); // vzero - tpcn
 879       fSubEventDPhiv2new4060->Fill(2.5, TMath::Cos(2.*(evPlAngTPCp-evPlAngTPCn))); // tpcp - tpcn
 880
 881       hEvPlaneTPCvsEvPVz4060 ->Fill(  evPlAngTPC,evPlAngV0);
 882     }
 883
 884     if(percentile>=0 || percentile<=7.5)
 885       hEvPlaneTPCvsEvPVz075->Fill(  evPlAngTPC,evPlAngV0);
 886     if(percentile>=15 || percentile<=30)
 887       hEvPlaneTPCvsEvPVz1530->Fill(  evPlAngTPC,evPlAngV0);
 888     if(percentile>=30 || percentile<=50)
 889       hEvPlaneTPCvsEvPVz3050->Fill(  evPlAngTPC,evPlAngV0);
 890
 891     //====================================================================================================================
 892
 893     // To remove auto-correlation
 894     TVector2 *q = 0x0;
 895     TVector2 *qsub1=0x0;
 896     TVector2 *qsub2=0x0;
 897     qsub1 = pl->GetQsub1();
 898     qsub2 = pl->GetQsub2();
 899
 900     q = pl->GetQVector();
 901
 902     trpangleTPC      = evPlAngTPC;
 903     trpangleVZERO[0] = evPlAngV0;
 904     trpangleVZERO[1] = evPlAngV0A;
 905     trpangleVZERO[2] = evPlAngV0C;
 906
 907     // cout<<"rpangle TPC: "<<rpangleTPC<<endl;
 908
 909     Int_t isTOF=0;
 910     Int_t isoutTPC=0;
 911     //  cout<<"TRack number MC "<<TrackNumber<<endl;
 912
 913     for (Int_t j=0; j<TrackNumber; j++) { //loop on tracks
 914
 915       AliESDtrack *esdtrack=lESDevent->GetTrack(j);
 916       if (!fESDtrackCuts->AcceptTrack(esdtrack)) continue;
 917
 918       status  = (ULong_t)esdtrack->GetStatus();
 919
 920       isTPC    = (((status) & AliESDtrack::kTPCin)   != 0);
 921       isTOF    = ((((status) & AliESDtrack::kTOFout) != 0) && (((status) & AliESDtrack::kTIME) != 0));
 922       isoutTPC = (((status) & AliESDtrack::kTPCout)  != 0);
 923
 924       TPCSignal=esdtrack->GetTPCsignal();
 925
 926       if (TPCSignal<10)continue;
 927       if (TPCSignal>1000)continue;
 928       if (!isTPC)continue;
 929
 930       if(!esdtrack->GetTPCInnerParam())continue;
 931       AliExternalTrackParam trackIn(*esdtrack->GetInnerParam());
 932       pinTPC= trackIn.GetP();
 933
 934       fhBB->Fill(pinTPC*esdtrack->GetSign(),TPCSignal);
 935
 936       if(isTOF){
 937         if(!esdtrack->GetOuterParam())continue;
 938         AliExternalTrackParam trackOut(*esdtrack->GetOuterParam());
 939         poutTPC = trackOut.GetP();
 940         fhTOF->Fill(poutTPC*esdtrack->GetSign(),(esdtrack->GetIntegratedLength()/esdtrack->GetTOFsignal())/2.99792458e-2);
 941
 942       }
 943
 944       Int_t   fIdxInt[200]; //dummy array
 945       Int_t   nClustersTPC = esdtrack->GetTPCclusters(fIdxInt);
 946       Float_t chi2PerClusterTPC = esdtrack->GetTPCchi2()/(Float_t)(nClustersTPC);
 947
 948       esdtrack->GetImpactParameters(impactXY, impactZ);
 949
 950       Float_t deutExp  = -999;
 951       Float_t tritExp  = -999;
 952       Float_t hel3Exp  = -999;
 953
 954       if(fDataType == "REAL"){
 955         deutExp  = AliExternalTrackParam::BetheBlochAleph(pinTPC/(0.938*2),1.45802,27.4992,4.00313e-15,2.48485,8.31768);
 956         tritExp  = AliExternalTrackParam::BetheBlochAleph(pinTPC/(0.938*3),1.45802,27.4992,4.00313e-15,2.48485,8.31768);
 957         hel3Exp  = 4*AliExternalTrackParam::BetheBlochAleph(2*pinTPC/(0.938*3),1.74962,27.4992,4.00313e-15,2.42485,8.31768);
 958       }
 959
 960       if(fDataType == "SIM"){
 961         Double_t parMC[5] = {1.17329, 27.4992, 4.00313e-15, 2.1204316, 4.1373729}; // NEW!!!
 962         deutExp = AliExternalTrackParam::BetheBlochAleph(pinTPC/(0.938*2),parMC[0],parMC[1],parMC[2],parMC[3],parMC[4]);
 963         tritExp = AliExternalTrackParam::BetheBlochAleph(pinTPC/(0.938*3),parMC[0],parMC[1],parMC[2],parMC[3],parMC[4]);
 964         hel3Exp = 4*AliExternalTrackParam::BetheBlochAleph(2*pinTPC/(0.938*3),parMC[0],parMC[1],parMC[2],parMC[3],parMC[4]);
 965       }
 966
 967       Double_t pullTPC     = (TPCSignal - deutExp)/(0.07*deutExp);
 968       Double_t pullTPChel3 = (TPCSignal - hel3Exp)/(0.07*hel3Exp);
 969       Double_t pullTPCtrit = (TPCSignal - tritExp)/(0.07*tritExp);
 970
 971       tPulls[0] = pullTPC;
 972       tPulls[1] = pullTPChel3;
 973       tPulls[2] = pullTPCtrit;
 974
 975       //Fill the tree
 976
 977       tCentrality[0] = percentile;
 978       tCentrality[1] = eventtype;
 979
 980
 981       tMomentum[0] = esdtrack->Px();
 982       tMomentum[1] = esdtrack->Py();
 983       tMomentum[2] = esdtrack->Pz();
 984
 985       //Corrected momentum from Alexander
 986       Double_t pT = esdtrack->Pt()/(1 - 0.333303/TMath::Power(esdtrack->Pt() + 0.651111, 5.27268));
 987       tMomentum[3] = pT;
 988
 989       tDCA[0]      = impactXY;
 990       tDCA[1]      = impactZ;
 991
 992       tisTOF[0] = isTOF;
 993       tisTOF[1] = isoutTPC;
 994
 995       Double_t p   = esdtrack->P();
 996       Double_t tof = esdtrack->GetTOFsignal()-fPIDResponse->GetTOFResponse().GetStartTime(p);
 997
 998       tTOFtrack[0] = poutTPC;
 999       tTOFtrack[1] = tof;                             //ps = Time
1000       tTOFtrack[2] = esdtrack->GetIntegratedLength(); //cm
1001
1002       tCharge =  esdtrack->Charge();
1003       tPhi    =  esdtrack->Phi();
1004
1005       Float_t beta = 0;
1006       Float_t gamma = 0;
1007       Float_t deltaMass = 0;
1008
1009       if(fDataType == "REAL"){
1010         if(TMath::Abs(pinTPC) < 6 && TMath::Abs(pullTPC) < 3){
1011
1012           fhBBDeu->Fill(pinTPC*esdtrack->GetSign(),TPCSignal);
1013
1014           if(pinTPC < 1.0 && fFillNtuple == kTRUE)
1015             fNtuple1->Fill();
1016
1017           if(tTOFtrack[1] > 0 ){
1018             beta = tTOFtrack[2]/(tTOFtrack[1] * 2.99792457999999984e-02);
1019             gamma = 1/TMath::Sqrt(1 - beta*beta);
1020             deltaMass = poutTPC/TMath::Sqrt(gamma*gamma - 1) - 1.8756;
1021             fhMassTOF->Fill(deltaMass);
1022           }
1023
1024           if(pinTPC >= 1.0 && tTOFtrack[1] > 0 && TMath::Abs(deltaMass)<0.5 && fFillNtuple == kTRUE )
1025             fNtuple1->Fill();
1026         }
1027
1028         if(pinTPC < 10. && TMath::Abs(pullTPChel3) < 3){
1029           if(fFillNtuple == kTRUE)
1030             fNtuple1->Fill();
1031           fhBBDeu->Fill(pinTPC*esdtrack->GetSign(),TPCSignal);
1032         }
1033
1034
1035         //Fill also things for flow package
1036
1037         if(TMath::Abs(pullTPC)<2){
1038
1039           fhPtDeu->Fill(esdtrack->Pt(),pT);
1040
1041           //Remove AutoCorrelation
1042           trpangleTPC = GetEventPlaneForCandidate(esdtrack,q,pl,qsub1,qsub2);
1043
1044           Float_t deltaphiTPC=2*GetPhi0Pi(tPhi-trpangleTPC);
1045           Float_t deltaphiV0 =2*GetPhi0Pi(tPhi-trpangleVZERO[0]);
1046
1047
1048           if(pinTPC < 1.0){
1049
1050             //Here
1051
1052             //==================================================================
1053             //----------------------Flow of Inclusive Particles --------------------------------------------------------
1054
1055             AliFlowTrack *sTrack = new AliFlowTrack();
1056             sTrack->Set(esdtrack);
1057             sTrack->SetID(esdtrack->GetID());
1058             sTrack->SetForRPSelection(kTRUE);
1059             sTrack->SetForPOISelection(kTRUE);
1060
1061             for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs)
1062               {
1063                 //   cout << " no of rps " << iRPs << endl;
1064                 AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
1065                 if (!iRP) continue;
1066                 if (!iRP->InRPSelection()) continue;
1067                 if( sTrack->GetID() == iRP->GetID())
1068                   {
1069                     if(fDebug) printf(" was in RP set");
1070                     //       cout << sTrack->GetID() <<"   ==  " << iRP->GetID() << " was in RP set" <<endl;
1071                     iRP->SetForRPSelection(kFALSE);
1072                     // fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
1073                   }
1074               } //end of for loop on RPs
1075             fFlowEvent->InsertTrack(((AliFlowTrack*) sTrack));
1076             fFlowEvent->SetNumberOfPOIs(fFlowEvent->GetNumberOfPOIs()+1);
1077
1078
1079             //=================================================================
1080
1081             //  if(tCharge > 0){
1082
1083             if(tCentrality[0]>0 && tCentrality[0]<7.5){
1084               hCos2DeltaPhivsPt075           ->Fill(TMath::Cos(deltaphiTPC),tMomentum[3]*tCharge);
1085               hCos2DeltaPhiVZEROvsPt075      ->Fill(TMath::Cos(deltaphiV0) ,tMomentum[3]*tCharge);
1086
1087             }
1088
1089             if(tCentrality[0]>0 && tCentrality[0]<5){
1090               hCos2DeltaPhivsPt05           ->Fill(TMath::Cos(deltaphiTPC),tMomentum[3]*tCharge);
1091               hCos2DeltaPhiVZEROvsPt05      ->Fill(TMath::Cos(deltaphiV0) ,tMomentum[3]*tCharge);
1092             }
1093
1094             if(tCentrality[0]>15 && tCentrality[0]<30){
1095               hCos2DeltaPhivsPt1530           ->Fill(TMath::Cos(deltaphiTPC),tMomentum[3]*tCharge);
1096               hCos2DeltaPhiVZEROvsPt1530      ->Fill(TMath::Cos(deltaphiV0) ,tMomentum[3]*tCharge);
1097             }
1098
1099             if(tCentrality[0]>20 && tCentrality[0]<40){
1100               hCos2DeltaPhivsPt2040           ->Fill(TMath::Cos(deltaphiTPC),tMomentum[3]*tCharge);
1101               hCos2DeltaPhiVZEROvsPt2040      ->Fill(TMath::Cos(deltaphiV0) ,tMomentum[3]*tCharge);
1102             }
1103
1104             if(tCentrality[0]>30 && tCentrality[0]<50){
1105               hCos2DeltaPhivsPt3050           ->Fill(TMath::Cos(deltaphiTPC),tMomentum[3]*tCharge);
1106               hCos2DeltaPhiVZEROvsPt3050      ->Fill(TMath::Cos(deltaphiV0) ,tMomentum[3]*tCharge);
1107             }
1108
1109             if(tCentrality[0]>40 && tCentrality[0]<60){
1110               hCos2DeltaPhivsPt4060           ->Fill(TMath::Cos(deltaphiTPC),tMomentum[3]*tCharge);
1111               hCos2DeltaPhiVZEROvsPt4060      ->Fill(TMath::Cos(deltaphiV0) ,tMomentum[3]*tCharge);
1112             }
1113           }
1114
1115           if(pinTPC > 1.0 && pinTPC  < 6.0 && tTOFtrack[1] > 0 && TMath::Abs(deltaMass)<0.5){
1116
1117
1118             if(tCentrality[0]>0 && tCentrality[0]<7.5){
1119               hCos2DeltaPhivsPt075           ->Fill(TMath::Cos(deltaphiTPC),tMomentum[3]*tCharge);
1120               hCos2DeltaPhiVZEROvsPt075      ->Fill(TMath::Cos(deltaphiV0) ,tMomentum[3]*tCharge);
1121
1122             }
1123
1124             if(tCentrality[0]>0 && tCentrality[0]<5){
1125               hCos2DeltaPhivsPt05           ->Fill(TMath::Cos(deltaphiTPC),tMomentum[3]*tCharge);
1126               hCos2DeltaPhiVZEROvsPt05      ->Fill(TMath::Cos(deltaphiV0) ,tMomentum[3]*tCharge);
1127             }
1128
1129             if(tCentrality[0]>15 && tCentrality[0]<30){
1130               hCos2DeltaPhivsPt1530           ->Fill(TMath::Cos(deltaphiTPC),tMomentum[3]*tCharge);
1131               hCos2DeltaPhiVZEROvsPt1530      ->Fill(TMath::Cos(deltaphiV0) ,tMomentum[3]*tCharge);
1132             }
1133
1134             if(tCentrality[0]>20 && tCentrality[0]<40){
1135               hCos2DeltaPhivsPt2040           ->Fill(TMath::Cos(deltaphiTPC),tMomentum[3]*tCharge);
1136               hCos2DeltaPhiVZEROvsPt2040      ->Fill(TMath::Cos(deltaphiV0) ,tMomentum[3]*tCharge);
1137             }
1138
1139             if(tCentrality[0]>30 && tCentrality[0]<50){
1140               hCos2DeltaPhivsPt3050           ->Fill(TMath::Cos(deltaphiTPC),tMomentum[3]*tCharge);
1141               hCos2DeltaPhiVZEROvsPt3050      ->Fill(TMath::Cos(deltaphiV0) ,tMomentum[3]*tCharge);
1142             }
1143
1144             if(tCentrality[0]>40 && tCentrality[0]<60){
1145               hCos2DeltaPhivsPt4060           ->Fill(TMath::Cos(deltaphiTPC),tMomentum[3]*tCharge);
1146               hCos2DeltaPhiVZEROvsPt4060      ->Fill(TMath::Cos(deltaphiV0) ,tMomentum[3]*tCharge);
1147             }
1148           }
1149         }
1150       }
1151
1152       if(fDataType == "SIM"){
1153
1154         Int_t  label = TMath::Abs(esdtrack->GetLabel());
1155         TParticle * part = stack->Particle(label);
1156         Int_t PDGCode=part->GetPdgCode();
1157
1158         Int_t motherPDG=0;
1159
1160         Int_t mumid = part->GetFirstMother();
1161         if(mumid>-1){
1162           TParticle *mother=(TParticle*)stack->Particle(mumid);
1163           motherPDG = mother->GetPdgCode();
1164         }
1165
1166         if( PDGCode == 1000010020 || PDGCode == -1000010020 || PDGCode == 1000020030 || PDGCode == -1000020030){
1167
1168           tPDGCode=PDGCode;
1169           tPDGCodeMum = motherPDG;
1170
1171           tIsPrimaryTr      = stack->IsPhysicalPrimary(label);
1172           tIsSecondaryTr[0] = stack->IsSecondaryFromMaterial(label);
1173           tIsSecondaryTr[1] = stack->IsSecondaryFromWeakDecay(label);
1174
1175           fNtuple1->Fill();
1176           fhPtDeu->Fill(esdtrack->Pt(),pT);
1177
1178           if(tTOFtrack[1] > 0){
1179             beta = tTOFtrack[2]/(tTOFtrack[1] * 2.99792457999999984e-02);
1180             gamma = 1/TMath::Sqrt(1 - beta*beta);
1181             fhMassTOF->Fill(poutTPC/TMath::Sqrt(gamma*gamma - 1) - 1.8756);
1182           }
1183         }
1184       }
1185
1186     }   //track
1187
1188     //==END RECONSTRUCTION==
1189
1190
1191     // MC truth
1192
1193     if(fDataType == "SIM"){
1194
1195       for (Int_t iMC=0; iMC<stack->GetNtrack(); iMC++){
1196
1197         const TParticle *tparticle = stack->Particle(iMC);
1198         Long_t PDGCode = tparticle->GetPdgCode();
1199
1200         Double_t eta = tparticle->Eta();
1201         Double_t pt  = tparticle->Pt();
1202         Double_t rap = tparticle->Y();
1203
1204         //check which particle it is
1205
1206         Float_t codemoth = 0;
1207
1208         Int_t indexMoth=tparticle->GetFirstMother();
1209
1210         if(indexMoth>=0){
1211           TParticle* moth = stack->Particle(indexMoth);
1212           codemoth = TMath::Abs(moth->GetPdgCode());
1213         }
1214
1215         //d, 3He
1216
1217         if( PDGCode == 1000010020 || PDGCode == -1000010020 || PDGCode == 1000020030 || PDGCode == -1000020030){
1218
1219
1220
1221           tCentralityMC = percentile;
1222
1223           tVertexCoordMC[0] = vtx->GetXv();
1224           tVertexCoordMC[1] = vtx->GetYv();
1225           tVertexCoordMC[2] = vtx->GetZv();
1226
1227           tVertexCoordMC[0] = tparticle->Px();
1228           tVertexCoordMC[1] = tparticle->Py();
1229           tVertexCoordMC[2] = tparticle->Pz();
1230
1231           tPDGCodeMC = PDGCode;
1232           tPDGCodeMumMC = codemoth;
1233           tEtaMC = eta;
1234
1235           tPtMC = pt;
1236           tYMC  = rap;
1237
1238           tIsPrimary      = stack->IsPhysicalPrimary(iMC);
1239           tIsSecondary[0] = stack->IsSecondaryFromMaterial(iMC);
1240           tIsSecondary[1] = stack->IsSecondaryFromWeakDecay(iMC);
1241
1242
1243           fNtuple2->Fill();
1244         }
1245       }
1246
1247     }
1248   }
1249
1250   PostData(1, fListHist);
1251   PostData(2,  fNtuple1);
1252   PostData(3,  fNtuple2);
1253
1254   PostData(4, fFlowEvent);
1255
1256 } //end userexec
1257
1258
1259 //________________________________________________________________________
1260
1261 void AliAnalysisTaskNucleiv2::Terminate(Option_t *)
1262 {
1263   // Draw result to the screen
1264   // Called once at the end of the query
1265 }
1266