PWG/FLOW/Tasks/AliAnalysisTaskFlowStrange.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008,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 // AliAnalysisTaskFlowStrange:
  18 // Analysis task to select K0/Lambda candidates for flow analysis.
  19 // Authors: Cristian Ivan (civan@cern.ch)
  20 //          Carlos Perez  (cperez@cern.ch)
  21 //          Pawel Debski  (pdebski@cern.ch)
  22 //////////////////////////////////////////////////////
  23
  24 #include "TChain.h"
  25 #include "TList.h"
  26 #include "TH1D.h"
  27 #include "TH2D.h"
  28 #include "TH3D.h"
  29 #include "TF1.h"
  30 #include "TProfile.h"
  31 #include "TProfile2D.h"
  32 #include "TVector3.h"
  33 #include "TStopwatch.h"
  34 #include "TFile.h"
  35
  36 #include "TRandom3.h"
  37
  38 #include "AliAnalysisManager.h"
  39 #include "AliInputEventHandler.h"
  40
  41 #include "AliVVertex.h"
  42 #include "AliVVZERO.h"
  43 #include "AliStack.h"
  44 #include "AliMCEvent.h"
  45
  46 #include "AliESDEvent.h"
  47 #include "AliESDtrack.h"
  48 #include "AliESDVertex.h"
  49 #include "AliESDv0.h"
  50 #include "AliESDtrackCuts.h"
  51
  52 #include "AliAODEvent.h"
  53 #include "AliAODTrack.h"
  54 #include "AliAODVertex.h"
  55 #include "AliAODv0.h"
  56 #include "AliAODTracklets.h"
  57 #include "AliAODHeader.h"
  58
  59 #include "AliAODMCHeader.h"
  60 #include "AliAODMCParticle.h"
  61 #include "TClonesArray.h"
  62 #include "TDatabasePDG.h"
  63 #include "TParticlePDG.h"
  64
  65 #include "TMath.h"
  66 #include "TObjArray.h"
  67 #include "AliFlowCandidateTrack.h"
  68
  69 #include "AliFlowTrackCuts.h"
  70 #include "AliFlowEventCuts.h"
  71 #include "AliFlowEvent.h"
  72 #include "AliFlowBayesianPID.h"
  73 #include "AliFlowCommonConstants.h"
  74 #include "AliFlowVector.h"
  75
  76 #include "AliAnalysisTaskFlowStrange.h"
  77
  78 ClassImp(AliAnalysisTaskFlowStrange)
  79
  80 //=======================================================================
  81 AliAnalysisTaskFlowStrange::AliAnalysisTaskFlowStrange() :
  82   AliAnalysisTaskSE(),
  83   fPIDResponse(NULL),
  84   fFB1(NULL),
  85   fFB1024(NULL),
  86   fTPCevent(NULL),
  87   fVZEevent(NULL),
  88   fCandidates(NULL),
  89   fList(NULL),
  90   fRunNumber(-1),
  91   fDebug(0),
  92   fQAlevel(0),
  93   fReadESD(kFALSE),
  94   fReadMC(kFALSE),
  95   fAvoidExec(kFALSE),
  96   fSkipSelection(kFALSE),
  97   fSkipFlow(kFALSE),
  98   fSkipDHcorr(kTRUE),
  99   fUseFP(kFALSE),
 100   fRunOnpA(kFALSE),
 101   fRunOnpp(kFALSE),
 102   fExtraEventRejection(kFALSE),
 103   fCentMethod("V0MTRK"),
 104   fCentPerMin(0),
 105   fCentPerMax(100),
 106   fThisCent(-1.0),
 107   fExcludeTPCEdges(kFALSE),
 108   fSpecie(0),
 109   fOnline(kFALSE),
 110   fHomemade(kFALSE),
 111   fWhichPsi(1),
 112   fVZEsave(kFALSE),
 113   fVZEload(NULL),
 114   fVZEResponse(NULL),
 115   fVZEmb(kFALSE),
 116   fVZEByDisk(kTRUE),
 117   fVZECa(0),
 118   fVZECb(3),
 119   fVZEAa(0),
 120   fVZEAb(3),
 121   fVZEQA(NULL),
 122   fPsi2(0.0),
 123   fMCEP(0.0),
 124   fMassBins(0),
 125   fMinMass(0.0),
 126   fMaxMass(0.0),
 127   fRFPFilterBit(1),
 128   fRFPminPt(0.2),
 129   fRFPmaxPt(5.0),
 130   fRFPminEta(-0.8),
 131   fRFPmaxEta(+0.8),
 132   fRFPTPCsignal(10.0),
 133   fRFPmaxIPxy(2.4),
 134   fRFPmaxIPz(3.2),
 135   fRFPTPCncls(70),
 136   fDecayMass(0.0),
 137   fDecayPhi(0.0),
 138   fDecayEta(0.0),
 139   fDecayPt(0.0),
 140   fDecayDCAdaughters(0.0),
 141   fDecayCosinePointingAngleXY(0.0),
 142   fDecayRadXY(0.0),
 143   fDecayDecayLength(0.0),
 144   fDecayQt(0.0),
 145   fDecayAlpha(0.0),
 146   fDecayRapidity(0.0),
 147   fDecayProductIPXY(0.0),
 148   fDecayIDneg(-1),
 149   fDecayIDpos(-1),
 150   fDecayID(-1),
 151   fDecayMinEta(0.0),
 152   fDecayMaxEta(0.0),
 153   fDecayMinPt(0.0),
 154   fDecayMaxDCAdaughters(0.0),
 155   fDecayMinCosinePointingAngleXY(0.0),
 156   fDecayMinQt(0.0),
 157   fDecayAPCutPie(kTRUE),
 158   fDecayMinRadXY(0.0),
 159   fDecayMaxDecayLength(0.0),
 160   fDecayMaxProductIPXY(0.0),
 161   fDecayMaxRapidity(0.0),
 162   fDaughterPhi(0.0),
 163   fDaughterEta(0.0),
 164   fDaughterPt(0.0),
 165   fDaughterNClsTPC(0),
 166   fDaughterCharge(0),
 167   fDaughterNFClsTPC(0),
 168   fDaughterNSClsTPC(0),
 169   fDaughterChi2PerNClsTPC(0.0),
 170   fDaughterXRows(0.0),
 171   fDaughterImpactParameterXY(0.0),
 172   fDaughterImpactParameterZ(0.0),
 173   fDaughterStatus(0),
 174   fDaughterNSigmaPID(0.0),
 175   fDaughterKinkIndex(0),
 176   fDaughterUnTag(kTRUE),
 177   fDaughterMinEta(0.0),
 178   fDaughterMaxEta(0.0),
 179   fDaughterMinPt(0.0),
 180   fDaughterMinNClsTPC(0),
 181   fDaughterMinXRows(0),
 182   fDaughterMaxChi2PerNClsTPC(0.0),
 183   fDaughterMinXRowsOverNClsFTPC(0.0),
 184   fDaughterMinImpactParameterXY(0.0),
 185   fDaughterMaxNSigmaPID(0.0) {
 186   //ctor
 187 }
 188 //=======================================================================
 189 AliAnalysisTaskFlowStrange::AliAnalysisTaskFlowStrange(const char *name) :
 190   AliAnalysisTaskSE(name),
 191   fPIDResponse(NULL),
 192   fFB1(NULL),
 193   fFB1024(NULL),
 194   fTPCevent(NULL),
 195   fVZEevent(NULL),
 196   fCandidates(NULL),
 197   fList(NULL),
 198   fRunNumber(-1),
 199   fDebug(0),
 200   fQAlevel(0),
 201   fReadESD(kFALSE),
 202   fReadMC(kFALSE),
 203   fAvoidExec(kFALSE),
 204   fSkipSelection(kFALSE),
 205   fSkipFlow(kFALSE),
 206   fSkipDHcorr(kTRUE),
 207   fUseFP(kFALSE),
 208   fRunOnpA(kFALSE),
 209   fRunOnpp(kFALSE),
 210   fExtraEventRejection(kFALSE),
 211   fCentMethod("V0MTRK"),
 212   fCentPerMin(0),
 213   fCentPerMax(100),
 214   fThisCent(-1.0),
 215   fExcludeTPCEdges(kFALSE),
 216   fSpecie(0),
 217   fOnline(kFALSE),
 218   fHomemade(kFALSE),
 219   fWhichPsi(1),
 220   fVZEsave(kFALSE),
 221   fVZEload(NULL),
 222   fVZEResponse(NULL),
 223   fVZEmb(kFALSE),
 224   fVZEByDisk(kTRUE),
 225   fVZECa(0),
 226   fVZECb(3),
 227   fVZEAa(0),
 228   fVZEAb(3),
 229   fVZEQA(NULL),
 230   fPsi2(0.0),
 231   fMCEP(0.0),
 232   fMassBins(0),
 233   fMinMass(0.0),
 234   fMaxMass(0.0),
 235   fRFPFilterBit(1),
 236   fRFPminPt(0.2),
 237   fRFPmaxPt(5.0),
 238   fRFPminEta(-0.8),
 239   fRFPmaxEta(+0.8),
 240   fRFPTPCsignal(10.0),
 241   fRFPmaxIPxy(2.4),
 242   fRFPmaxIPz(3.2),
 243   fRFPTPCncls(70),
 244   fDecayMass(0.0),
 245   fDecayPhi(0.0),
 246   fDecayEta(0.0),
 247   fDecayPt(0.0),
 248   fDecayDCAdaughters(0.0),
 249   fDecayCosinePointingAngleXY(0.0),
 250   fDecayRadXY(0.0),
 251   fDecayDecayLength(0.0),
 252   fDecayQt(0.0),
 253   fDecayAlpha(0.0),
 254   fDecayRapidity(0.0),
 255   fDecayProductIPXY(0.0),
 256   fDecayIDneg(-1),
 257   fDecayIDpos(-1),
 258   fDecayID(-1),
 259   fDecayMinEta(0.0),
 260   fDecayMaxEta(0.0),
 261   fDecayMinPt(0.0),
 262   fDecayMaxDCAdaughters(0.0),
 263   fDecayMinCosinePointingAngleXY(0.0),
 264   fDecayMinQt(0.0),
 265   fDecayAPCutPie(kTRUE),
 266   fDecayMinRadXY(0.0),
 267   fDecayMaxDecayLength(0.0),
 268   fDecayMaxProductIPXY(0.0),
 269   fDecayMaxRapidity(0.0),
 270   fDaughterPhi(0.0),
 271   fDaughterEta(0.0),
 272   fDaughterPt(0.0),
 273   fDaughterNClsTPC(0),
 274   fDaughterCharge(0),
 275   fDaughterNFClsTPC(0),
 276   fDaughterNSClsTPC(0),
 277   fDaughterChi2PerNClsTPC(0.0),
 278   fDaughterXRows(0.0),
 279   fDaughterImpactParameterXY(0.0),
 280   fDaughterImpactParameterZ(0.0),
 281   fDaughterStatus(0),
 282   fDaughterNSigmaPID(0.0),
 283   fDaughterKinkIndex(0),
 284   fDaughterUnTag(kTRUE),
 285   fDaughterMinEta(0.0),
 286   fDaughterMaxEta(0.0),
 287   fDaughterMinPt(0.0),
 288   fDaughterMinNClsTPC(0),
 289   fDaughterMinXRows(0),
 290   fDaughterMaxChi2PerNClsTPC(0.0),
 291   fDaughterMinXRowsOverNClsFTPC(0.0),
 292   fDaughterMinImpactParameterXY(0.0),
 293   fDaughterMaxNSigmaPID(0.0) {
 294   //ctor
 295   DefineInput( 0,TChain::Class());
 296   DefineOutput(1,TList::Class());
 297   DefineOutput(2,AliFlowEventSimple::Class()); // TPC object
 298   DefineOutput(3,AliFlowEventSimple::Class()); // VZE object
 299 }
 300 //=======================================================================
 301 AliAnalysisTaskFlowStrange::~AliAnalysisTaskFlowStrange() {
 302   //dtor
 303   if (fCandidates) delete fCandidates;
 304   if (fTPCevent)   delete fTPCevent;
 305   if (fVZEevent)   delete fVZEevent;
 306   if (fList)       delete fList;
 307 }
 308 //=======================================================================
 309 void AliAnalysisTaskFlowStrange::UserCreateOutputObjects() {
 310   //UserCreateOutputObjects
 311   fList=new TList();
 312   fList->SetOwner();
 313   AddQAEvents();
 314   AddQACandidates();
 315
 316   if(fReadESD) MakeFilterBits();
 317
 318   AliFlowCommonConstants *cc = AliFlowCommonConstants::GetMaster();
 319   cc->SetNbinsMult(3000); cc->SetMultMin(0);   cc->SetMultMax(30000);
 320   cc->SetNbinsPt(200); cc->SetPtMin(0.0);   cc->SetPtMax(20.0);
 321   cc->SetNbinsPhi(100);  cc->SetPhiMin(0.0);  cc->SetPhiMax(TMath::TwoPi());
 322   cc->SetNbinsEta(100);  cc->SetEtaMin(-5.0); cc->SetEtaMax(+5.0);
 323   cc->SetNbinsQ(100);    cc->SetQMin(0.0);    cc->SetQMax(3.0);
 324   cc->SetNbinsMass(fMassBins);
 325   cc->SetMassMin(fMinMass);
 326   cc->SetMassMax(fMaxMass);
 327
 328   //loading pid response
 329   AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
 330   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 331   fPIDResponse = inputHandler->GetPIDResponse();
 332
 333   fTPCevent = new AliFlowEvent(100);
 334   fVZEevent = new AliFlowEvent(100);
 335
 336   //array of candidates
 337   fCandidates = new TObjArray(100);
 338   fCandidates->SetOwner();
 339
 340   PostData(1,fList);
 341   if(fUseFP) { // for connection to the flow package
 342     PostData(2,fTPCevent);
 343     PostData(3,fVZEevent);
 344   }
 345
 346   gRandom->SetSeed();
 347 }
 348 //=======================================================================
 349 void AliAnalysisTaskFlowStrange::AddQAEvents() {
 350   // function to add event qa
 351   TH1D *tH1D;
 352   TProfile *tProfile;
 353   TList *tQAEvents=new TList();
 354   tQAEvents->SetName("Event");
 355   tQAEvents->SetOwner();
 356   tH1D = new TH1D("Events","Number of Events",6,0,6); tQAEvents->Add(tH1D);
 357   tH1D->GetXaxis()->SetBinLabel(1,"exec");
 358   tH1D->GetXaxis()->SetBinLabel(2,"userexec");
 359   tH1D->GetXaxis()->SetBinLabel(3,"reached");
 360   tH1D->GetXaxis()->SetBinLabel(4,"selected");
 361   tH1D->GetXaxis()->SetBinLabel(5,"rejectedByLowQw");
 362   tH1D->GetXaxis()->SetBinLabel(6,"rejectedByErrorLoadVZEcal");
 363   tProfile = new TProfile("Configuration","Configuration",20,0,20); tQAEvents->Add(tProfile);
 364   tProfile->Fill( 0.5,fCentPerMin,1); tProfile->GetXaxis()->SetBinLabel( 1,"fCentPerMin");
 365   tProfile->Fill( 1.5,fCentPerMax,1); tProfile->GetXaxis()->SetBinLabel( 2,"fCentPerMax");
 366   tProfile->Fill( 2.5,fDaughterMinEta,1);               tProfile->GetXaxis()->SetBinLabel( 3,"fDaughterMinEta");
 367   tProfile->Fill( 3.5,fDaughterMaxEta,1);               tProfile->GetXaxis()->SetBinLabel( 4,"fDaughterMaxEta");
 368   tProfile->Fill( 4.5,fDaughterMinPt,1);                tProfile->GetXaxis()->SetBinLabel( 5,"fDaughterMinPt");
 369   tProfile->Fill( 5.5,fDaughterMinNClsTPC,1);           tProfile->GetXaxis()->SetBinLabel( 6,"fDaughterMinNClsTPC");
 370   tProfile->Fill( 6.5,fDaughterMaxChi2PerNClsTPC,1);    tProfile->GetXaxis()->SetBinLabel( 7,"fDaughterMaxChi2PerNClsTPC");
 371   tProfile->Fill( 7.5,fDaughterMinXRowsOverNClsFTPC,1); tProfile->GetXaxis()->SetBinLabel( 8,"fDaughterMinXRowsOverNClsFTPC");
 372   tProfile->Fill( 8.5,fDaughterMinImpactParameterXY,1); tProfile->GetXaxis()->SetBinLabel( 9,"fDaughterMinImpactParameterXY");
 373   tProfile->Fill( 9.5,fDaughterMaxNSigmaPID,1);         tProfile->GetXaxis()->SetBinLabel(10,"fDaughterMaxNSigmaPID");
 374   tProfile->Fill(10.5,fDecayMaxDCAdaughters,1);          tProfile->GetXaxis()->SetBinLabel(11,"fDecayMaxDCAdaughters");
 375   tProfile->Fill(11.5,fDecayMinCosinePointingAngleXY,1); tProfile->GetXaxis()->SetBinLabel(12,"fDecayMinCosinePointingAngleXY");
 376   tProfile->Fill(12.5,fDecayMinQt,1);                    tProfile->GetXaxis()->SetBinLabel(13,"fDecayMinQt");
 377   tProfile->Fill(13.5,fDecayMinRadXY,1);                 tProfile->GetXaxis()->SetBinLabel(14,"fDecayMinRadXY");
 378   tProfile->Fill(14.5,fDecayMaxDecayLength,1);           tProfile->GetXaxis()->SetBinLabel(15,"fDecayMaxDecayLength");
 379   tProfile->Fill(15.5,fDecayMaxProductIPXY,1);           tProfile->GetXaxis()->SetBinLabel(16,"fDecayMaxProductIPXY");
 380   tProfile->Fill(16.5,fDecayMaxRapidity,1);              tProfile->GetXaxis()->SetBinLabel(17,"fDecayMaxRapidity");
 381   tProfile->Fill(17.5,fDecayMinEta,1);                   tProfile->GetXaxis()->SetBinLabel(18,"fDecayMinEta");
 382   tProfile->Fill(18.5,fDecayMaxEta,1);                   tProfile->GetXaxis()->SetBinLabel(19,"fDecayMaxEta");
 383   tProfile->Fill(19.5,fDecayMinPt,1);                    tProfile->GetXaxis()->SetBinLabel(20,"fDecayMinPt");
 384
 385   tH1D = new TH1D("POI","POIs;multiplicity",800,0,800);         tQAEvents->Add(tH1D);
 386   tH1D = new TH1D("UNTAG","UNTAG;Untagged Daughters",800,0,800);tQAEvents->Add(tH1D);
 387   tH1D = new TH1D("RealTime","RealTime;LogT sec",2000,-10,+10); tQAEvents->Add(tH1D);
 388   fList->Add(tQAEvents);
 389   AddEventSpy();
 390   AddMakeQSpy();
 391 }
 392 //=======================================================================
 393 void AliAnalysisTaskFlowStrange::AddEventSpy() {
 394   TH1D *tH1D;
 395   TH2D *tH2D;
 396   TList *tList=new TList();
 397   tList->SetName("EventSpy");
 398   tList->SetOwner();
 399   if(fQAlevel>0) {
 400     tH2D = new TH2D("VTXZ","VTXZ;Global||SPD;SPD",60,-25,+25,60,-25,+25); tList->Add( tH2D );
 401     tH2D = new TH2D("CCCC","CCCC;V0M;TRK",60,-10,110,60,-10,110);         tList->Add( tH2D );
 402     tH2D = new TH2D("REFM","REFM;TPC;GLOBAL",100,0,3000,100,0,3000);      tList->Add( tH2D );
 403     if(fReadMC) {
 404       tH1D = new TH1D("MCCC","MCCC;Xsection",100,-10,110); tList->Add( tH1D );
 405       tH1D = new TH1D("MCEP","MCEP;MCEP",100,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
 406     }
 407   }
 408   fList->Add(tList);
 409 }
 410 //=======================================================================
 411 void AliAnalysisTaskFlowStrange::AddMakeQSpy() {
 412   if(fSkipFlow) return;
 413   TH1D *tH1D;
 414   TH2D *tH2D;
 415   TProfile *tPF1;
 416   TList *tList=new TList();
 417   tList->SetName("MakeQSpy");
 418   tList->SetOwner();
 419   tH1D = new TH1D("RFPTPC","TPC Refrence Multiplicity;multiplicity",3000,0,3000);     tList->Add( tH1D );
 420   tH1D = new TH1D("RFPVZE","VZERO Reference Multiplicity;multiplicity",3000,0,30000); tList->Add( tH1D );
 421   tH2D = new TH2D("TPCAllPhiEta","TPCall;Phi;Eta",180,0,TMath::TwoPi(),80,-0.9,+0.9); tList->Add( tH2D );
 422   tH2D = new TH2D("VZEAllPhiEta","VZEall;Phi;Eta",20,0,TMath::TwoPi(),40,-4.0,+6.0);  tList->Add( tH2D );
 423   tH1D = new TH1D("TPCPSI","TPCPSI;PSI",72,0,TMath::Pi()); tList->Add( tH1D );
 424   tH1D = new TH1D("TPCPSIA","TPCPSIA;PSIA",72,0,TMath::Pi()); tList->Add( tH1D );
 425   tH1D = new TH1D("TPCPSIB","TPCPSIB;PSIB",72,0,TMath::Pi()); tList->Add( tH1D );
 426   tH1D = new TH1D("VZEPSI","VZEPSI;PSI",72,0,TMath::Pi()); tList->Add( tH1D );
 427   tH1D = new TH1D("VZEPSIA","VZEPSIA;PSIA",72,0,TMath::Pi()); tList->Add( tH1D );
 428   tH1D = new TH1D("VZEPSIB","VZEPSIB;PSIB",72,0,TMath::Pi()); tList->Add( tH1D );
 429   tPF1 = new TProfile("TPCQ","TPCQ",6,0.5,6.5,"s");
 430   tPF1->GetXaxis()->SetBinLabel(1,"Qay"); tPF1->GetXaxis()->SetBinLabel(2,"Qax");
 431   tPF1->GetXaxis()->SetBinLabel(3,"Qby"); tPF1->GetXaxis()->SetBinLabel(4,"Qbx");
 432   tPF1->GetXaxis()->SetBinLabel(5,"Qy");  tPF1->GetXaxis()->SetBinLabel(6,"Qx");
 433   tList->Add( tPF1 );
 434   tPF1 = new TProfile("VZEQ","VZEQ",6,0.5,6.5,"s");
 435   tPF1->GetXaxis()->SetBinLabel(1,"Qay"); tPF1->GetXaxis()->SetBinLabel(2,"Qax");
 436   tPF1->GetXaxis()->SetBinLabel(3,"Qby"); tPF1->GetXaxis()->SetBinLabel(4,"Qbx");
 437   tPF1->GetXaxis()->SetBinLabel(5,"Qy");  tPF1->GetXaxis()->SetBinLabel(6,"Qx");
 438   tList->Add( tPF1 );
 439
 440   fList->Add(tList);
 441   if(!fSkipFlow) {
 442     tList=new TList(); tList->SetName("TPCRFPall"); tList->SetOwner(); AddTPCRFPSpy(tList); fList->Add(tList);
 443     tList=new TList(); tList->SetName("TPCRFPsel"); tList->SetOwner(); AddTPCRFPSpy(tList); fList->Add(tList);
 444   }
 445 }
 446 //=======================================================================
 447 void AliAnalysisTaskFlowStrange::AddQACandidates() {
 448   // function to add histogramming for candidates
 449   if(fSkipSelection) return;
 450
 451   TList *tList;
 452   TH1D *tH1D;
 453   TH2D *tH2D;
 454   TH3D *tH3D;
 455
 456   //reconstruction
 457   if(fReadESD) {
 458     tList=new TList(); tList->SetName("TrkAll"); tList->SetOwner(); AddTracksSpy(tList); fList->Add(tList);
 459     tList=new TList(); tList->SetName("TrkSel"); tList->SetOwner(); AddTracksSpy(tList); fList->Add(tList);
 460     tH2D = new TH2D("NPAIR", "NPAIR;NPOS;NNEG",1000,0,5000,1000,0,5000); tList->Add(tH2D);
 461     tH2D = new TH2D("PtIPXY","PtIPXY;Pt;IPxy", 100,0,10,200,-10,+10); tList->Add(tH2D);
 462   }
 463   //candidates
 464   tList=new TList(); tList->SetName("RecAll"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
 465   tH2D = new TH2D("V0SADC","V0S AFTER DAUGHTER CUTS;V0ALL;V0IMW",100,0,1000,100,0,1000); tList->Add(tH2D);
 466   tList=new TList(); tList->SetName("RecSel"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
 467   //daughters
 468   tList=new TList(); tList->SetName("TrkDau"); tList->SetOwner(); AddTracksSpy(tList); fList->Add(tList);
 469   if(!fSkipDHcorr) {
 470     //corr
 471     tList=new TList(); tList->SetName("DHCORR"); tList->SetOwner();
 472     tH3D = new TH3D("DPHI","DPHI;dPT;dPHI;dETA", 20, -1, +1, 120, -TMath::TwoPi(), TMath::TwoPi(), 16, -1.6, +1.6 ); tList->Add(tH3D);
 473     fList->Add(tList);
 474   }
 475   if(fQAlevel>1) {
 476     // IN-OUT
 477     tList=new TList(); tList->SetName("RecAllIP"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
 478     tList=new TList(); tList->SetName("RecAllOP"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
 479     tList=new TList(); tList->SetName("RecSelIP"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
 480     tList=new TList(); tList->SetName("RecSelOP"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
 481   }
 482   //match
 483   if(fReadMC) {
 484     tList=new TList(); tList->SetName("RecMth"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
 485     tH1D = new TH1D("MCOrigin", "MCOrigin;Rad2",1000,0,100); tList->Add(tH1D);
 486     tH2D = new TH2D("PTRes", "PTRes;MCPt;DAT-MC/MC",  100,   0,  20,100,-0.2,+0.2); tList->Add(tH2D);
 487     tH2D = new TH2D("ETARes","ETARes;MCETA;DAT-MC/MC", 16,   0, 0.8,100,-0.5,+0.5); tList->Add(tH2D);
 488     tH2D = new TH2D("RXYRes","RXYRes;MCRXY;DAT-MC/MC",100,   0,  20,100,-1,1);      tList->Add(tH2D);
 489     tH2D = new TH2D("DLERes","DLERes;MCDLE;DAT-MC/MC",100,   0,  20,100,-1,1);      tList->Add(tH2D);
 490     tH2D = new TH2D("RAPRes","RAPRes;MCRAP;DAT-MC/MC", 10,   0, 0.5,100,-0.5,+0.5); tList->Add(tH2D);
 491     tH2D = new TH2D("APARes","APARes;MCAPA;DAT-MC/MC", 24,-1.2, 1.2,100,-0.5,+0.5); tList->Add(tH2D);
 492     tH2D = new TH2D("APQRes","APQRes;MCAPQ;DAT-MC/MC", 25,   0,0.25,100,-0.3,+0.3); tList->Add(tH2D);
 493
 494     tList=new TList(); tList->SetName("TrkMth"); tList->SetOwner(); AddTracksSpy(tList); fList->Add(tList);
 495     tList=new TList(); tList->SetName("MthFDW"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
 496     tH1D = new TH1D("MCOrigin", "MCOrigin;Rad2",1000,0,100); tList->Add(tH1D);
 497   }
 498   //stack
 499   if(fReadMC) {
 500     tList=new TList(); tList->SetName("GenTru"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
 501     tList=new TList(); tList->SetName("MCTK0sGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
 502     tList=new TList(); tList->SetName("MCTLdaGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
 503     tList=new TList(); tList->SetName("MCTPhiGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
 504     tList=new TList(); tList->SetName("MCTXiGenAcc");  tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
 505     tList=new TList(); tList->SetName("MCTK0s"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
 506     tList=new TList(); tList->SetName("MCTLda"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
 507   }
 508 }
 509 //=======================================================================
 510 void AliAnalysisTaskFlowStrange::Exec(Option_t* option) {
 511   // bypassing ::exec (needed because of AMPT)
 512   ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(0);
 513   if(fAvoidExec) {
 514     AliAnalysisTaskFlowStrange::UserExec(option);
 515   } else {
 516     AliAnalysisTaskSE::Exec(option);
 517   }
 518 }
 519 //=======================================================================
 520 void AliAnalysisTaskFlowStrange::UserExec(Option_t *option) {
 521   // bridge
 522   ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(1);
 523   AliAnalysisTaskFlowStrange::MyUserExec(option);
 524 }
 525 //=======================================================================
 526 void AliAnalysisTaskFlowStrange::MyNotifyRun() {
 527   if(fQAlevel>5 && !fReadESD) AddVZEQA();
 528   if(fVZEsave) AddVZEROResponse();
 529 }
 530 //=======================================================================
 531 Bool_t AliAnalysisTaskFlowStrange::CalibrateEvent() {
 532   if(fVZEsave) SaveVZEROResponse();
 533   if(fQAlevel>5 && !fReadESD) SaveVZEROQA(); // 2BIMPROVED
 534   Bool_t okay=kTRUE;
 535   if(fVZEload) {
 536     LoadVZEROResponse();
 537     if(!fVZEResponse) okay = kFALSE;
 538   }
 539   return okay;
 540 }
 541 //=======================================================================
 542 Bool_t AliAnalysisTaskFlowStrange::AcceptAAEvent(AliESDEvent *tESD) {
 543   Double_t acceptEvent=kTRUE;
 544   Double_t tTPCVtxZ = tESD->GetPrimaryVertexTPC()->GetZ();
 545   if(tESD->GetPrimaryVertexTPC()->GetNContributors()<=0) return kFALSE;
 546   Double_t tSPDVtxZ = tESD->GetPrimaryVertexSPD()->GetZ();
 547   if(tESD->GetPrimaryVertexSPD()->GetNContributors()<=0) return kFALSE;
 548   // EventCuts
 549   AliCentrality *cent = tESD->GetCentrality();
 550   Double_t cc1, cc2;
 551   cc1 = cent->GetCentralityPercentile("V0M");
 552   cc2 = cent->GetCentralityPercentile("TRK");
 553   TString mycent = fCentMethod;
 554   if(fCentMethod.Contains("V0MTRK")) {
 555     acceptEvent = TMath::Abs(cc1-cc2)>5.0?kFALSE:acceptEvent; // a la Alex
 556     mycent = "V0M";
 557   }
 558   fThisCent = cent->GetCentralityPercentile( mycent );
 559   acceptEvent = (fThisCent<fCentPerMin||fThisCent>fCentPerMax)?kFALSE:acceptEvent;
 560   acceptEvent = TMath::Abs(tTPCVtxZ-tSPDVtxZ)>0.5?kFALSE:acceptEvent;
 561   acceptEvent = TMath::Abs(tTPCVtxZ)>10.0?kFALSE:acceptEvent;
 562   if(fQAlevel>0) ((TH2D*)((TList*)fList->FindObject("EventSpy"))->FindObject("VTXZ"))->Fill( tTPCVtxZ, tSPDVtxZ );
 563   if(fQAlevel>0) ((TH2D*)((TList*)fList->FindObject("EventSpy"))->FindObject("CCCC"))->Fill( cc1, cc2 );
 564   // EndOfCuts
 565   return acceptEvent;
 566 }
 567 //=======================================================================
 568 Bool_t AliAnalysisTaskFlowStrange::AcceptAAEvent(AliAODEvent *tAOD) {
 569   Double_t acceptEvent=kTRUE;
 570   //=>Pile-up rejection (hardcoded)
 571   Double_t tVtxZ = tAOD->GetPrimaryVertex()->GetZ();
 572   if(tAOD->GetPrimaryVertex()->GetNContributors()<=0) return kFALSE;
 573   Double_t tSPDVtxZ = tAOD->GetPrimaryVertexSPD()->GetZ();
 574   if(tAOD->GetPrimaryVertexSPD()->GetNContributors()<=0) return kFALSE;
 575   Int_t tpc = RefMultTPC();
 576   Int_t glo = RefMultGlobal();
 577   if(fExtraEventRejection) {
 578     TString name = tAOD->GetPrimaryVertex()->GetTitle();
 579     if( !name.Contains("VertexerTracks") ) return kFALSE;
 580     if( ( Float_t(tpc) < -40.3+1.22*glo ) || ( Float_t(tpc)>(32.1+1.59*glo) ) ) return kFALSE;
 581   }
 582   // EventCuts
 583   AliCentrality *cent = tAOD->GetHeader()->GetCentralityP();
 584   Double_t cc1, cc2;
 585   cc1 = cent->GetCentralityPercentile("V0M");
 586   cc2 = cent->GetCentralityPercentile("TRK");
 587   TString mycent = fCentMethod;
 588   if(fCentMethod.Contains("V0MTRK")) {
 589     acceptEvent = TMath::Abs(cc1-cc2)>5.0?kFALSE:acceptEvent;
 590     mycent = "V0M";
 591   }
 592   fThisCent = cent->GetCentralityPercentile( mycent );
 593
 594   Double_t xsec=0;
 595   if(fReadMC) {
 596     AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(tAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 597     if (!mcHeader) {
 598       return kFALSE;
 599     }
 600     xsec = mcHeader->GetCrossSection();
 601     //fMCEP = mcHeader->GetReactionPlaneAngle();
 602     fMCEP = mcHeader->GetReactionPlaneAngle() + (gRandom->Rndm()>0.5)*TMath::Pi();
 603   }
 604
 605   acceptEvent = (fThisCent<fCentPerMin||fThisCent>fCentPerMax)?kFALSE:acceptEvent;
 606   acceptEvent = TMath::Abs(tVtxZ-tSPDVtxZ)>0.5?kFALSE:acceptEvent;
 607   acceptEvent = TMath::Abs(tVtxZ)>10.0?kFALSE:acceptEvent;
 608   if(fQAlevel>0) {
 609     ((TH2D*)((TList*)fList->FindObject("EventSpy"))->FindObject("VTXZ"))->Fill( tVtxZ, tSPDVtxZ );
 610     ((TH2D*)((TList*)fList->FindObject("EventSpy"))->FindObject("CCCC"))->Fill( cc1, cc2 );
 611     ((TH2D*)((TList*)fList->FindObject("EventSpy"))->FindObject("REFM"))->Fill( tpc, glo );
 612     if(fReadMC) {
 613       ((TH1D*)((TList*)fList->FindObject("EventSpy"))->FindObject("MCCC"))->Fill( xsec );
 614       ((TH1D*)((TList*)fList->FindObject("EventSpy"))->FindObject("MCEP"))->Fill( fMCEP );
 615     }
 616   }
 617   // EndOfCuts
 618   return acceptEvent;
 619 }
 620 //=======================================================================
 621 Bool_t AliAnalysisTaskFlowStrange::AcceptPPEvent(AliAODEvent *tAOD) {
 622   Double_t acceptEvent=kTRUE;
 623   Double_t tVtxZ = tAOD->GetPrimaryVertex()->GetZ();
 624   if(tAOD->GetPrimaryVertex()->GetNContributors()<=0) return kFALSE;
 625   Double_t tSPDVtxZ = tAOD->GetPrimaryVertexSPD()->GetZ();
 626   if(tAOD->GetPrimaryVertexSPD()->GetNContributors()<=0) return kFALSE;
 627   // EventCuts
 628   AliCentrality *cent = tAOD->GetHeader()->GetCentralityP();
 629   Double_t cc1, cc2;
 630   cc1 = cent->GetCentralityPercentile("V0M");
 631   cc2 = cent->GetCentralityPercentile("TRK");
 632   fThisCent = GetReferenceMultiplicity();
 633   //for pp i use fCentPerXXX to select on multiplicity
 634   acceptEvent = (fThisCent<fCentPerMin||fThisCent>fCentPerMax)?kFALSE:acceptEvent;
 635   acceptEvent = TMath::Abs(tVtxZ-tSPDVtxZ)>0.5?kFALSE:acceptEvent;
 636   acceptEvent = TMath::Abs(tVtxZ)>10.0?kFALSE:acceptEvent;
 637   if(fQAlevel>0) ((TH2D*)((TList*)fList->FindObject("EventSpy"))->FindObject("VTXZ"))->Fill( tVtxZ, tSPDVtxZ );
 638   if(fQAlevel>0) ((TH2D*)((TList*)fList->FindObject("EventSpy"))->FindObject("CCCC"))->Fill( cc1, cc2 );
 639   // EndOfCuts
 640   return acceptEvent;
 641 }
 642 //=======================================================================
 643 Int_t AliAnalysisTaskFlowStrange::GetReferenceMultiplicity() { //toberefined
 644   AliAODEvent *tAOD = (AliAODEvent *) InputEvent();
 645   if(!tAOD) return -1;
 646   AliAODTrack *track;
 647   Int_t rawN = tAOD->GetNumberOfTracks();
 648   Int_t ref=0;
 649   for(Int_t id=0; id!=rawN; ++id) {
 650     track = tAOD->GetTrack(id);
 651     if(!track->TestFilterBit(fRFPFilterBit)) continue;
 652     ++ref;
 653   }
 654   return ref;
 655 }
 656 //=======================================================================
 657 Bool_t AliAnalysisTaskFlowStrange::AcceptPAEvent(AliAODEvent *tAOD) {
 658   //if(aod->GetHeader()->GetEventNumberESDFile() == 0) return; //rejecting first chunk NOT NEEDED ANYMORE
 659   Int_t bc2 = tAOD->GetHeader()->GetIRInt2ClosestInteractionMap();
 660   if(bc2!=0) return kFALSE;
 661   Int_t bc1 = tAOD->GetHeader()->GetIRInt1ClosestInteractionMap();
 662   if(bc1!=0) return kFALSE;
 663   Short_t isPileup = tAOD->IsPileupFromSPD(5);
 664   if(isPileup!=0) return kFALSE;
 665   if(tAOD->GetHeader()->GetRefMultiplicityComb08()<0) return kFALSE;
 666
 667   const AliAODVertex* spdVtx = tAOD->GetPrimaryVertexSPD();
 668   if(!spdVtx) return kFALSE;
 669   if(spdVtx->GetNContributors()<=0) return kFALSE;
 670
 671   const AliAODVertex* tpcVtx=NULL;
 672   Int_t nVertices = tAOD->GetNumberOfVertices();
 673   for(Int_t iVertices = 0; iVertices < nVertices; iVertices++){
 674     const AliAODVertex* vertex = tAOD->GetVertex(iVertices);
 675     if (vertex->GetType() != AliAODVertex::kMainTPC) continue;
 676     tpcVtx = vertex;
 677   }
 678   if(!tpcVtx) return kFALSE;
 679   if(tpcVtx->GetNContributors()<=0) return kFALSE;
 680   Double_t tTPCVtxZ = tpcVtx->GetZ();
 681   Double_t tSPDVtxZ = spdVtx->GetZ();
 682   if (TMath::Abs(tSPDVtxZ - tTPCVtxZ)>2.0) return kFALSE;
 683   if(plpMV(tAOD)) return kFALSE;
 684
 685   Double_t acceptEvent=kTRUE;
 686   // EventCuts
 687   AliCentrality *cent = tAOD->GetHeader()->GetCentralityP();
 688   Double_t cc1, cc2;
 689   cc1 = cent->GetCentralityPercentile("V0M");
 690   cc2 = cent->GetCentralityPercentile("TRK");
 691   if(fCentMethod.Contains("V0MTRK")) fCentMethod = "V0M";
 692   fThisCent = cent->GetCentralityPercentile( fCentMethod );
 693   acceptEvent = (fThisCent<fCentPerMin||fThisCent>fCentPerMax)?kFALSE:acceptEvent;
 694   acceptEvent = TMath::Abs(tTPCVtxZ)>10.0?kFALSE:acceptEvent;
 695   // EndOfCuts
 696   if(fQAlevel>0) ((TH2D*)((TList*)fList->FindObject("EventSpy"))->FindObject("VTXZ"))->Fill( tTPCVtxZ, tSPDVtxZ );
 697   if(fQAlevel>0) ((TH2D*)((TList*)fList->FindObject("EventSpy"))->FindObject("CCCC"))->Fill( cc1, cc2 );
 698   return acceptEvent;
 699 }
 700 //=======================================================================
 701 void AliAnalysisTaskFlowStrange::MyUserExec(Option_t *) {
 702   // MAIN ROUTINE
 703   TStopwatch tTime;
 704   tTime.Start();
 705   fCandidates->SetLast(-1);
 706   AliESDEvent *tESD=dynamic_cast<AliESDEvent*>(InputEvent());
 707   AliAODEvent *tAOD=dynamic_cast<AliAODEvent*>(InputEvent());
 708   Int_t thisRun = fRunNumber;
 709   //=>check event
 710   Bool_t acceptEvent=kFALSE;
 711   if(fReadESD) {
 712     if(!tESD) {ResetContainers(); Publish(); return;}
 713     acceptEvent = fRunOnpp?kFALSE:fRunOnpA?kFALSE:AcceptAAEvent(tESD);
 714     thisRun = tESD->GetRunNumber();
 715   } else {
 716     if(!tAOD) {ResetContainers(); Publish(); return;}
 717     acceptEvent = fRunOnpp?AcceptPPEvent(tAOD):fRunOnpA?AcceptPAEvent(tAOD):AcceptAAEvent(tAOD);
 718     thisRun = tAOD->GetRunNumber();
 719   }
 720   if(thisRun!=fRunNumber) {
 721     fRunNumber = thisRun;
 722     MyNotifyRun();
 723   }
 724   ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(2);
 725   //=>does the event clear?
 726   if(!acceptEvent) {ResetContainers(); Publish(); return;}
 727   // healthy event incomming
 728   if( !CalibrateEvent() ) { // saves/retrieves/qas VZEROCAL
 729     ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(5);
 730     ResetContainers(); Publish(); return; // issue retrieving callibration
 731   }
 732   if(!fSkipFlow) {
 733     MakeQVectors();
 734     if(fPsi2<-0.1) {
 735       ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(4);
 736       ResetContainers(); Publish(); return;
 737     }
 738   }
 739   //=>great, lets do our stuff!
 740   ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(3);
 741   //=>load candidates
 742   if(!fSkipSelection) {
 743     if(fReadESD) {
 744       ReadFromESD(tESD);
 745     } else {
 746       if(fSpecie<10) ReadFromAODv0(tAOD);
 747       else ChargeParticles(tAOD);
 748     }
 749     if(fUseFP) {
 750       if(!fSkipDHcorr) MakeDHcorr();
 751       AddCandidates();
 752     }
 753     //=>flow
 754     //=>done
 755   }
 756   tTime.Stop();
 757   ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("RealTime"))->Fill( TMath::Log( tTime.RealTime() ) );
 758   Publish();
 759 }
 760 //=======================================================================
 761 void AliAnalysisTaskFlowStrange::Publish() {
 762   PostData(1,fList);
 763   if(fUseFP) {
 764     PostData(2,fTPCevent);
 765     PostData(3,fVZEevent);
 766   }
 767 }
 768 //=======================================================================
 769 void AliAnalysisTaskFlowStrange::ReadFromESD(AliESDEvent *tESD) {
 770   AliStack *stack=NULL;
 771   if(fReadMC) {
 772     AliMCEvent *mcevent=NULL;
 773     mcevent = MCEvent();
 774     if(mcevent) stack = mcevent->Stack();
 775   }
 776
 777   Int_t num = tESD->GetNumberOfTracks();
 778   AliESDtrack *myTrack;
 779   Int_t plist[3000], nlist[3000], np=0, nn=0;
 780   Double_t pd0[3000], nd0[3000];
 781   for (Int_t i=0; i!=num; ++i) {
 782     myTrack = (AliESDtrack*) tESD->GetTrack(i);
 783     if(!myTrack) continue;
 784     LoadTrack(myTrack);
 785     FillTrackSpy("TrkAll");
 786     if(!AcceptDaughter()) continue;
 787     FillTrackSpy("TrkSel");
 788     ((TH2D*)((TList*)fList->FindObject("TrkSel"))->FindObject("PtIPXY" ))->Fill( myTrack->Pt(), fDaughterImpactParameterXY );
 789     if( myTrack->Charge()>0 ) {
 790       pd0[np] = fDaughterImpactParameterXY;
 791       plist[np++] = i;
 792     } else {
 793       nd0[nn] = fDaughterImpactParameterXY;
 794       nlist[nn++] = i;
 795     }
 796   }
 797   ((TH1D*)((TList*)fList->FindObject("TrkSel"))->FindObject("NPAIR" ))->Fill( np,nn );
 798   const AliESDVertex *vtx = tESD->GetPrimaryVertex();
 799   AliESDtrack *pT, *nT;
 800   for(int p=0; p!=np; ++p) {
 801     pT = (AliESDtrack*) tESD->GetTrack( plist[p] );
 802     for(int n=0; n!=nn; ++n) {
 803       nT = (AliESDtrack*) tESD->GetTrack( nlist[n] );
 804       fDecayProductIPXY = pd0[p]*nd0[n];
 805       AliExternalTrackParam pETP(*pT), nETP(*nT);
 806       Double_t xa, xb;
 807       pETP.GetDCA(&nETP,tESD->GetMagneticField(),xa,xb);
 808       fDecayDCAdaughters = pETP.PropagateToDCA(&nETP,tESD->GetMagneticField());
 809       AliESDv0 vertex( nETP,nlist[n], pETP,plist[p] );
 810       fDecayCosinePointingAngleXY = CosThetaPointXY( &vertex, vtx );
 811       fDecayRadXY = DecayLengthXY( &vertex, vtx );
 812       fDecayPt = vertex.Pt();
 813       fDecayPhi = vertex.Phi();
 814       fDecayEta = vertex.Eta();
 815       Double_t pmx, pmy, pmz, nmx, nmy, nmz;
 816       vertex.GetNPxPyPz(nmx,nmy,nmz);
 817       vertex.GetPPxPyPz(pmx,pmy,pmz);
 818       TVector3 mom1(pmx,pmy,pmz), mom2(nmx,nmy,nmz), mom(vertex.Px(),vertex.Py(),vertex.Pz());
 819       Double_t qlpos = mom1.Dot(mom)/mom.Mag();
 820       Double_t qlneg = mom2.Dot(mom)/mom.Mag();
 821       fDecayQt = mom1.Perp(mom);
 822       fDecayAlpha = (qlpos-qlneg)/(qlpos+qlneg);
 823       Double_t mpi = 0.13957018;
 824       if(fSpecie==0) {
 825         Double_t eppi = TMath::Sqrt( mpi*mpi + pmx*pmx + pmy*pmy + pmz*pmz );
 826         Double_t enpi = TMath::Sqrt( mpi*mpi + nmx*nmx + nmy*nmy + nmz*nmz );
 827         fDecayMass = TMath::Sqrt( mpi*mpi + mpi*mpi + 2*(eppi*enpi - pmx*nmx - pmy*nmy - pmz*nmz ) );
 828         fDecayRapidity = vertex.RapK0Short();
 829       } else {
 830         Double_t mpr = 0.938272013;
 831         Double_t epi, epr;
 832         if(fDecayAlpha>0) {
 833           epr = TMath::Sqrt( mpr*mpr + pmx*pmx + pmy*pmy + pmz*pmz );
 834           epi = TMath::Sqrt( mpi*mpi + nmx*nmx + nmy*nmy + nmz*nmz );
 835         } else {
 836           epi = TMath::Sqrt( mpi*mpi + pmx*pmx + pmy*pmy + pmz*pmz );
 837           epr = TMath::Sqrt( mpr*mpr + nmx*nmx + nmy*nmy + nmz*nmz );
 838         }
 839         fDecayMass = TMath::Sqrt( mpi*mpi + mpr*mpr + 2*(epi*epr - pmx*nmx - pmy*nmy - pmz*nmz ) );
 840         fDecayRapidity = vertex.RapLambda();
 841       }
 842       Double_t energy = TMath::Sqrt( fDecayMass*fDecayMass + vertex.Px()*vertex.Px() + vertex.Py()*vertex.Py() + vertex.Pz()*vertex.Pz() );
 843       Double_t gamma = energy/fDecayMass;
 844       fDecayDecayLength = DecayLength( &vertex, vtx )/gamma;
 845       Double_t dPHI = fDecayPhi;
 846       Double_t dDPHI = dPHI - fPsi2;
 847       if( dDPHI < 0 ) dDPHI += TMath::TwoPi();
 848       if( dDPHI > TMath::Pi() ) dDPHI = TMath::TwoPi()-dDPHI;
 849       if(fQAlevel>1) {
 850         if( (dDPHI>TMath::PiOver4()) && (dDPHI<3*TMath::PiOver4()) ) FillCandidateSpy("RecAllOP");
 851         else FillCandidateSpy("RecAllIP");
 852       }
 853       FillCandidateSpy("RecAll");
 854       ((TH2D*)((TList*)fList->FindObject("RecAll"))->FindObject("D0PD0N"))->Fill( pd0[p],nd0[n] );
 855       ((TH2D*)((TList*)fList->FindObject("RecAll"))->FindObject("XPOSXNEG"))->Fill( xa, xb );
 856       if(!AcceptCandidate()) continue;
 857       if(fDecayMass<fMinMass) continue;
 858       if(fDecayMass>fMaxMass) continue;
 859       // PID missing
 860       if(fQAlevel>1) {
 861         if( (dDPHI>TMath::PiOver4()) && (dDPHI<3*TMath::PiOver4()) ) FillCandidateSpy("RecSelOP");
 862         else FillCandidateSpy("RecSelIP");
 863       }
 864       FillCandidateSpy("RecSel");
 865       ((TH2D*)((TList*)fList->FindObject("RecSel"))->FindObject("D0PD0N"))->Fill( pd0[p],nd0[n] );
 866       ((TH2D*)((TList*)fList->FindObject("RecSel"))->FindObject("XPOSXNEG"))->Fill( xa, xb );
 867
 868       fDecayIDneg = nT->GetID();
 869       fDecayIDpos = pT->GetID();
 870       MakeTrack();
 871       LoadTrack(pT); FillTrackSpy("TrkDau");
 872       LoadTrack(nT); FillTrackSpy("TrkDau");
 873
 874       //===== BEGIN OF MCMATCH
 875       if(stack) {
 876         bool matched = false;
 877         Int_t labelpos = pT->GetLabel();
 878         Int_t labelneg = nT->GetLabel();
 879         Double_t rOri=-1;
 880         if( labelpos>0 && labelneg>0 ) {
 881           TParticle *mcpos = stack->Particle( labelpos );
 882           TParticle *mcneg = stack->Particle( labelneg );
 883           Int_t pdgRecPos = mcpos->GetPdgCode();
 884           Int_t pdgRecNeg = mcneg->GetPdgCode();
 885           if( pdgRecPos==211&&pdgRecNeg==-211 ) if(mcpos->GetMother(0)>0) {
 886             if( mcpos->GetMother(0)==mcneg->GetMother(0) ) {
 887               TParticle *mcmot = stack->Particle( mcpos->GetMother(0) );
 888               rOri = TMath::Sqrt( mcmot->Vx()*mcmot->Vx() + mcmot->Vy()*mcmot->Vy() );
 889               if( TMath::Abs(mcmot->GetPdgCode())==310) {
 890                 if(mcmot->GetNDaughters()==2) matched=true;
 891               }
 892             }
 893           }
 894         }
 895         if(matched) {
 896           FillCandidateSpy("RecMth");
 897           ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("D0PD0N"))->Fill( pd0[p],nd0[n] );
 898           ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("XPOSXNEG"))->Fill( xa, xb );
 899           ((TH1D*)((TList*)fList->FindObject("RecMth"))->FindObject("MCOrigin"))->Fill( rOri );
 900           LoadTrack(pT); FillTrackSpy("TrkMth");
 901           LoadTrack(nT); FillTrackSpy("TrkMth");
 902         }
 903       }
 904       //===== END OF MCMATCH
 905     }
 906   }
 907 }
 908 //=======================================================================
 909 void AliAnalysisTaskFlowStrange::ReadStack(TClonesArray* mcArray) {
 910   if(!mcArray) return;
 911   AliAODMCParticle *myMCTrack, *iMCDau, *jMCDau;
 912   for(int i=0; i!=mcArray->GetEntriesFast(); ++i) {
 913     myMCTrack = dynamic_cast<AliAODMCParticle*>(mcArray->At( i ));
 914     if(!myMCTrack) continue;
 915     int tPDG=310;
 916     if(fSpecie>0) tPDG = 3122;
 917     if( TMath::Abs(myMCTrack->PdgCode())==tPDG )
 918       if( myMCTrack->GetNDaughters() == 2 ) {
 919         Int_t iDau = myMCTrack->GetDaughter(0);
 920         Int_t jDau = myMCTrack->GetDaughter(1);
 921         AliAODMCParticle *posDau=NULL;
 922         AliAODMCParticle *negDau=NULL;
 923         if(iDau>0&&jDau>0) {
 924           iMCDau = dynamic_cast<AliAODMCParticle*>(mcArray->At( iDau ));
 925           jMCDau = dynamic_cast<AliAODMCParticle*>(mcArray->At( jDau ));
 926           if(iMCDau) {
 927             if(iMCDau->Charge()>0) posDau=iMCDau;
 928             else negDau=iMCDau;
 929           }
 930           if(jMCDau) {
 931             if(jMCDau->Charge()>0) posDau=jMCDau;
 932             else negDau=jMCDau;
 933           }
 934         } //got two daughters
 935         if(posDau&&negDau) {
 936           Double_t dx = myMCTrack->Xv() - posDau->Xv();
 937           Double_t dy = myMCTrack->Yv() - posDau->Yv();
 938           Double_t dz = myMCTrack->Zv() - posDau->Zv();
 939           fDecayRadXY = TMath::Sqrt( dx*dx + dy*dy );
 940           TVector3 momPos(posDau->Px(),posDau->Py(),posDau->Pz());
 941           TVector3 momNeg(negDau->Px(),negDau->Py(),negDau->Pz());
 942           TVector3 momTot(myMCTrack->Px(),myMCTrack->Py(),myMCTrack->Pz());
 943           Double_t qlpos = momPos.Dot(momTot)/momTot.Mag();
 944           Double_t qlneg = momNeg.Dot(momTot)/momTot.Mag();
 945           fDecayQt = momPos.Perp(momTot);
 946           fDecayAlpha = 1.-2./(1.+qlpos/qlneg);
 947           fDecayMass = myMCTrack->GetCalcMass();
 948           Double_t energy = myMCTrack->E();
 949           Double_t gamma = energy/fDecayMass;
 950           fDecayDecayLength = TMath::Sqrt(dx*dx+dy*dy+dz*dz)/gamma;
 951           fDecayPt = myMCTrack->Pt();
 952           fDecayPhi = myMCTrack->Phi();
 953           fDecayEta = myMCTrack->Eta();
 954           fDecayRapidity = myMCTrack->Y();
 955           fDecayDCAdaughters = 0;
 956           fDecayCosinePointingAngleXY = 1;
 957           fDecayProductIPXY = -1;
 958           if(AcceptCandidate()) FillCandidateSpy("GenTru");
 959         }
 960       } // k0/lda with two daughters
 961     //==== BEGIN TRACK CUTS
 962     if(myMCTrack->Eta()<-0.8) continue;
 963     if(myMCTrack->Eta()>+0.8) continue;
 964     //==== END TRACK CUTS
 965     switch( TMath::Abs(myMCTrack->PdgCode()) ) {
 966     case (310): //k0s
 967       FillMCParticleSpy( "MCTK0s", myMCTrack );
 968       if( myMCTrack->IsPrimary() )
 969         FillMCParticleSpy( "MCTK0sGenAcc", myMCTrack );
 970       break;
 971     case (3122): //lda
 972       FillMCParticleSpy( "MCTLda", myMCTrack );
 973       if( myMCTrack->IsPrimary() )
 974         FillMCParticleSpy( "MCTLdaGenAcc", myMCTrack );
 975       break;
 976     case (333): //phi
 977       if( myMCTrack->IsPrimary() )
 978         FillMCParticleSpy( "MCTPhiGenAcc", myMCTrack );
 979       break;
 980     case (3312): //xi
 981       if( myMCTrack->IsPrimary() )
 982         FillMCParticleSpy( "MCTXiGenAcc", myMCTrack );
 983       break;
 984     }
 985
 986   }
 987 }
 988 //=======================================================================
 989 Double_t AliAnalysisTaskFlowStrange::CosThetaPointXY(AliESDv0 *me, const AliVVertex *vtx) {
 990   TVector3 mom( me->Px(), me->Py(), 0 );
 991   TVector3 fli( me->Xv()-vtx->GetX(), me->Yv()-vtx->GetY(), 0 );
 992   Double_t ctp = mom.Dot(fli) / mom.Mag() / fli.Mag();
 993   return ctp;
 994 }
 995 //=======================================================================
 996 Double_t AliAnalysisTaskFlowStrange::CosThetaPointXY(AliAODv0 *me, const AliVVertex *vtx) {
 997   TVector3 mom( me->Px(), me->Py(), 0 );
 998   TVector3 fli( me->Xv()-vtx->GetX(), me->Yv()-vtx->GetY(), 0 );
 999   Double_t ctp = mom.Dot(fli) / mom.Mag() / fli.Mag();
1000   return ctp;
1001 }
1002 //=======================================================================
1003 Double_t AliAnalysisTaskFlowStrange::DecayLengthXY(AliESDv0 *me, const AliVVertex *vtx) {
1004   Double_t dx = me->Xv()-vtx->GetX();
1005   Double_t dy = me->Yv()-vtx->GetY();
1006   Double_t dxy = TMath::Sqrt( dx*dx + dy*dy );
1007   return dxy;
1008 }
1009 //=======================================================================
1010 Double_t AliAnalysisTaskFlowStrange::DecayLengthXY(AliAODv0 *me, const AliVVertex *vtx) {
1011   Double_t dx = me->Xv()-vtx->GetX();
1012   Double_t dy = me->Yv()-vtx->GetY();
1013   Double_t dxy = TMath::Sqrt( dx*dx + dy*dy );
1014   return dxy;
1015 }
1016 //=======================================================================
1017 Double_t AliAnalysisTaskFlowStrange::DecayLength(AliESDv0 *me, const AliVVertex *vtx) {
1018   Double_t dx = me->Xv()-vtx->GetX();
1019   Double_t dy = me->Yv()-vtx->GetY();
1020   Double_t dz = me->Zv()-vtx->GetZ();
1021   Double_t dxy = TMath::Sqrt( dx*dx + dy*dy + dz*dz );
1022   return dxy;
1023 }
1024 //=======================================================================
1025 Double_t AliAnalysisTaskFlowStrange::DecayLength(AliAODv0 *me, const AliVVertex *vtx) {
1026   Double_t dx = me->Xv()-vtx->GetX();
1027   Double_t dy = me->Yv()-vtx->GetY();
1028   Double_t dz = me->Zv()-vtx->GetZ();
1029   Double_t dxy = TMath::Sqrt( dx*dx + dy*dy + dz*dz );
1030   return dxy;
1031 }
1032 //=======================================================================
1033 void AliAnalysisTaskFlowStrange::ReadFromAODv0(AliAODEvent *tAOD) {
1034   TClonesArray* mcArray=NULL;
1035   if(fReadMC) {
1036     mcArray = dynamic_cast<TClonesArray*>(tAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1037     ReadStack(mcArray);
1038   }
1039
1040   Int_t nV0s = tAOD->GetNumberOfV0s();
1041   AliAODv0 *myV0;
1042   Int_t v0all=0, v0imw=0;
1043   for (Int_t i=0; i!=nV0s; ++i) {
1044     myV0 = (AliAODv0*) tAOD->GetV0(i);
1045     if(!myV0) continue;
1046     if(!fOnline) if(myV0->GetOnFlyStatus() ) continue;
1047     if(fOnline) if(!myV0->GetOnFlyStatus() ) continue;
1048     AliAODTrack *iT, *jT;
1049     AliAODVertex *vtx = tAOD->GetPrimaryVertex();
1050     Double_t pos[3],cov[6];
1051     vtx->GetXYZ(pos);
1052     vtx->GetCovarianceMatrix(cov);
1053     const AliESDVertex vESD(pos,cov,100.,100);
1054     // TESTING CHARGE
1055     int iPos, iNeg;
1056     iT=(AliAODTrack*) myV0->GetDaughter(0);
1057     if(iT->Charge()>0) {
1058       iPos = 0; iNeg = 1;
1059     } else {
1060       iPos = 1; iNeg = 0;
1061     }
1062     // END OF TEST
1063
1064     iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1065     AliESDtrack ieT( iT );
1066     ieT.SetTPCClusterMap( iT->GetTPCClusterMap() );
1067     ieT.SetTPCSharedMap( iT->GetTPCSharedMap() );
1068     ieT.SetTPCPointsF( iT->GetTPCNclsF() );
1069     ieT.PropagateToDCA(&vESD, tAOD->GetMagneticField(), 100);
1070     LoadTrack(&ieT,iT->Chi2perNDF());
1071     Float_t ip[2];
1072     ieT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1073     fDaughterImpactParameterXY = ip[0];
1074     fDaughterImpactParameterZ = ip[1];
1075     Double_t dD0P = fDaughterImpactParameterXY; //ieT.GetD(pos[0], pos[1], tAOD->GetMagneticField());
1076     if(!AcceptDaughter()) continue;
1077
1078     jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1079     AliESDtrack jeT( jT );
1080     jeT.SetTPCClusterMap( jT->GetTPCClusterMap() );
1081     jeT.SetTPCSharedMap( jT->GetTPCSharedMap() );
1082     jeT.SetTPCPointsF( jT->GetTPCNclsF() );
1083     jeT.PropagateToDCA(&vESD, tAOD->GetMagneticField(), 100);
1084     LoadTrack(&jeT,jT->Chi2perNDF());
1085     jeT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1086     fDaughterImpactParameterXY = ip[0];
1087     fDaughterImpactParameterZ = ip[1];
1088     Double_t dD0N = fDaughterImpactParameterXY; //jeT.GetD(pos[0], pos[1], tAOD->GetMagneticField());
1089     if(!AcceptDaughter()) continue;
1090     if( fExcludeTPCEdges ) {
1091       if( IsAtTPCEdge(iT->Phi(),iT->Pt(),+1,tAOD->GetMagneticField()) ) continue;
1092       if( IsAtTPCEdge(jT->Phi(),jT->Pt(),-1,tAOD->GetMagneticField()) ) continue;
1093     }
1094     Double_t xa, xb;
1095     ieT.GetDCA(&jeT,tAOD->GetMagneticField(),xa,xb);
1096     /*
1097     // cutting out population close to TPC edges :: strange excess saw in 2010
1098     if( fExcludeTPCEdges ) {
1099     Double_t phimod = myV0->Phi();
1100     int sectors[6] = {5,6,9,10,11,12};
1101     for(int ii=0; ii!=6; ++ii)
1102     if( (phimod<(sectors[ii]+1)*TMath::Pi()/9.0) && (phimod>sectors[ii]*TMath::Pi()/9.0) )
1103     return 0;
1104     }
1105     */
1106     if(fSpecie==0)
1107       fDecayRapidity = myV0->RapK0Short();
1108     else
1109       fDecayRapidity = myV0->RapLambda();
1110     fDecayDCAdaughters = myV0->DcaV0Daughters();
1111     fDecayCosinePointingAngleXY = CosThetaPointXY( myV0, vtx );
1112     fDecayRadXY = DecayLengthXY( myV0, vtx );
1113     fDecayPt = myV0->Pt();
1114     fDecayPhi = myV0->Phi();
1115     fDecayEta = myV0->Eta();
1116     fDecayProductIPXY = dD0P*dD0N;
1117     fDecayQt = myV0->PtArmV0();
1118     fDecayAlpha = myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1119     if(myV0->ChargeProng(iPos)<0) fDecayAlpha = -fDecayAlpha; // protects for a change in convention
1120     fDecayPt = myV0->Pt();
1121     fDecayEta = myV0->Eta();
1122     if( fSpecie==0 ) {
1123       fDecayMass = myV0->MassK0Short();
1124     } else {
1125       if(fDecayAlpha>0) fDecayMass = myV0->MassLambda();
1126       else fDecayMass = myV0->MassAntiLambda();
1127     }
1128     v0all++;
1129     if(fDecayMass<fMinMass) continue;
1130     if(fDecayMass>fMaxMass) continue;
1131     v0imw++;
1132     Double_t energy = TMath::Sqrt( fDecayMass*fDecayMass + myV0->Px()*myV0->Px() + myV0->Py()*myV0->Py() + myV0->Pz()*myV0->Pz() );
1133     Double_t gamma = energy/fDecayMass;
1134     fDecayDecayLength = DecayLength( myV0, vtx )/gamma;
1135     Double_t dPHI = fDecayPhi;
1136     Double_t dDPHI = dPHI - fPsi2;
1137     if( dDPHI < 0 ) dDPHI += TMath::TwoPi();
1138     if( dDPHI > TMath::Pi() ) dDPHI = TMath::TwoPi()-dDPHI;
1139     if(fQAlevel>1) {
1140       if( (dDPHI>TMath::PiOver4()) && (dDPHI<3*TMath::PiOver4()) ) FillCandidateSpy("RecAllOP");
1141       else FillCandidateSpy("RecAllIP");
1142     }
1143     FillCandidateSpy("RecAll");
1144     ((TH2D*)((TList*)fList->FindObject("RecAll"))->FindObject("D0PD0N"))->Fill( dD0P, dD0N );
1145     ((TH2D*)((TList*)fList->FindObject("RecAll"))->FindObject("XPOSXNEG"))->Fill( xa, xb );
1146     if(!AcceptCandidate()) continue;
1147     //PID for lambda::proton
1148     if( fSpecie>0 )
1149       if(fDecayPt<1.2) {
1150         if(fDecayAlpha>0) {
1151           if( !PassesPIDCuts(&ieT) ) continue; //positive track
1152         } else {
1153           if( !PassesPIDCuts(&jeT) ) continue; //negative track
1154         }
1155       }
1156     if(fQAlevel>1) {
1157       if( (dDPHI>TMath::PiOver4()) && (dDPHI<3*TMath::PiOver4()) ) FillCandidateSpy("RecSelOP");
1158       else FillCandidateSpy("RecSelIP");
1159     }
1160     FillCandidateSpy("RecSel");
1161     ((TH2D*)((TList*)fList->FindObject("RecSel"))->FindObject("D0PD0N"))->Fill( dD0P, dD0N );
1162     ((TH2D*)((TList*)fList->FindObject("RecSel"))->FindObject("XPOSXNEG"))->Fill( xa, xb );
1163     fDecayIDneg = iT->GetID();
1164     fDecayIDpos = jT->GetID();
1165     MakeTrack();
1166     LoadTrack(&ieT,iT->Chi2perNDF());
1167     ieT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1168     fDaughterImpactParameterXY = ip[0];
1169     fDaughterImpactParameterZ = ip[1];
1170     FillTrackSpy("TrkDau");
1171     LoadTrack(&jeT,jT->Chi2perNDF());
1172     jeT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1173     fDaughterImpactParameterXY = ip[0];
1174     fDaughterImpactParameterZ = ip[1];
1175     FillTrackSpy("TrkDau");
1176     //===== BEGIN OF MCMATCH
1177     if(mcArray) {
1178       bool matched = false;
1179       bool feeddown = false;
1180       Int_t labelpos = iT->GetLabel();
1181       Int_t labelneg = jT->GetLabel();
1182       Double_t rOri=-1;
1183       Double_t mcPt=-1, mcEta=-1, mcRap=-1, mcRxy=-1, mcDle=-1, mcApa=-100, mcApq=-1;
1184       Double_t resPt=-1, resEta=-1, resRap=-1, resRxy=-1, resDle=-1, resApa=-1, resApq=-1;
1185       if( labelpos>0 && labelneg>0 ) {
1186         AliAODMCParticle *mcpos = (AliAODMCParticle*) mcArray->At( labelpos );
1187         AliAODMCParticle *mcneg = (AliAODMCParticle*) mcArray->At( labelneg );
1188         Int_t pdgRecPos = mcpos->GetPdgCode();
1189         Int_t pdgRecNeg = mcneg->GetPdgCode();
1190         int pospdg=211, negpdg=211;
1191         int mompdg=310, fdwpdg=333;
1192         if(fSpecie>0) {
1193           mompdg=3122;
1194           fdwpdg=3312;
1195           if(fDecayAlpha>0) {
1196             pospdg=2212; negpdg=211;
1197           } else {
1198             negpdg=2212; pospdg=211;
1199           }
1200         }
1201         if( TMath::Abs(pdgRecPos)==pospdg&&TMath::Abs(pdgRecNeg)==negpdg )
1202           if(mcpos->GetMother()>0)
1203             if( mcpos->GetMother()==mcneg->GetMother() ) {
1204               AliAODMCParticle *mcmot = (AliAODMCParticle*) mcArray->At( mcpos->GetMother() );
1205               rOri = TMath::Sqrt( mcmot->Xv()*mcmot->Xv() + mcmot->Yv()*mcmot->Yv() );
1206               mcPt = mcmot->Pt();
1207               mcEta = TMath::Abs( mcmot->Eta() );
1208               mcRap = TMath::Abs( mcmot->Y() );
1209               if(!TMath::AreEqualAbs(mcPt,0,1e-6))  resPt = (fDecayPt - mcPt) / mcPt;
1210               if(!TMath::AreEqualAbs(mcEta,0,1e-6)) resEta = (fDecayEta - mcEta) / mcEta;
1211               if(!TMath::AreEqualAbs(mcRap,0,1e-6)) resRap = (fDecayRapidity - mcRap) / mcRap;
1212               if( TMath::Abs(mcmot->GetPdgCode())==mompdg) {
1213                 if(mcmot->GetNDaughters()==2) {
1214                   matched=true;
1215                   Double_t dx = mcmot->Xv() - mcpos->Xv();
1216                   Double_t dy = mcmot->Yv() - mcpos->Yv();
1217                   Double_t dz = mcmot->Zv() - mcpos->Zv();
1218                   Double_t mcGamma = mcmot->E() / mcmot->GetCalcMass();
1219                   mcRxy = TMath::Sqrt( dx*dx + dy*dy );
1220                   mcDle = TMath::Sqrt(dx*dx+dy*dy+dz*dz)/mcGamma;
1221                   if(!TMath::AreEqualAbs(mcRxy,0,1e-6)) resRxy = (fDecayRadXY - mcRxy) / mcRxy;
1222                   if(!TMath::AreEqualAbs(mcDle,0,1e-6)) resDle = (fDecayDecayLength - mcDle) / mcDle;
1223                   TVector3 momPos(mcpos->Px(),mcpos->Py(),mcpos->Pz());
1224                   TVector3 momNeg(mcneg->Px(),mcneg->Py(),mcneg->Pz());
1225                   TVector3 momTot(mcmot->Px(),mcmot->Py(),mcmot->Pz());
1226                   Double_t qlpos = momPos.Dot(momTot)/momTot.Mag();
1227                   Double_t qlneg = momNeg.Dot(momTot)/momTot.Mag();
1228                   mcApq = momPos.Perp(momTot);
1229                   mcApa = 1.-2./(1.+qlpos/qlneg);
1230                   if(!TMath::AreEqualAbs(mcApq,0,1e-6)) resApq = (fDecayQt - mcApq) / mcApq;
1231                   if(!TMath::AreEqualAbs(mcApa,0,1e-6)) resApa = (fDecayAlpha - mcApa) / mcApa;
1232                 }
1233                 if(mcmot->GetMother()>0) {
1234                   AliAODMCParticle *mcfdw = (AliAODMCParticle*) mcArray->At( mcmot->GetMother() );
1235                   if( TMath::Abs(mcfdw->GetPdgCode())==fdwpdg)
1236                     feeddown=true;
1237                 } // k0/lda have mother
1238               } // mother matches k0/lda
1239             } // both have same mother
1240       } // match to MCparticles
1241       if(matched) {
1242         FillCandidateSpy("RecMth");
1243         ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("D0PD0N"))->Fill( dD0P,dD0N );
1244         ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("XPOSXNEG"))->Fill( xa, xb );
1245         ((TH1D*)((TList*)fList->FindObject("RecMth"))->FindObject("MCOrigin"))->Fill( rOri );
1246         ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("PTRes"))->Fill( mcPt, resPt );
1247         ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("ETARes"))->Fill( mcEta, resEta );
1248         ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("RAPRes"))->Fill( mcRap, resRap );
1249         ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("RXYRes"))->Fill( mcRxy, resRxy );
1250         ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("DLERes"))->Fill( mcDle, resDle );
1251         ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("APARes"))->Fill( mcApa, resApa );
1252         ((TH2D*)((TList*)fList->FindObject("RecMth"))->FindObject("APQRes"))->Fill( mcApq, resApq );
1253         LoadTrack(&ieT,iT->Chi2perNDF());
1254         ieT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1255         fDaughterImpactParameterXY = ip[0];
1256         fDaughterImpactParameterZ = ip[1];
1257         FillTrackSpy("TrkMth");
1258         LoadTrack(&jeT,jT->Chi2perNDF());
1259         jeT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1260         fDaughterImpactParameterXY = ip[0];
1261         fDaughterImpactParameterZ = ip[1];
1262         FillTrackSpy("TrkMth");
1263       }
1264       if(feeddown) {
1265         FillCandidateSpy("MthFDW");
1266         ((TH2D*)((TList*)fList->FindObject("MthFDW"))->FindObject("D0PD0N"))->Fill( dD0P,dD0N );
1267         ((TH1D*)((TList*)fList->FindObject("MthFDW"))->FindObject("MCOrigin"))->Fill( rOri );
1268       }
1269     }
1270     //===== END OF MCMATCH
1271   }
1272   ((TH2D*)((TList*)fList->FindObject("RecAll"))->FindObject("V0SADC"))->Fill( v0all,v0imw );
1273   return;
1274 }
1275 //=======================================================================
1276 Bool_t AliAnalysisTaskFlowStrange::PassesPIDCuts(AliESDtrack *myTrack, AliPID::EParticleType pid) {
1277   Bool_t pass=kTRUE;
1278   if(fPIDResponse) {
1279     fDaughterNSigmaPID = fPIDResponse->NumberOfSigmasTPC(myTrack,pid);
1280     if( TMath::Abs(fDaughterNSigmaPID) > fDaughterMaxNSigmaPID )
1281       pass = kFALSE;
1282   }
1283   return pass;
1284 }
1285 //=======================================================================
1286 void AliAnalysisTaskFlowStrange::ChargeParticles(AliAODEvent *tAOD) {
1287   //benchmark purposes (for the ultimate measurement for LHC10h see Alex, Francesco)
1288   if(!tAOD) return;
1289   // (for the moment) only compatible with SPVZE (no untagging)
1290   for(int i=0; i!=tAOD->GetNumberOfTracks(); ++i) {
1291     AliAODTrack *t = tAOD->GetTrack( i );
1292     if(!t) continue;
1293     if( !t->TestFilterBit(1024) ) continue;
1294     fDecayMass=0.0; // using mass as nsigmas control plot
1295     if(fPIDResponse) { // PID
1296       switch(fSpecie) { // TPC PID only
1297       case(kPION):
1298         fDecayMass = fPIDResponse->NumberOfSigmasTPC(t,AliPID::kPion);
1299         break;
1300       case(kKAON):
1301         fDecayMass = fPIDResponse->NumberOfSigmasTPC(t,AliPID::kKaon);
1302         break;
1303       case(kPROTON):
1304         fDecayMass = fPIDResponse->NumberOfSigmasTPC(t,AliPID::kProton);
1305         break;
1306       }
1307     }
1308     Bool_t pass = kTRUE;
1309     if( TMath::Abs(fDecayMass) > 3.0 ) pass=kFALSE;
1310     if( t->Eta()<-0.8 || t->Eta()>+0.8 ) pass=kFALSE;
1311     if( t->Pt()<0.2 || t->Pt()>20.0 ) pass=kFALSE;
1312     if( t->GetTPCsignal()<10.0 ) pass=kFALSE;
1313     fDecayPt=t->Pt();
1314     fDecayPhi=t->Phi();
1315     fDecayEta=t->Eta();
1316     fDecayID=t->GetID();
1317
1318     FillCandidateSpy("RecAll");
1319     if(!pass) continue;
1320
1321     AliESDtrack et( t );
1322     et.SetTPCClusterMap( t->GetTPCClusterMap() );
1323     et.SetTPCSharedMap( t->GetTPCSharedMap() );
1324     et.SetTPCPointsF( t->GetTPCNclsF() );
1325     Float_t ip[2];
1326     LoadTrack(&et,t->Chi2perNDF());
1327     AliAODVertex *vtx = tAOD->GetPrimaryVertex();
1328     Double_t pos[3];
1329     vtx->GetXYZ(pos);
1330     et.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1331     fDaughterImpactParameterXY = ip[0];
1332     fDaughterImpactParameterZ = ip[1];
1333     FillTrackSpy("TrkDau");
1334     FillCandidateSpy("RecSel");
1335
1336     MakeTrack();
1337   }
1338   return;
1339 }
1340 //=======================================================================
1341 void AliAnalysisTaskFlowStrange::Terminate(Option_t *) {
1342   //terminate
1343 }
1344 //=======================================================================
1345 void AliAnalysisTaskFlowStrange::MakeTrack() {
1346   // create track for flow tasks
1347   if(fCandidates->GetLast()+5>fCandidates->GetSize()) {
1348     fCandidates->Expand( fCandidates->GetSize()+20 );
1349   }
1350   Bool_t overwrite = kTRUE;
1351   AliFlowCandidateTrack *oTrack = (static_cast<AliFlowCandidateTrack*> (fCandidates->At( fCandidates->GetLast()+1 )));
1352   if( !oTrack ) { // creates new
1353     oTrack = new AliFlowCandidateTrack();
1354     overwrite = kFALSE;
1355   } else { // overwrites
1356     oTrack->ClearMe();
1357   }
1358   oTrack->SetMass(fDecayMass);
1359   oTrack->SetPt(fDecayPt);
1360   oTrack->SetPhi(fDecayPhi);
1361   oTrack->SetEta(fDecayEta);
1362   if(fSpecie<10) {
1363     oTrack->AddDaughter(fDecayIDpos);
1364     oTrack->AddDaughter(fDecayIDneg);
1365   } else {
1366     oTrack->SetID( fDecayID );
1367   }
1368   oTrack->SetForPOISelection(kTRUE);
1369   oTrack->SetForRPSelection(kFALSE);
1370   if(overwrite) {
1371     fCandidates->SetLast( fCandidates->GetLast()+1 );
1372   } else {
1373     fCandidates->AddLast(oTrack);
1374   }
1375   return;
1376 }
1377 //=======================================================================
1378 void AliAnalysisTaskFlowStrange::AddCandidates() {
1379   // adds candidates to flow events (untaging if necessary)
1380   if(fDebug) printf("FlowEventTPC %d tracks | %d RFP | %d POI\n",fTPCevent->NumberOfTracks(),fTPCevent->GetNumberOfRPs(),fTPCevent->GetNumberOfPOIs());
1381   if(fDebug) printf("FlowEventVZE %d tracks | %d RFP | %d POI\n",fVZEevent->NumberOfTracks(),fVZEevent->GetNumberOfRPs(),fVZEevent->GetNumberOfPOIs());
1382   if(fDebug) printf("I received %d candidates\n",fCandidates->GetEntriesFast());
1383   Int_t untagged=0;
1384   Int_t poi=0;
1385   for(int iCand=0; iCand!=fCandidates->GetEntriesFast(); ++iCand ) {
1386     AliFlowCandidateTrack *cand = static_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
1387     if(!cand) continue;
1388     cand->SetForPOISelection(kTRUE);
1389     cand->SetForRPSelection(kFALSE);
1390     poi++;
1391     if(fDebug) printf(" >Checking at candidate %d with %d daughters: mass %f\n",
1392                       iCand,cand->GetNDaughters(),cand->Mass());
1393     if(fSpecie<10) { // DECAYS
1394       // untagging ===>
1395       if(fDaughterUnTag) {
1396         for(int iDau=0; iDau!=cand->GetNDaughters(); ++iDau) {
1397           if(fDebug) printf("  >Daughter %d with fID %d", iDau, cand->GetIDDaughter(iDau));
1398           for(int iRPs=0; iRPs!=fTPCevent->NumberOfTracks(); ++iRPs ) {
1399             AliFlowTrack *iRP = static_cast<AliFlowTrack*>(fTPCevent->GetTrack( iRPs ));
1400             if(!iRP) continue;
1401             if(!iRP->InRPSelection()) continue;
1402             if(cand->GetIDDaughter(iDau) == iRP->GetID()) {
1403               if(fDebug) printf(" was in RP set");
1404               ++untagged;
1405               iRP->SetForRPSelection(kFALSE);
1406               fTPCevent->SetNumberOfRPs( fTPCevent->GetNumberOfRPs() -1 );
1407             }
1408           }
1409           if(fDebug) printf("\n");
1410         }
1411       }
1412       // <=== untagging
1413       fTPCevent->InsertTrack( ((AliFlowTrack*) cand) );
1414     } else {  // CHARGED
1415       // adding only new tracks and tagging accordingly ===>
1416       Bool_t found=kFALSE;
1417       for(int iRPs=0; iRPs!=fTPCevent->NumberOfTracks(); ++iRPs ) {
1418         AliFlowTrack *iRP = static_cast<AliFlowTrack*>(fTPCevent->GetTrack( iRPs ));
1419         if(!iRP) continue;
1420         if(!iRP->InRPSelection()) continue;
1421         if(cand->GetID() == iRP->GetID()) {
1422           if(fDebug) printf("  >charged track (%d) was also found in RP set (adding poi tag)\n",cand->GetID());
1423           iRP->SetForPOISelection(kTRUE);
1424           found = kTRUE;
1425         }
1426       }
1427       if(!found) // not found adding track
1428         fTPCevent->InsertTrack( ((AliFlowTrack*) cand) );
1429     }
1430     fVZEevent->InsertTrack( ((AliFlowTrack*) cand) );
1431   } //END OF LOOP
1432   fTPCevent->SetNumberOfPOIs( poi );
1433   fVZEevent->SetNumberOfPOIs( poi );
1434   ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("POI"))->Fill( poi );
1435   ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("UNTAG"))->Fill( untagged );
1436   if(fDebug) printf("FlowEventTPC %d tracks | %d RFP | %d POI\n",fTPCevent->NumberOfTracks(),fTPCevent->GetNumberOfRPs(),fTPCevent->GetNumberOfPOIs());
1437   if(fDebug) printf("FlowEventVZE %d tracks | %d RFP | %d POI\n",fVZEevent->NumberOfTracks(),fVZEevent->GetNumberOfRPs(),fVZEevent->GetNumberOfPOIs());
1438 }
1439 //=======================================================================
1440 void AliAnalysisTaskFlowStrange::PushBackFlowTrack(AliFlowEvent *flowevent, Double_t pt, Double_t phi, Double_t eta, Double_t w, Int_t id) {
1441   AliFlowTrack rfp;
1442   rfp.SetPt(pt);
1443   rfp.SetPhi(phi);
1444   rfp.SetEta(eta);
1445   rfp.SetWeight(w);
1446   rfp.SetForRPSelection(kTRUE);
1447   rfp.SetForPOISelection(kFALSE);
1448   rfp.SetID(id);
1449   flowevent->InsertTrack( &rfp );
1450 }
1451 //=======================================================================
1452 Bool_t AliAnalysisTaskFlowStrange::IsAtTPCEdge(Double_t phi,Double_t pt,Int_t charge,Double_t b) {
1453   // Origin: Alex Dobrin
1454   // Implemented by Carlos Perez
1455   TF1 cutLo("cutLo", "-0.01/x+pi/18.0-0.015", 0, 100);
1456   TF1 cutHi("cutHi", "0.55/x/x+pi/18.0+0.03", 0, 100);
1457   Double_t phimod = phi;
1458   if(b<0) phimod = TMath::TwoPi()-phimod;  //for negatve polarity field
1459   if(charge<0) phimod = TMath::TwoPi()-phimod; //for negatve charge
1460   phimod += TMath::Pi()/18.0;
1461   phimod = fmod(phimod, TMath::Pi()/9.0);
1462   if( phimod<cutHi.Eval(pt) && phimod>cutLo.Eval(pt) )
1463     return kTRUE;
1464
1465   return kFALSE;
1466 }
1467 //=======================================================================
1468 void AliAnalysisTaskFlowStrange::MakeQVectors() {
1469   //computes event plane and updates fPsi2
1470   //if there is a problem fPsi->-1
1471   fPsi2=-1;
1472   Double_t vzeqax, vzeqay, vzeqaw, vzeqbx, vzeqby, vzeqbw;
1473   Double_t tpcqax, tpcqay, tpcqaw, tpcqbx, tpcqby, tpcqbw;
1474   //=>loading event
1475   MakeQVZE(InputEvent(),vzeqax,vzeqay,vzeqaw,vzeqbx,vzeqby,vzeqbw);
1476   MakeQTPC(InputEvent(),tpcqax,tpcqay,tpcqaw,tpcqbx,tpcqby,tpcqbw);
1477   //=>computing psi
1478   //VZERO
1479   Double_t vqx, vqy;
1480   vqx=vqy=0;
1481   Double_t psivzea,psivzeb,psivze,vzew;
1482   psivzea = ( TMath::Pi()+TMath::ATan2(-vzeqay,-vzeqax) )/2.0;
1483   psivzeb = ( TMath::Pi()+TMath::ATan2(-vzeqby,-vzeqbx) )/2.0;
1484   vqx = vzeqax + vzeqbx;
1485   vqy = vzeqay + vzeqby;
1486   vzew = vzeqaw + vzeqbw;
1487   psivze = ( TMath::Pi()+TMath::ATan2(-vqy,-vqx) )/2.0;
1488   //TPC
1489   Double_t tqx, tqy;
1490   tqx=tqy=0;
1491   Double_t psitpca,psitpcb,psitpc,tpcw;
1492   psitpca = ( TMath::Pi()+TMath::ATan2(-tpcqay,-tpcqax) )/2.0;
1493   psitpcb = ( TMath::Pi()+TMath::ATan2(-tpcqby,-tpcqbx) )/2.0;
1494   tqx = tpcqax + tpcqbx;
1495   tqy = tpcqay + tpcqby;
1496   tpcw = tpcqaw + tpcqbw;
1497   psitpc = ( TMath::Pi()+TMath::ATan2(-tqy,-tqx) )/2.0;
1498   //=>does the event clear?
1499   switch(fWhichPsi) {
1500   case(1): //VZERO
1501     if( (vzeqaw<2)||(vzeqbw<2) ) return;
1502     fPsi2 = psivze;
1503     break;
1504   case(2): //TPC
1505     if( (tpcqaw<2)||(tpcqbw<2) ) return;
1506     fPsi2 = psitpc;
1507     break;
1508   }
1509   //=>great! recording
1510   ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEPSI"))->Fill( psivze );
1511   ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEPSIA"))->Fill( psivzea );
1512   ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEPSIB"))->Fill( psivzeb );
1513   ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("RFPVZE"))->Fill( vzew );
1514   //------
1515   ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCPSI"))->Fill( psitpc );
1516   ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCPSIA"))->Fill( psitpca );
1517   ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCPSIB"))->Fill( psitpcb );
1518   ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("RFPTPC"))->Fill( tpcw );
1519   //------
1520   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQ"))->Fill( 1., tpcqay/tpcqaw, tpcqaw );
1521   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQ"))->Fill( 2., tpcqax/tpcqaw, tpcqaw );
1522   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQ"))->Fill( 3., tpcqby/tpcqbw, tpcqbw );
1523   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQ"))->Fill( 4., tpcqbx/tpcqbw, tpcqbw );
1524   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQ"))->Fill( 5., tqy/tpcw, tpcw );
1525   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQ"))->Fill( 6., tqx/tpcw, tpcw );
1526   //------
1527   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQ"))->Fill( 1., vzeqay/vzeqaw, vzeqaw );
1528   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQ"))->Fill( 2., vzeqax/vzeqaw, vzeqaw );
1529   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQ"))->Fill( 3., vzeqby/vzeqbw, vzeqbw );
1530   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQ"))->Fill( 4., vzeqbx/vzeqbw, vzeqbw );
1531   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQ"))->Fill( 5., vqy/vzew, vzew );
1532   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQ"))->Fill( 6., vqx/vzew, vzew );
1533
1534   return;
1535 }
1536 //=======================================================================
1537 void AliAnalysisTaskFlowStrange::MakeQVZE(AliVEvent *tevent,
1538                                           Double_t &qxa,Double_t &qya,Double_t &qwa,
1539                                           Double_t &qxb,Double_t &qyb,Double_t &qwb) {
1540   //=>cleaning
1541   if(fUseFP) fVZEevent->ClearFast();
1542   //=>external weights
1543   Double_t extW[64];
1544   for(int i=0;i!=64;++i) extW[i]=1;
1545   if((!fVZEsave)&&(fVZEResponse)) {
1546     Double_t minC = fCentPerMin, maxC = fCentPerMax;
1547     if(fVZEmb) {
1548       minC = 0;
1549       maxC = 80;
1550     }
1551     Int_t ybinmin = fVZEResponse->GetYaxis()->FindBin(minC+1e-6);
1552     Int_t ybinmax = fVZEResponse->GetYaxis()->FindBin(maxC-1e-6);
1553     for(int i=0;i!=64;++i) extW[i] = fVZEResponse->Integral(i+1,i+1,ybinmin,ybinmax)/(maxC-minC);
1554     //ring-wise normalization
1555     Double_t ring[8];
1556     for(int j=0; j!=8; ++j) {
1557       ring[j]=0;
1558       for(int i=0;i!=8;++i) ring[j] += extW[j*8+i]/8;
1559     }
1560     //disk-wise normalization
1561     Double_t disk[2];
1562     int xbinmin, xbinmax;
1563     xbinmin = 1+8*fVZECa;
1564     xbinmax = 8+8*fVZECb;
1565     disk[0] = fVZEResponse->Integral(xbinmin,xbinmax,ybinmin,ybinmax)/(maxC-minC)/(xbinmax-xbinmin+1);
1566     xbinmin = 33+8*fVZEAa;
1567     xbinmax = 40+8*fVZEAb;
1568     disk[1] = fVZEResponse->Integral(xbinmin,xbinmax,ybinmin,ybinmax)/(maxC-minC)/(xbinmax-xbinmin+1);
1569     //for(int i=0;i!=64;++i) printf("CELL %d -> W = %f ||",i,extW[i]);
1570
1571     if(fVZEByDisk) {
1572       for(int i=0;i!=64;++i) extW[i] = disk[i/32]/extW[i];
1573     } else {
1574       for(int i=0;i!=64;++i) extW[i] = ring[i/8]/extW[i];
1575     }
1576     //for(int i=0;i!=64;++i) printf(" W = %f \n",extW[i]);
1577   }
1578   //=>computing
1579   qxa=qya=qwa=qxb=qyb=qwb=0;
1580   Int_t rfp=0;
1581   Double_t eta, phi, w;
1582   //v0c -> qa
1583   for(int id=fVZECa*8;id!=8+fVZECb*8;++id) {
1584     eta = -3.45+0.5*(id/8);
1585     phi = TMath::PiOver4()*(0.5+id%8);
1586     w = tevent->GetVZEROEqMultiplicity(id)*extW[id];
1587     qxa += w*TMath::Cos(2*phi);
1588     qya += w*TMath::Sin(2*phi);
1589     qwa += w;
1590     ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEAllPhiEta"))->Fill( phi, eta, w );
1591     rfp++;
1592     if(fUseFP) PushBackFlowTrack(fVZEevent,0,phi,eta,w,0);
1593   }
1594   //v0a -> qb
1595   for(int id=32+fVZEAa*8;id!=40+fVZEAb*8;++id) {
1596     eta = +4.8-0.6*((id/8)-4);
1597     phi = TMath::PiOver4()*(0.5+id%8);
1598     w = tevent->GetVZEROEqMultiplicity(id)*extW[id];
1599     qxb += w*TMath::Cos(2*phi);
1600     qyb += w*TMath::Sin(2*phi);
1601     qwb += w;
1602     ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEAllPhiEta"))->Fill( phi, eta, w );
1603     rfp++;
1604     if(fUseFP) PushBackFlowTrack(fVZEevent,0,phi,eta,w,0);
1605   }
1606   if(fUseFP) fVZEevent->SetNumberOfRPs(rfp);
1607   if(fDebug>0&&fUseFP) printf("Inserted tracks in FlowEventVZE %d ==> %.1f\n",fVZEevent->NumberOfTracks(),qwa+qwb);
1608 }
1609 //=======================================================================
1610 void AliAnalysisTaskFlowStrange::AddTPCRFPSpy(TList *me) {
1611   TH1D *tH1D;
1612   tH1D = new TH1D("PT",   "PT",           50,0,5);     me->Add(tH1D);
1613   tH1D = new TH1D("PHI",  "PHI", 90,0,TMath::TwoPi()); me->Add(tH1D);
1614   tH1D = new TH1D("ETA",  "ETA",          40,-1,+1);   me->Add(tH1D);
1615   tH1D = new TH1D("TPCS", "TPC Signal",   100,0,500);  me->Add(tH1D);
1616   tH1D = new TH1D("IPXY", "IPXY",         100,-2,+2);  me->Add(tH1D);
1617   tH1D = new TH1D("IPZ",  "IPZ",          100,-2,+2);  me->Add(tH1D);
1618   // TPC
1619   tH1D = new TH1D("TPCNCLS", "NCLS", 170,-0.5,+169.5);   me->Add(tH1D);
1620   tH1D = new TH1D("TPCSHCL", "NSCLS / NCLS", 100,0,1);   me->Add(tH1D);
1621   tH1D = new TH1D("TPCFICL", "NCLS1I / NCLS",100,0,1);   me->Add(tH1D);
1622   tH1D = new TH1D("TPCXRNF", "XROW / NFCLS", 100,0,1.5); me->Add(tH1D);
1623   tH1D = new TH1D("TPCRCHI", "CHI2 / NCLS",  50,0,5);    me->Add(tH1D);
1624   // ITS
1625   tH1D = new TH1D("ITSNCLS", "NCLS",   7,-0.5,+6.5); me->Add(tH1D);
1626   tH1D = new TH1D("ITSRCHI", "CHI2 / NCLS", 50,0,5); me->Add(tH1D);
1627
1628 }
1629 //=======================================================================
1630 Bool_t AliAnalysisTaskFlowStrange::PassesRFPTPCCuts(AliESDtrack *track, Double_t aodchi2cls, Float_t aodipxy, Float_t aodipz) {
1631   if(track->GetKinkIndex(0)>0) return kFALSE;
1632   if( (track->GetStatus()&AliESDtrack::kTPCrefit)==0 ) return kFALSE;
1633   Double_t pt = track->Pt();
1634   Double_t phi = track->Phi();
1635   Double_t eta = track->Eta();
1636   Double_t tpcs = track->GetTPCsignal();
1637   Float_t ipxy, ipz;
1638   track->GetImpactParameters(ipxy,ipz);
1639   Int_t cls = track->GetTPCclusters(0);
1640   Double_t xrows, findcls, chi2;
1641   findcls = track->GetTPCNclsF();
1642   xrows = track->GetTPCCrossedRows();
1643   chi2 = track->GetTPCchi2();
1644   Double_t rchi2 = chi2/cls;
1645   if(!fReadESD) {
1646     rchi2 = aodchi2cls;
1647     ipxy = aodipxy;
1648     ipz = aodipz;
1649   }
1650   Double_t xrnfcls = xrows/findcls;
1651   Double_t scls, cls1i, itschi2;
1652   Int_t itscls;
1653   cls1i = track->GetTPCNclsIter1();
1654   scls = track->GetTPCnclsS();
1655   itscls = track->GetITSclusters(0);
1656   itschi2 = track->GetITSchi2();
1657   Double_t shcl = scls/cls;
1658   Double_t ficl = cls1i/cls;
1659   Double_t itsrchi2 = itscls/itschi2;
1660   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("PT"))->Fill( pt );
1661   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("PHI"))->Fill( phi );
1662   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("ETA"))->Fill( eta );
1663   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCS"))->Fill( tpcs );
1664   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("IPXY"))->Fill( ipxy );
1665   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("IPZ"))->Fill( ipz );
1666   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCNCLS"))->Fill( cls );
1667   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCSHCL"))->Fill( shcl );
1668   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCFICL"))->Fill( ficl );
1669   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCXRNF"))->Fill( xrnfcls );
1670   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCRCHI"))->Fill( rchi2 );
1671   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("ITSNCLS"))->Fill( itscls );
1672   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("ITSRCHI"))->Fill( itsrchi2 );
1673   if(pt<fRFPminPt) return kFALSE; //0.2
1674   if(pt>fRFPmaxPt) return kFALSE; //5.0
1675   if(eta<fRFPminEta) return kFALSE; //-0.8
1676   if(eta>fRFPmaxEta) return kFALSE; //+0.8
1677   if(tpcs<fRFPTPCsignal) return kFALSE; //10.0
1678   if( TMath::Sqrt(ipxy*ipxy/fRFPmaxIPxy/fRFPmaxIPxy+ipz*ipz/fRFPmaxIPz/fRFPmaxIPz)>1 ) return kFALSE; //2.4 3.2
1679   if(cls<fRFPTPCncls) return kFALSE; //70
1680   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("PT"))->Fill( pt );
1681   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("PHI"))->Fill( phi );
1682   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("ETA"))->Fill( eta );
1683   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCS"))->Fill( tpcs );
1684   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("IPXY"))->Fill( ipxy );
1685   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("IPZ"))->Fill( ipz );
1686   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCNCLS"))->Fill( cls );
1687   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCSHCL"))->Fill( shcl );
1688   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCFICL"))->Fill( ficl );
1689   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCXRNF"))->Fill( xrnfcls );
1690   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCRCHI"))->Fill( rchi2 );
1691   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("ITSNCLS"))->Fill( itscls );
1692   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("ITSRCHI"))->Fill( itsrchi2 );
1693   return kTRUE;
1694 }
1695 //=======================================================================
1696 void AliAnalysisTaskFlowStrange::MakeQTPC(AliVEvent *tevent,
1697                                           Double_t &qxa,Double_t &qya,Double_t &qwa,
1698                                           Double_t &qxb,Double_t &qyb,Double_t &qwb) {
1699   AliESDEvent *tESD = (AliESDEvent*) (tevent);
1700   AliAODEvent *tAOD = (AliAODEvent*) (tevent);
1701   if(fReadESD) {
1702     if(!tESD) return;
1703     MakeQTPC(tESD,qxa,qya,qwa,qxb,qyb,qwb);
1704   } else {
1705     if(!tAOD) return;
1706     MakeQTPC(tAOD,qxa,qya,qwa,qxb,qyb,qwb);
1707   }
1708 }
1709 //=======================================================================
1710 void AliAnalysisTaskFlowStrange::MakeQTPC(AliAODEvent *tAOD,
1711                                           Double_t &qxa,Double_t &qya,Double_t &qwa,
1712                                           Double_t &qxb,Double_t &qyb,Double_t &qwb) {
1713   //=>cleaning
1714   if(fUseFP) fTPCevent->ClearFast();
1715   qxa=qya=qwa=qxb=qyb=qwb=0;
1716   Int_t rfp=0;
1717   Double_t eta, phi, w;
1718   //=>aod stuff
1719   AliAODVertex *vtx = tAOD->GetPrimaryVertex();
1720   Double_t pos[3],cov[6];
1721   vtx->GetXYZ(pos);
1722   vtx->GetCovarianceMatrix(cov);
1723   const AliESDVertex vESD(pos,cov,100.,100);
1724   AliAODTrack *track;
1725   //=>looping
1726   Int_t rawN = tAOD->GetNumberOfTracks();
1727   for(Int_t id=0; id!=rawN; ++id) {
1728     track = tAOD->GetTrack(id);
1729     //=>cuts
1730     if(!track->TestFilterBit(fRFPFilterBit)) continue;
1731     if( fExcludeTPCEdges )
1732       if( IsAtTPCEdge( track->Phi(), track->Pt(), track->Charge(), tAOD->GetMagneticField() ) ) continue;
1733     AliESDtrack etrack( track );
1734     etrack.SetTPCClusterMap( track->GetTPCClusterMap() );
1735     etrack.SetTPCSharedMap( track->GetTPCSharedMap() );
1736     etrack.SetTPCPointsF( track->GetTPCNclsF() );
1737     Float_t ip[2];
1738     etrack.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1739     if(!PassesRFPTPCCuts(&etrack,track->Chi2perNDF(),ip[0],ip[1])) continue;
1740     //=>collecting info
1741     phi = track->Phi();
1742     eta = track->Eta();
1743     w = 1;
1744     //=>subevents
1745     if(eta<0) {
1746       qxa += w*TMath::Cos(2*phi);
1747       qya += w*TMath::Sin(2*phi);
1748       qwa += w;
1749     } else {
1750       qxb += w*TMath::Cos(2*phi);
1751       qyb += w*TMath::Sin(2*phi);
1752       qwb += w;
1753     }
1754     ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCAllPhiEta"))->Fill( phi, eta, w );
1755     rfp++;
1756     if(fUseFP) PushBackFlowTrack(fTPCevent,track->Pt(),phi,eta,w,track->GetID());
1757   }
1758   if(fUseFP) fTPCevent->SetNumberOfRPs(rfp);
1759   if(fDebug) printf("Inserted tracks in FlowEventTPC %d ==> %.1f\n",fTPCevent->NumberOfTracks(),qwa+qwb);
1760 }
1761 //=======================================================================
1762 void AliAnalysisTaskFlowStrange::MakeQTPC(AliESDEvent *tESD,
1763                                           Double_t &qxa,Double_t &qya,Double_t &qwa,
1764                                           Double_t &qxb,Double_t &qyb,Double_t &qwb) {
1765   //=>cleaning
1766   if(fUseFP) fTPCevent->ClearFast();
1767   qxa=qya=qwa=qxb=qyb=qwb=0;
1768   Int_t rfp=0;
1769   Double_t eta, phi, w;
1770   //=>looping
1771   AliESDtrack *track;
1772   Int_t rawN = tESD->GetNumberOfTracks();
1773   for(Int_t id=0; id!=rawN; ++id) {
1774     track = tESD->GetTrack(id);
1775     //=>cuts
1776     if( fExcludeTPCEdges )
1777       if( IsAtTPCEdge( track->Phi(), track->Pt(), track->Charge(), tESD->GetMagneticField() ) ) continue;
1778     if(!PassesFilterBit(track)) continue;
1779     if(!PassesRFPTPCCuts(track)) continue;
1780     //=>collecting info
1781     phi = track->Phi();
1782     eta = track->Eta();
1783     w = 1;
1784     //=>subevents
1785     if(eta<0) {
1786       qxa += w*TMath::Cos(2*phi);
1787       qya += w*TMath::Sin(2*phi);
1788       qwa += w;
1789     } else {
1790       qxb += w*TMath::Cos(2*phi);
1791       qyb += w*TMath::Sin(2*phi);
1792       qwb += w;
1793     }
1794     ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCAllPhiEta"))->Fill( phi, eta, w );
1795     rfp++;
1796     if(fUseFP) PushBackFlowTrack(fTPCevent,track->Pt(),phi,eta,w,track->GetID());
1797   }
1798   if(fUseFP) fTPCevent->SetNumberOfRPs(rfp);
1799   if(fDebug) printf("Inserted tracks in FlowEventTPC %d ==> %.1f\n",fTPCevent->NumberOfTracks(),qwa+qwb);
1800 }
1801 //=======================================================================
1802 void AliAnalysisTaskFlowStrange::AddMCParticleSpy(TList *me) {
1803   TH1D *tH1D;
1804   TH2D *tH2D;
1805   tH1D = new TH1D("Pt",   "Pt",   100,0.0,20);  me->Add(tH1D);
1806   tH1D = new TH1D("Phi",  "Phi",  100,0,TMath::TwoPi()); me->Add(tH1D);
1807   tH1D = new TH1D("Eta",  "Eta",  100,-1,+1);   me->Add(tH1D);
1808   tH1D = new TH1D("Rad2", "Rad2", 1000,0,+100); me->Add(tH1D);
1809   tH2D = new TH2D("Dphi", "phi-MCEP;pt;dphi",100,0,20, 72,0,TMath::Pi()); me->Add(tH2D);
1810   return;
1811 }
1812 //=======================================================================
1813 void AliAnalysisTaskFlowStrange::FillMCParticleSpy(TString listName, AliAODMCParticle *p) {
1814   ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Pt" ))->Fill( p->Pt() );
1815   ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Eta" ))->Fill( p->Eta() );
1816   ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Phi" ))->Fill( p->Phi() );
1817   ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Rad2" ))->Fill( TMath::Sqrt( p->Xv()*p->Xv() +
1818                                                                                                  p->Yv()*p->Yv() ) );
1819   ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Dphi" ))->Fill( p->Pt(), GetMCDPHI(p->Phi()) );
1820   return;
1821 }
1822 //=======================================================================
1823 Double_t AliAnalysisTaskFlowStrange::GetMCDPHI(Double_t phi) {
1824   Double_t dDPHI = phi - fMCEP;
1825   if( dDPHI < 0 ) dDPHI += TMath::TwoPi();
1826   if( dDPHI > TMath::Pi() ) dDPHI = TMath::TwoPi()-dDPHI;
1827   return dDPHI;
1828 }
1829 //=======================================================================
1830 void AliAnalysisTaskFlowStrange::FillMCParticleSpy(TString listName, TParticle *p) {
1831   ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Pt" ))->Fill( p->Pt() );
1832   ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Eta" ))->Fill( p->Eta() );
1833   ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Phi" ))->Fill( p->Phi() );
1834   ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Rad2" ))->Fill( TMath::Sqrt( p->Vx()*p->Vx() +
1835                                                                                                  p->Vy()*p->Vy() ) );
1836   return;
1837 }
1838 //=======================================================================
1839 void AliAnalysisTaskFlowStrange::AddCandidatesSpy(TList *me) {
1840   TH2D *tH2D;
1841   tH2D = new TH2D("PhiEta",  "PhiEta;Phi;Eta", 100,0,TMath::TwoPi(),100,-1,+1); me->Add(tH2D);
1842   tH2D = new TH2D("PtRAP",   "PtRAP;Pt;Y",          100,0,20,100,-2.0,+2.0);  me->Add(tH2D);
1843   tH2D = new TH2D("PtDCA",   "PtDCA;Pt;DCA",        100,0,20,100,0,10);       me->Add(tH2D);
1844   tH2D = new TH2D("PtCTP",   "PtCTP;Pt;CTP",        100,0,20,100,-1,+1);      me->Add(tH2D);
1845   //tH2D = new TH2D("PtCTP",   "PtCTP;Pt;CTP",        100,0,10,100,0.90,+1);    me->Add(tH2D);
1846   tH2D = new TH2D("PtD0D0",  "PtD0D0;Pt;D0D0",      100,0,20,100,-5,+5);      me->Add(tH2D);
1847   tH2D = new TH2D("PtRad2",  "PtRad2;Pt;RadXY",     100,0,20,100,0,+50);      me->Add(tH2D);
1848   tH2D = new TH2D("PtDL",    "PtDL;Pt;DL",          100,0,20,100,0,+50);      me->Add(tH2D);
1849   tH2D = new TH2D("PtMASS",  "PtMASS;Pt;MASS", 100,0,20,fMassBins,fMinMass,fMaxMass); me->Add(tH2D);
1850   tH2D = new TH2D("APPOS",   "APPOS;alphaPOS;QtPOS",100,-2,+2,100,0,0.3);     me->Add(tH2D);
1851   tH2D = new TH2D("D0PD0N",  "D0PD0N;D0P;D0N",      200,-10,+10,200,-10,+10); me->Add(tH2D);
1852   tH2D = new TH2D("XPOSXNEG","XPOSXNEG;XPOS;XNEG",  200,-50,+50,200,-50,+50); me->Add(tH2D);
1853   tH2D = new TH2D("PTDPHIMC","PtDPHIMC;Pt;PHI-MCEP",100,0,20,72,0,TMath::Pi()); me->Add(tH2D);
1854   return;
1855 }
1856 //=======================================================================
1857 void AliAnalysisTaskFlowStrange::FillCandidateSpy(TString listName) {
1858   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PhiEta"))->Fill( fDecayPhi, fDecayEta );
1859   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtRAP" ))->Fill( fDecayPt, fDecayRapidity );
1860   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtDCA" ))->Fill( fDecayPt, fDecayDCAdaughters );
1861   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtCTP" ))->Fill( fDecayPt, fDecayCosinePointingAngleXY );
1862   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtD0D0"))->Fill( fDecayPt, fDecayProductIPXY );
1863   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtRad2"))->Fill( fDecayPt, fDecayRadXY );
1864   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtDL"  ))->Fill( fDecayPt, fDecayDecayLength );
1865   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtMASS"))->Fill( fDecayPt, fDecayMass );
1866   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("APPOS" ))->Fill( fDecayAlpha, fDecayQt );
1867   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PTDPHIMC" ))->Fill( fDecayPt, GetMCDPHI( fDecayPhi ) );
1868 }
1869 //=======================================================================
1870 Bool_t AliAnalysisTaskFlowStrange::AcceptCandidate() {
1871   if(fDecayEta<fDecayMinEta) return kFALSE;
1872   if(fDecayEta>fDecayMaxEta) return kFALSE;
1873   if(fDecayPt<fDecayMinPt) return kFALSE;
1874   if(fDecayProductIPXY>fDecayMaxProductIPXY) return kFALSE;
1875   if(fDecayDCAdaughters>fDecayMaxDCAdaughters) return kFALSE;
1876   if(fDecayCosinePointingAngleXY<fDecayMinCosinePointingAngleXY) return kFALSE;
1877   if(fDecayRadXY<fDecayMinRadXY) return kFALSE;
1878   if(fDecayDecayLength>fDecayMaxDecayLength*2.6842) return kFALSE;
1879   if(TMath::Abs(fDecayRapidity)>fDecayMaxRapidity) return kFALSE;
1880   if(fSpecie==0) {
1881     if(fDecayAPCutPie) {
1882       if(fDecayQt/TMath::Abs(fDecayAlpha)<fDecayMinQt) return kFALSE;
1883     } else {
1884       if(fDecayQt<fDecayMinQt) return kFALSE;
1885     }
1886   }
1887   if(fSpecie==1) if(fDecayAlpha>0) return kFALSE;
1888   if(fSpecie==2) if(fDecayAlpha<0) return kFALSE;
1889   return kTRUE;
1890 }
1891 //=======================================================================
1892 void AliAnalysisTaskFlowStrange::AddTracksSpy(TList *me) {
1893   TH2D *tH2D;
1894   tH2D = new TH2D("PHIETA",       "PHIETA;PHI;ETA",               100,0,TMath::TwoPi(),100,-2,2); me->Add(tH2D);
1895   tH2D = new TH2D("IPXYIPZ",      "IPXYIPZ;IPXY;IPZ",             1000,-20,+20,1000,-20,+20); me->Add(tH2D);
1896   tH2D = new TH2D("PTTPCNCLS",    "PTTPCNCLS;PT;NCLS",            100,0,20,170,0,170);  me->Add(tH2D);
1897   tH2D = new TH2D("POSTPCNCLCHI2","POSTPCNCLCHI2;NCLS;CHI2/NCLS", 170,0,170,100,0,8);   me->Add(tH2D);
1898   tH2D = new TH2D("POSTPCNFCLNXR","POSTPCNFCLNXR;NFCLS;NXR",      170,0,170,170,0,170); me->Add(tH2D);
1899   tH2D = new TH2D("POSTPCNCLNFCL","POSTPCNCLNFCL;NCLS;NFCLS",     170,0,170,170,0,170); me->Add(tH2D);
1900   tH2D = new TH2D("POSTPCNCLNSCL","POSTPCNCLNSCL;NCLS;NSCLS",     170,0,170,170,0,170); me->Add(tH2D);
1901   tH2D = new TH2D("NEGTPCNCLCHI2","NEGTPCNCLCHI2;NCLS;CHI2/NCLS", 170,0,170,100,0,8);   me->Add(tH2D);
1902   tH2D = new TH2D("NEGTPCNFCLNXR","NEGTPCNFCLNXR;NFCLS;NXR",      170,0,170,170,0,170); me->Add(tH2D);
1903   tH2D = new TH2D("NEGTPCNCLNFCL","NEGTPCNCLNFCL;NCLS;NFCLS",     170,0,170,170,0,170); me->Add(tH2D);
1904   tH2D = new TH2D("NEGTPCNCLNSCL","NEGTPCNCLNSCL;NCLS;NSCLS",     170,0,170,170,0,170); me->Add(tH2D);
1905 }
1906 //=======================================================================
1907 void AliAnalysisTaskFlowStrange::FillTrackSpy(TString listName) {
1908   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PHIETA" ))->Fill( fDaughterPhi, fDaughterEta );
1909   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "IPXYIPZ" ))->Fill( fDaughterImpactParameterXY, fDaughterImpactParameterZ );
1910   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTTPCNCLS" ))->Fill( fDaughterPt, fDaughterNClsTPC );
1911   TString ch="NEG";
1912   if(fDaughterCharge>0) ch="POS";
1913   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( Form("%sTPCNCLCHI2",ch.Data()) ))->Fill( fDaughterNClsTPC, fDaughterChi2PerNClsTPC );
1914   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( Form("%sTPCNFCLNXR",ch.Data()) ))->Fill( fDaughterNFClsTPC, fDaughterXRows );
1915   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( Form("%sTPCNCLNFCL",ch.Data()) ))->Fill( fDaughterNClsTPC, fDaughterNFClsTPC );
1916   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( Form("%sTPCNCLNSCL",ch.Data()) ))->Fill( fDaughterNClsTPC, fDaughterNSClsTPC );
1917 }
1918 //=======================================================================
1919 void AliAnalysisTaskFlowStrange::LoadTrack(AliESDtrack *myTrack, Double_t aodChi2NDF) {
1920   fDaughterCharge = myTrack->Charge();
1921   fDaughterXRows = myTrack->GetTPCCrossedRows();
1922   fDaughterNFClsTPC = myTrack->GetTPCNclsF();
1923   fDaughterNSClsTPC = myTrack->GetTPCnclsS();
1924   fDaughterNClsTPC = myTrack->GetTPCclusters(0);
1925   if(fReadESD) {
1926     if(fDaughterNClsTPC>0) fDaughterChi2PerNClsTPC = myTrack->GetTPCchi2()/fDaughterNClsTPC;
1927   } else {
1928     fDaughterChi2PerNClsTPC = aodChi2NDF;
1929   }
1930   myTrack->GetImpactParameters(fDaughterImpactParameterXY,fDaughterImpactParameterZ);
1931   fDaughterStatus = myTrack->GetStatus();
1932   fDaughterPhi = myTrack->Phi();
1933   fDaughterEta = myTrack->Eta();
1934   fDaughterPt = myTrack->Pt();
1935   fDaughterKinkIndex = myTrack->GetKinkIndex(0);
1936 }
1937 //=======================================================================
1938 Bool_t AliAnalysisTaskFlowStrange::AcceptDaughter() {
1939   if(fDaughterKinkIndex>0) return kFALSE;
1940   if( (fDaughterStatus&AliESDtrack::kTPCrefit)==0 ) return kFALSE;
1941   if(fDaughterNFClsTPC<1) return kFALSE;
1942   if(fDaughterPt<fDaughterMinPt) return kFALSE;
1943   if(fDaughterEta<fDaughterMinEta) return kFALSE;
1944   if(fDaughterEta>fDaughterMaxEta) return kFALSE;
1945   if(fDaughterNClsTPC<fDaughterMinNClsTPC) return kFALSE;
1946   if(fDaughterXRows<fDaughterMinXRows) return kFALSE;
1947   if(fDaughterChi2PerNClsTPC>fDaughterMaxChi2PerNClsTPC) return kFALSE;
1948   if(TMath::Abs(fDaughterImpactParameterXY)<fDaughterMinImpactParameterXY) return kFALSE;
1949   if(fDaughterXRows<fDaughterMinXRowsOverNClsFTPC*fDaughterNFClsTPC) return kFALSE;
1950   return kTRUE;
1951 }
1952 //=======================================================================
1953 Double_t AliAnalysisTaskFlowStrange::GetWDist(const AliVVertex* v0, const AliVVertex* v1) {
1954   // calculate sqrt of weighted distance to other vertex
1955   if (!v0 || !v1) {
1956     printf("One of vertices is not valid\n");
1957     return 0;
1958   }
1959   static TMatrixDSym vVb(3);
1960   double dist = -1;
1961   double dx = v0->GetX()-v1->GetX();
1962   double dy = v0->GetY()-v1->GetY();
1963   double dz = v0->GetZ()-v1->GetZ();
1964   double cov0[6],cov1[6];
1965   v0->GetCovarianceMatrix(cov0);
1966   v1->GetCovarianceMatrix(cov1);
1967   vVb(0,0) = cov0[0]+cov1[0];
1968   vVb(1,1) = cov0[2]+cov1[2];
1969   vVb(2,2) = cov0[5]+cov1[5];
1970   vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
1971   vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
1972   vVb.InvertFast();
1973   if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
1974   dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
1975     +    2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
1976   return dist>0 ? TMath::Sqrt(dist) : -1;
1977 }
1978 //=======================================================================
1979 Bool_t AliAnalysisTaskFlowStrange::plpMV(const AliVEvent *event) {
1980   // check for multi-vertexer pile-up
1981   const AliAODEvent *aod = (const AliAODEvent*)event;
1982   const AliESDEvent *esd = (const AliESDEvent*)event;
1983   //
1984   const int    kMinPlpContrib = 5;
1985   const double kMaxPlpChi2 = 5.0;
1986   const double kMinWDist = 15;
1987   //
1988   if (!aod && !esd) {
1989     printf("Event is neither of AOD nor ESD\n");
1990     exit(1);
1991   }
1992   //
1993   const AliVVertex* vtPrm = 0;
1994   const AliVVertex* vtPlp = 0;
1995   int nPlp = 0;
1996   //
1997   if (aod) {
1998     if ( !(nPlp=aod->GetNumberOfPileupVerticesTracks()) ) return kFALSE;
1999     vtPrm = aod->GetPrimaryVertex();
2000     if (vtPrm == aod->GetPrimaryVertexSPD()) return kTRUE; // there are pile-up vertices but no primary
2001   }
2002   else {
2003     if ( !(nPlp=esd->GetNumberOfPileupVerticesTracks())) return kFALSE;
2004     vtPrm = esd->GetPrimaryVertexTracks();
2005     if (((AliESDVertex*)vtPrm)->GetStatus()!=1) return kTRUE; // there are pile-up vertices but no primary
2006   }
2007   //int bcPrim = vtPrm->GetBC();
2008   //
2009   for (int ipl=0;ipl<nPlp;ipl++) {
2010     vtPlp = aod ? (const AliVVertex*)aod->GetPileupVertexTracks(ipl) : (const AliVVertex*)esd->GetPileupVertexTracks(ipl);
2011     //
2012     if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
2013     if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
2014     //  int bcPlp = vtPlp->GetBC();
2015     //  if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other BC
2016     //
2017     double wDst = GetWDist(vtPrm,vtPlp);
2018     if (wDst<kMinWDist) continue;
2019     //
2020     return kTRUE; // pile-up: well separated vertices
2021   }
2022   //
2023   return kFALSE;
2024   //
2025 }
2026 //=======================================================================
2027 void AliAnalysisTaskFlowStrange::MakeFilterBits() {
2028   //FilterBit 1
2029   fFB1 = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
2030   //FilterBit1024
2031   fFB1024 = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE,1);
2032   fFB1024->SetMinNCrossedRowsTPC(120);
2033   fFB1024->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
2034   fFB1024->SetMaxChi2PerClusterITS(36);
2035   fFB1024->SetMaxFractionSharedTPCClusters(0.4);
2036   fFB1024->SetMaxChi2TPCConstrainedGlobal(36);
2037   fFB1024->SetEtaRange(-0.9,0.9);
2038   fFB1024->SetPtRange(0.15, 1e10);
2039 }
2040 //=======================================================================
2041 Bool_t AliAnalysisTaskFlowStrange::PassesFilterBit(AliESDtrack *track) {
2042   Bool_t ret=kFALSE;
2043   switch(fRFPFilterBit) {
2044     case(1024):
2045       ret = fFB1024->AcceptTrack(track);
2046       break;
2047     default:
2048       ret = fFB1->AcceptTrack(track);
2049   }
2050   return ret;
2051 }
2052 //=======================================================================
2053 void AliAnalysisTaskFlowStrange::LoadVZEROResponse() {
2054   if(fVZEResponse) {
2055     TString run = fVZEResponse->GetTitle();
2056     if( run.Atoi() == fRunNumber ) return;
2057     fVZEResponse = NULL;
2058   }
2059   //==>loading
2060   fVZEResponse = dynamic_cast<TH2D*> (fVZEload->FindObject( Form("%d",fRunNumber) ));
2061   if(fVZEResponse) {
2062     printf("New VZE calibration: run %d || %s -> Entries %.0f\n",fRunNumber, fVZEResponse->GetTitle(),fVZEResponse->GetEntries());
2063   } else {
2064     printf("VZE calibration: request but not found!!!\n");
2065   }
2066 }
2067 //=======================================================================
2068 void AliAnalysisTaskFlowStrange::AddVZEQA() {
2069   fVZEQA = new TList();
2070   fVZEQA->SetName( Form("VZEQA%d",fRunNumber) );
2071   fVZEQA->SetOwner();
2072   if(fQAlevel>0) {
2073     TProfile2D *prof = new TProfile2D("LINP","LINP;VZEcell;VZEmult;SPDtrkl", 64,0,64,500,0,700,0,10000); fVZEQA->Add( prof );
2074     prof = new TProfile2D("MULP","MULP;VZEcell;CENTR;VZEmult", 64,0,64,100,0,100,0,10000); fVZEQA->Add( prof );
2075     TH3D *tH3D = new TH3D("EQU","EQU;VZEeqmult;VZEmult",100,0,700,100,0,700,64,0,64); fVZEQA->Add( tH3D );
2076   }
2077   fList->Add(fVZEQA);
2078 }
2079 //=======================================================================
2080 void AliAnalysisTaskFlowStrange::SaveVZEROQA() {
2081   AliAODEvent *event = dynamic_cast<AliAODEvent*> (InputEvent());
2082   if(!event) return;
2083   AliVVZERO *vzero = event->GetVZEROData();
2084   AliAODTracklets *tracklets = event->GetTracklets();
2085   if(!vzero) return;
2086   if(!tracklets) return;
2087   Double_t mult, eqmult;
2088   Int_t trkl;
2089   for(int id=0; id!=64; ++id) {
2090     trkl = tracklets->GetNumberOfTracklets();
2091     mult = vzero->GetMultiplicity(id);
2092     eqmult = event->GetVZEROEqMultiplicity(id);
2093     ((TProfile2D*) fVZEQA->FindObject( "LINP" ))->Fill(id,mult,trkl,1);
2094     ((TProfile2D*) fVZEQA->FindObject( "MULP" ))->Fill(id,fThisCent,mult,1);
2095     ((TH3D*) fVZEQA->FindObject("EQU"))->Fill(eqmult,mult,id);
2096   }
2097 }
2098 //=======================================================================
2099 void AliAnalysisTaskFlowStrange::AddVZEROResponse() {
2100   fVZEResponse = NULL;
2101   AliVEvent *event = InputEvent();
2102   if(!event) return;
2103   Int_t thisrun = event->GetRunNumber();
2104   fVZEResponse = new TH2D( Form("%d",thisrun), Form("%d;cell;CC",thisrun), 64,0,64, 100, 0, 100);
2105   fList->Add(fVZEResponse);
2106 }
2107 //=======================================================================
2108 void AliAnalysisTaskFlowStrange::SaveVZEROResponse() {
2109   if(!fVZEResponse) return;
2110   AliVEvent *event = InputEvent();
2111   if(!event) return;
2112   Double_t w;
2113   for(int id=0; id!=64; ++id) {
2114     w = event->GetVZEROEqMultiplicity(id);
2115     fVZEResponse->Fill(id,fThisCent,w);
2116   }
2117 }
2118 //=======================================================================
2119 Int_t AliAnalysisTaskFlowStrange::RefMultTPC() {
2120   AliAODEvent *ev = dynamic_cast<AliAODEvent*> (InputEvent());
2121   if(!ev) return -1;
2122   Int_t found = 0;
2123   for(int i=0; i!=ev->GetNumberOfTracks(); ++i) {
2124     AliAODTrack *t = ev->GetTrack( i );
2125     if(!t) continue;
2126     if( !t->TestFilterBit(1) ) continue;
2127     if( t->Eta()<-0.8 || t->Eta()>+0.8 ) continue;
2128     if( t->Pt()<0.2 || t->Pt()>5.0 ) continue;
2129     if( t->GetTPCNcls()<70 ) continue;
2130     if( t->GetTPCsignal()<10.0 ) continue;
2131     if( t->Chi2perNDF()<0.2 ) continue;
2132     ++found;
2133   }
2134   return found;
2135 }
2136 //=======================================================================
2137 Int_t AliAnalysisTaskFlowStrange::RefMultGlobal() {
2138   AliAODEvent *ev = dynamic_cast<AliAODEvent*> (InputEvent());
2139   if(!ev) return -1;
2140   Int_t found = 0;
2141   for(int i=0; i!=ev->GetNumberOfTracks(); ++i) {
2142     AliAODTrack *t = ev->GetTrack( i );
2143     if(!t) continue;
2144     if( !t->TestFilterBit(16) ) continue;
2145     if( t->Eta()<-0.8 || t->Eta()>+0.8 ) continue;
2146     if( t->Pt()<0.2 || t->Pt()>5.0 ) continue;
2147     if( t->GetTPCNcls()<70 ) continue;
2148     if( t->GetTPCsignal()<10.0 ) continue;
2149     if( t->Chi2perNDF()<0.1 ) continue;
2150     Double_t b[3], bcov[3];
2151     if( !t->PropagateToDCA(ev->GetPrimaryVertex(),ev->GetMagneticField(),100,b,bcov) ) continue;
2152     if( b[0]>+0.3 || b[0]<-0.3 || b[1]>+0.3 || b[1]<-0.3) continue;
2153     ++found;
2154   }
2155   return found;
2156 }
2157 //=======================================================================
2158 void AliAnalysisTaskFlowStrange::MakeDHcorr() {
2159   // Adds DH corr
2160   for(int iCand=0; iCand!=fCandidates->GetEntriesFast(); ++iCand ) {
2161     AliFlowCandidateTrack *cand = static_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
2162     if(!cand) continue;
2163     for(int iRPs=0; iRPs!=fTPCevent->NumberOfTracks(); ++iRPs ) {
2164       AliFlowTrack *iRP = static_cast<AliFlowTrack*>(fTPCevent->GetTrack( iRPs ));
2165       if(!iRP) continue;
2166       if(!iRP->InRPSelection()) continue;
2167       if(cand->GetID() == iRP->GetID()) continue; //avoid autocorr
2168       for(int iDau=0; iDau!=cand->GetNDaughters(); ++iDau) //if it is a decay
2169         if(cand->GetIDDaughter(iDau) == iRP->GetID()) continue;
2170       //corr
2171       Double_t dDPHI = iRP->Phi() - cand->Phi();
2172       Double_t dDETA = iRP->Eta() - cand->Eta();
2173       Double_t dDPT  = iRP->Pt()  - cand->Pt();
2174       ((TH3D*)((TList*)fList->FindObject("DHCORR"))->FindObject("DPHI"))->Fill( dDPT, dDPHI, dDETA );
2175       //end of corr
2176     }
2177   }
2178 }
2179 //=======================================================================
2180 void AliAnalysisTaskFlowStrange::ResetContainers() {
2181   fTPCevent->ClearFast();
2182   fVZEevent->ClearFast();
2183 }