PWG/FLOW/Tasks/AliAnalysisTaskFlowStrange.cxx

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