AliAODEvent::GetHeader() returns AliVHeader
[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliAnalysisTaskFlowStrange.cxx
1 /*************************************************************************
2 * Copyright(c) 1998-2008,ALICE Experiment at CERN,All rights reserved. *
3 *                                                                        *
4 * Author: The ALICE Off-line Project.                                    *
5 * Contributors are mentioned in the code where appropriate.              *
6 *                                                                        *
7 * Permission to use,copy,modify and distribute this software and its   *
8 * documentation strictly for non-commercial purposes is hereby granted   *
9 * without fee,provided that the above copyright notice appears in all   *
10 * copies and that both the copyright notice and this permission notice   *
11 * appear in the supporting documentation. The authors make no claims     *
12 * about the suitability of this software for any purpose. It is          *
13 * provided "as is" without express or implied warranty.                  *
14 **************************************************************************/
15
16 /////////////////////////////////////////////////////
17 // AliAnalysisTaskFlowStrange:
18 // Analysis task to select K0/Lambda candidates for flow analysis.
19 // Authors: Cristian Ivan (civan@cern.ch)
20 //          Carlos Perez  (cperez@cern.ch)
21 //          Pawel Debski  (pdebski@cern.ch)
22 //////////////////////////////////////////////////////
23
24 #include "TChain.h"
25 #include "TList.h"
26 #include "TH1D.h"
27 #include "TH2D.h"
28 #include "TH3D.h"
29 #include "TF1.h"
30 #include "TProfile.h"
31 #include "TProfile2D.h"
32 #include "TVector3.h"
33 #include "TStopwatch.h"
34 #include "TFile.h"
35
36 #include "TRandom3.h"
37
38 #include "AliAnalysisManager.h"
39 #include "AliInputEventHandler.h"
40
41 #include "AliVVertex.h"
42 #include "AliVVZERO.h"
43 #include "AliStack.h"
44 #include "AliMCEvent.h"
45
46 #include "AliESDEvent.h"
47 #include "AliESDtrack.h"
48 #include "AliESDVertex.h"
49 #include "AliESDv0.h"
50 #include "AliESDtrackCuts.h"
51
52 #include "AliAODEvent.h"
53 #include "AliAODTrack.h"
54 #include "AliAODVertex.h"
55 #include "AliAODv0.h"
56 #include "AliAODTracklets.h"
57 #include "AliAODHeader.h"
58
59 #include "AliAODMCHeader.h"
60 #include "AliAODMCParticle.h"
61 #include "TClonesArray.h"
62 #include "TDatabasePDG.h"
63 #include "TParticlePDG.h"
64
65 #include "TMath.h"
66 #include "TObjArray.h"
67 #include "AliFlowCandidateTrack.h"
68
69 #include "AliFlowTrackCuts.h"
70 #include "AliFlowEventCuts.h"
71 #include "AliFlowEvent.h"
72 #include "AliFlowBayesianPID.h"
73 #include "AliFlowCommonConstants.h"
74 #include "AliFlowVector.h"
75
76 #include "AliAnalysisTaskFlowStrange.h"
77
78 ClassImp(AliAnalysisTaskFlowStrange)
79
80 //=======================================================================
81 AliAnalysisTaskFlowStrange::AliAnalysisTaskFlowStrange() :
82   AliAnalysisTaskSE(),
83   fPIDResponse(NULL),
84   fFB1(NULL),
85   fFB1024(NULL),
86   fTPCevent(NULL),
87   fVZEevent(NULL),
88   fCandidates(NULL),
89   fList(NULL),
90   fRunNumber(-1),
91   fDebug(0),
92   fQAlevel(0),
93   fReadESD(kFALSE),
94   fReadMC(kFALSE),
95   fAddPiToMCReactionPlane(kTRUE),
96   fPostMatched(0),
97   fAvoidExec(kFALSE),
98   fSkipSelection(kFALSE),
99   fSkipVn(kFALSE),
100   fUseFP(kFALSE),
101   fRunOnpA(kFALSE),
102   fRunOnpp(kFALSE),
103   fExtraEventRejection(kFALSE),
104   fSkipCentralitySelection(kFALSE),
105   fCentMethod("V0MTRK"),
106   fCentPerMin(0),
107   fCentPerMax(100),
108   fThisCent(-1.0),
109   fV0M(0.0),
110   fTRK(0.0),
111   fPriVtxZ(0.0),
112   fSPDVtxZ(0.0),
113   fSPDtracklets(0),
114   fVZETotM(0.0),
115   fRefMultTPC(0),
116   fRefMultHyb(0),
117   fVertexZcut(10.0),
118   fExcludeTPCEdges(kFALSE),
119   fSpecie(0),
120   fOnline(kFALSE),
121   fHomemade(kFALSE),
122   fWhichPsi(1),
123   fVZEsave(kFALSE),
124   fVZEload(NULL),
125   fVZEResponse(NULL),
126   fVZEmb(kFALSE),
127   fVZEByDisk(kTRUE),
128   fVZECa(0),
129   fVZECb(3),
130   fVZEAa(0),
131   fVZEAb(3),
132   fVZEQA(NULL),
133   fHarmonic(2),
134   fPsi2(0.0),
135   fMCEP(0.0),
136   fQVZEACos(0.0),
137   fQVZEASin(0.0),
138   fQVZECCos(0.0),
139   fQVZECSin(0.0),
140   fQVZEA(0.0),
141   fQVZEC(0.0),
142   fVZEWarning(kFALSE),
143   fQTPCACos(0.0),
144   fQTPCASin(0.0),
145   fQTPCCCos(0.0),
146   fQTPCCSin(0.0),
147   fQTPC2hCos(0.0),
148   fQTPC2hSin(0.0),
149   fQTPCA(0.0),
150   fQTPCC(0.0),
151   fQTPCA_nTracks(0),
152   fQTPCC_nTracks(0),
153   fSkipTerminate(kTRUE),
154   fMassBins(0),
155   fMinMass(0.0),
156   fMaxMass(0.0),
157   fPtBins(0),
158   fRFPFilterBit(1),
159   fRFPminPt(0.2),
160   fRFPmaxPt(5.0),
161   fRFPAminEta(0.0),
162   fRFPAmaxEta(+0.8),
163   fRFPCminEta(-0.8),
164   fRFPCmaxEta(0.0),
165   fRFPTPCsignal(10.0),
166   fRFPmaxIPxy(2.4),
167   fRFPmaxIPz(3.2),
168   fRFPTPCncls(70),
169   fDecayMass(0.0),
170   fDecayPhi(0.0),
171   fDecayEta(0.0),
172   fDecayPt(0.0),
173   fDecayDCAdaughters(0.0),
174   fDecayCosinePointingAngleXY(0.0),
175   fDecayRadXY(0.0),
176   fDecayDecayLength(0.0),
177   fDecayDecayLengthLab(0.0),
178   fDecayQt(0.0),
179   fDecayAlpha(0.0),
180   fDecayRapidity(0.0),
181   fDecayProductIPXY(0.0),
182   fDecayIPneg(0.0),
183   fDecayIPpos(0.0),
184   fDecayXneg(0.0),
185   fDecayXpos(0.0),
186   fDecayIDneg(-1),
187   fDecayIDpos(-1),
188   fDecayID(-1),
189   fDecayMatchOrigin(0.0),
190   fDecayMatchPhi(0.0),
191   fDecayMatchEta(0.0),
192   fDecayMatchPt(0.0),
193   fDecayMatchRadXY(0.0),
194   fDecayMinEta(0.0),
195   fDecayMaxEta(0.0),
196   fDecayMinPt(0.0),
197   fDecayMaxDCAdaughters(0.0),
198   fDecayMinCosinePointingAngleXY(0.0),
199   fDecayMinQt(0.0),
200   fDecayAPCutPie(kTRUE),
201   fDecayStopPIDAtPt(3.0),
202   fDecayMinRadXY(0.0),
203   fDecayMaxDecayLength(0.0),
204   fDecayMaxProductIPXY(0.0),
205   fDecayMaxRapidity(0.0),
206   fDaughterPhi(0.0),
207   fDaughterEta(0.0),
208   fDaughterPt(0.0),
209   fDaughterNClsTPC(0),
210   fDaughterNClsITS(0),
211   fDaughterCharge(0),
212   fDaughterNFClsTPC(0),
213   fDaughterNSClsTPC(0),
214   fDaughterChi2PerNClsTPC(0.0),
215   fDaughterXRows(0.0),
216   fDaughterImpactParameterXY(0.0),
217   fDaughterImpactParameterZ(0.0),
218   fDaughterStatus(0),
219   fDaughterITScm(0),
220   fDaughterNSigmaPID(0.0),
221   fDaughterKinkIndex(0),
222   fDaughterAtSecPhi(0.0),
223   fDaughterAtSecEta(0.0),
224   fDaughterAtSecPt(0.0),
225   fDaughterMatchPhi(0.0),
226   fDaughterMatchEta(0.0),
227   fDaughterMatchPt(0.0),
228   fDaughterMatchImpactParameterXY(0.0),
229   fDaughterMatchImpactParameterZ(0.0),
230   fDaughterUnTag(kTRUE),
231   fDaughterMinEta(0.0),
232   fDaughterMaxEta(0.0),
233   fDaughterMinPt(0.0),
234   fDaughterMinNClsTPC(0),
235   fDaughterMinNClsITS(-1),
236   fDaughterMinXRows(0),
237   fDaughterMaxChi2PerNClsTPC(0.0),
238   fDaughterMinXRowsOverNClsFTPC(0.0),
239   fDaughterMinImpactParameterXY(0.0),
240   fDaughterMaxNSigmaPID(0.0),
241   fDaughterSPDRequireAny(kFALSE),
242   fDaughterITSrefit(kFALSE) {
243   //ctor
244   for(Int_t i=0; i!=100; ++i) fPtBinEdge[i]=0;
245   for(Int_t i=0; i!=6; ++i) fDaughterITSConfig[i]=-1;
246   for(Int_t i=0; i!=2000; ++i) fQTPCA_fID[i]=-1;
247   for(Int_t i=0; i!=2000; ++i) fQTPCC_fID[i]=-1;
248   for(Int_t i=0; i!=64; ++i) fVZEextW[i]=1;
249 }
250 //=======================================================================
251 AliAnalysisTaskFlowStrange::AliAnalysisTaskFlowStrange(const char *name) :
252   AliAnalysisTaskSE(name),
253   fPIDResponse(NULL),
254   fFB1(NULL),
255   fFB1024(NULL),
256   fTPCevent(NULL),
257   fVZEevent(NULL),
258   fCandidates(NULL),
259   fList(NULL),
260   fRunNumber(-1),
261   fDebug(0),
262   fQAlevel(0),
263   fReadESD(kFALSE),
264   fReadMC(kFALSE),
265   fAddPiToMCReactionPlane(kTRUE),
266   fPostMatched(0),
267   fAvoidExec(kFALSE),
268   fSkipSelection(kFALSE),
269   fSkipVn(kFALSE),
270   fUseFP(kFALSE),
271   fRunOnpA(kFALSE),
272   fRunOnpp(kFALSE),
273   fExtraEventRejection(kFALSE),
274   fSkipCentralitySelection(kFALSE),
275   fCentMethod("V0MTRK"),
276   fCentPerMin(0),
277   fCentPerMax(100),
278   fThisCent(-1.0),
279   fV0M(0.0),
280   fTRK(0.0),
281   fPriVtxZ(0.0),
282   fSPDVtxZ(0.0),
283   fSPDtracklets(0),
284   fVZETotM(0.0),
285   fRefMultTPC(0),
286   fRefMultHyb(0),
287   fVertexZcut(10.0),
288   fExcludeTPCEdges(kFALSE),
289   fSpecie(0),
290   fOnline(kFALSE),
291   fHomemade(kFALSE),
292   fWhichPsi(1),
293   fVZEsave(kFALSE),
294   fVZEload(NULL),
295   fVZEResponse(NULL),
296   fVZEmb(kFALSE),
297   fVZEByDisk(kTRUE),
298   fVZECa(0),
299   fVZECb(3),
300   fVZEAa(0),
301   fVZEAb(3),
302   fVZEQA(NULL),
303   fHarmonic(2),
304   fPsi2(0.0),
305   fMCEP(0.0),
306   fQVZEACos(0.0),
307   fQVZEASin(0.0),
308   fQVZECCos(0.0),
309   fQVZECSin(0.0),
310   fQVZEA(0.0),
311   fQVZEC(0.0),
312   fVZEWarning(kFALSE),
313   fQTPCACos(0.0),
314   fQTPCASin(0.0),
315   fQTPCCCos(0.0),
316   fQTPCCSin(0.0),
317   fQTPC2hCos(0.0),
318   fQTPC2hSin(0.0),
319   fQTPCA(0.0),
320   fQTPCC(0.0),
321   fQTPCA_nTracks(0),
322   fQTPCC_nTracks(0),
323   fSkipTerminate(kTRUE),
324   fMassBins(0),
325   fMinMass(0.0),
326   fMaxMass(0.0),
327   fPtBins(0),
328   fRFPFilterBit(1),
329   fRFPminPt(0.2),
330   fRFPmaxPt(5.0),
331   fRFPAminEta(0.0),
332   fRFPAmaxEta(+0.8),
333   fRFPCminEta(-0.8),
334   fRFPCmaxEta(0.0),
335   fRFPTPCsignal(10.0),
336   fRFPmaxIPxy(2.4),
337   fRFPmaxIPz(3.2),
338   fRFPTPCncls(70),
339   fDecayMass(0.0),
340   fDecayPhi(0.0),
341   fDecayEta(0.0),
342   fDecayPt(0.0),
343   fDecayDCAdaughters(0.0),
344   fDecayCosinePointingAngleXY(0.0),
345   fDecayRadXY(0.0),
346   fDecayDecayLength(0.0),
347   fDecayDecayLengthLab(0.0),
348   fDecayQt(0.0),
349   fDecayAlpha(0.0),
350   fDecayRapidity(0.0),
351   fDecayProductIPXY(0.0),
352   fDecayIPneg(0.0),
353   fDecayIPpos(0.0),
354   fDecayXneg(0.0),
355   fDecayXpos(0.0),
356   fDecayIDneg(-1),
357   fDecayIDpos(-1),
358   fDecayID(-1),
359   fDecayMatchOrigin(0.0),
360   fDecayMatchPhi(0.0),
361   fDecayMatchEta(0.0),
362   fDecayMatchPt(0.0),
363   fDecayMatchRadXY(0.0),
364   fDecayMinEta(0.0),
365   fDecayMaxEta(0.0),
366   fDecayMinPt(0.0),
367   fDecayMaxDCAdaughters(0.0),
368   fDecayMinCosinePointingAngleXY(0.0),
369   fDecayMinQt(0.0),
370   fDecayAPCutPie(kTRUE),
371   fDecayStopPIDAtPt(3.0),
372   fDecayMinRadXY(0.0),
373   fDecayMaxDecayLength(0.0),
374   fDecayMaxProductIPXY(0.0),
375   fDecayMaxRapidity(0.0),
376   fDaughterPhi(0.0),
377   fDaughterEta(0.0),
378   fDaughterPt(0.0),
379   fDaughterNClsTPC(0),
380   fDaughterNClsITS(0),
381   fDaughterCharge(0),
382   fDaughterNFClsTPC(0),
383   fDaughterNSClsTPC(0),
384   fDaughterChi2PerNClsTPC(0.0),
385   fDaughterXRows(0.0),
386   fDaughterImpactParameterXY(0.0),
387   fDaughterImpactParameterZ(0.0),
388   fDaughterStatus(0),
389   fDaughterITScm(0),
390   fDaughterNSigmaPID(0.0),
391   fDaughterKinkIndex(0),
392   fDaughterAtSecPhi(0.0),
393   fDaughterAtSecEta(0.0),
394   fDaughterAtSecPt(0.0),
395   fDaughterMatchPhi(0.0),
396   fDaughterMatchEta(0.0),
397   fDaughterMatchPt(0.0),
398   fDaughterMatchImpactParameterXY(0.0),
399   fDaughterMatchImpactParameterZ(0.0),
400   fDaughterUnTag(kTRUE),
401   fDaughterMinEta(0.0),
402   fDaughterMaxEta(0.0),
403   fDaughterMinPt(0.0),
404   fDaughterMinNClsTPC(0),
405   fDaughterMinNClsITS(-1),
406   fDaughterMinXRows(0),
407   fDaughterMaxChi2PerNClsTPC(0.0),
408   fDaughterMinXRowsOverNClsFTPC(0.0),
409   fDaughterMinImpactParameterXY(0.0),
410   fDaughterMaxNSigmaPID(0.0),
411   fDaughterSPDRequireAny(kFALSE),
412   fDaughterITSrefit(kFALSE) {
413   //ctor
414   for(Int_t i=0; i!=100; ++i) fPtBinEdge[i]=0;
415   for(Int_t i=0; i!=6; ++i) fDaughterITSConfig[i]=-1;
416   for(Int_t i=0; i!=2000; ++i) fQTPCA_fID[i]=-1;
417   for(Int_t i=0; i!=2000; ++i) fQTPCC_fID[i]=-1;
418   for(Int_t i=0; i!=64; ++i) fVZEextW[i]=1;
419   DefineInput( 0,TChain::Class());
420   DefineOutput(1,TList::Class());
421   DefineOutput(2,AliFlowEventSimple::Class()); // TPC object
422   DefineOutput(3,AliFlowEventSimple::Class()); // VZE object
423 }
424 //=======================================================================
425 AliAnalysisTaskFlowStrange::~AliAnalysisTaskFlowStrange() {
426   //dtor
427   if (fCandidates) delete fCandidates;
428   if (fTPCevent)   delete fTPCevent;
429   if (fVZEevent)   delete fVZEevent;
430   if (fList)       delete fList;
431 }
432 //=======================================================================
433 void AliAnalysisTaskFlowStrange::SetPtEdges(Int_t n, Double_t *p) {
434   fPtBins = n;
435   for(int i=0;i!=n+1;++i) fPtBinEdge[i] = p[i];
436 }
437 //=======================================================================
438 TList* AliAnalysisTaskFlowStrange::RunTerminateAgain(TList *lst) {
439   if(!lst) return NULL;
440   fList = lst;
441   fSpecie = Int_t( ((TProfile*)((TList*)fList->FindObject("Event"))->FindObject("Configuration"))->GetBinContent(kSpecie) );
442   fSkipSelection = ((TProfile*)((TList*)fList->FindObject("Event"))->FindObject("Configuration"))->GetBinContent(kSkipSelection);
443   fReadMC = ((TProfile*)((TList*)fList->FindObject("Event"))->FindObject("Configuration"))->GetBinContent(kReadMC);
444   Terminate(NULL);
445   return fList;
446 }
447 //=======================================================================
448 void AliAnalysisTaskFlowStrange::PrintConfig() {
449   //DUMP for main task
450   printf("******************************\n");
451   printf("<TASK Configuration> %s\n",GetName());
452   printf("  fDebug %d\n",fDebug);
453   printf("  fQAlevel %d\n",fQAlevel);
454   printf("  fExtraEventRejection %s\n",fExtraEventRejection?"kTRUE":"kFALSE");
455   printf("  fCentMethod %s\n",fCentMethod.Data());
456   printf("    fCentPerMin %d\n",fCentPerMin);
457   printf("    fCentPerMax %d\n",fCentPerMax);
458   printf("  fVextexZcut %f\n",fVertexZcut);
459   printf("  fRunOnpA %s\n",fRunOnpA?"kTRUE":"kFALSE");
460   printf("  fRunOnpp %s\n",fRunOnpp?"kTRUE":"kFALSE");
461   printf("  fReadESD %s\n",fReadESD?"kTRUE":"kFALSE");
462   printf("  fReadMC %s\n",fReadMC?"kTRUE":"kFALSE");
463   if(fReadMC) {
464     printf("    fAddPiToMCReactionPlane %s\n",fAddPiToMCReactionPlane?"kTRUE":"kFALSE");
465     printf("    fPostMatched %d\n",fPostMatched);
466     printf("    fAvoidExec %s\n",fAvoidExec?"kTRUE":"kFALSE");
467     printf("    fSkipCentralitySelection %s\n",fSkipCentralitySelection?"kTRUE":"kFALSE");
468   }
469   printf("  fVZEsave %s\n",fVZEsave?"kTRUE":"kFALSE");
470   if(fVZEload) {
471     printf("  fVZEload %d runs\n",fVZEload->GetEntries());
472     printf("    fVZEmb %s\n",fVZEmb?"kTRUE":"kFALSE");
473     printf("    fVZEByDisk %s\n",fVZEByDisk?"kTRUE":"kFALSE");
474   }
475   printf("  fHarmonic %d\n",fHarmonic);
476   printf("    fWhichPsi %d\n",fWhichPsi);
477   printf("    fVZECa %d\n",fVZECa);
478   printf("    fVZECb %d\n",fVZECb);
479   printf("    fVZEAa %d\n",fVZEAa);
480   printf("    fVZEAb %d\n",fVZEAb);
481   printf("    fRFPFilterBit %d\n",fRFPFilterBit);
482   printf("    fRFPminPt %f\n",fRFPminPt);
483   printf("    fRFPmaxPt %f\n",fRFPmaxPt);
484   printf("    fRFPAminEta %f\n",fRFPAminEta);
485   printf("    fRFPAmaxEta %f\n",fRFPAmaxEta);
486   printf("    fRFPCminEta %f\n",fRFPCminEta);
487   printf("    fRFPCmaxEta %f\n",fRFPCmaxEta);
488   printf("    fRFPmaxIPxy %f\n",fRFPmaxIPxy);
489   printf("    fRFPmaxIPz %f\n",fRFPmaxIPz);
490   printf("    fRFPTPCsignal %f\n",fRFPTPCsignal);
491   printf("    fRFPTPCncls %d\n",fRFPTPCncls);
492   printf("    fExcludeTPCEdges %s\n",fExcludeTPCEdges?"kTRUE":"kFALSE");
493   printf("  fSkipSelection %s\n",fSkipSelection?"kTRUE":"kFALSE");
494   if(!fSkipSelection) {
495     printf("    fSpecie %d\n",fSpecie);
496     printf("    fPtBins %d\n      |",fPtBins);
497     for(int i=0; i!=fPtBins+1; ++i) printf("%f|",fPtBinEdge[i]); printf("\n");
498     if(fSpecie<90) {
499       printf("    fMassBins %d\n",fMassBins);
500       printf("    fMinMass %f\n",fMinMass);
501       printf("    fMaxMass %f\n",fMaxMass);
502     }
503   }
504   printf("  fSkipVn %s\n",fSkipVn?"kTRUE":"kFALSE");
505   if(!fSkipVn) {
506     printf("    fUseFP %s\n",fUseFP?"kTRUE":"kFALSE");
507   }
508   MyPrintConfig();
509 }
510 //=======================================================================
511 void AliAnalysisTaskFlowStrange::MyPrintConfig() {
512   // Dump for derived task
513   printf("==================================\n");
514   printf("<FlowStrange> \n");
515   if(!fSkipSelection) {
516     if(fReadESD) {
517       printf("  fOnline %s\n",fOnline?"kTRUE":"kFALSE");
518       printf("  fHomemade %s\n",fHomemade?"kTRUE":"kFALSE");
519     }
520     printf("  fDecayMinEta %f\n",fDecayMinEta);
521     printf("  fDecayMaxEta %f\n",fDecayMaxEta);
522     printf("  fDecayMinPt %f\n",fDecayMinPt);
523     printf("  fDecayMaxDCAdaughters %f\n",fDecayMaxDCAdaughters);
524     printf("  fDecayMinCosinePointingAngleXY %f\n",fDecayMinCosinePointingAngleXY);
525     printf("  fDecayMinQt %f\n",fDecayMinQt);
526     printf("  fDecayAPCutPie %s\n",fDecayAPCutPie?"kTRUE":"kFALSE");
527     printf("  fDecayStopPIDAtPt %f\n",fDecayStopPIDAtPt);
528     printf("  fDecayMinRadXY %f\n",fDecayMinRadXY);
529     printf("  fDecayMaxDecayLength %f\n",fDecayMaxDecayLength);
530     printf("  fDecayMaxProductIPXY %f\n",fDecayMaxProductIPXY);
531     printf("  fDecayMaxRapidity %f\n",fDecayMaxRapidity);
532   }
533   printf("  fDaughterUnTag %s\n",fDaughterUnTag?"kTRUE":"kFALSE");
534   printf("  fDaughterMinEta %f\n",fDaughterMinEta);
535   printf("  fDaughterMaxEta %f\n",fDaughterMaxEta);
536   printf("  fDaughterMinPt %f\n",fDaughterMinPt);
537   printf("  fDaughterMinNClsTPC %d\n",fDaughterMinNClsTPC);
538   printf("  fDaughterMinXRows %d\n",fDaughterMinXRows);
539   printf("  fDaughterMaxChi2PerNClsTPC %f\n",fDaughterMaxChi2PerNClsTPC);
540   printf("  fDaughterMinXRowsOverNClsFTPC %f\n",fDaughterMinXRowsOverNClsFTPC);
541   printf("  fDaughterMinImpactParameterXY %f\n",fDaughterMinImpactParameterXY);
542   printf("  fDaughterMaxNSigmaPID %f\n",fDaughterMaxNSigmaPID);
543 }
544 //=======================================================================
545 void AliAnalysisTaskFlowStrange::UserCreateOutputObjects() {
546   //UserCreateOutputObjects
547   if(fDebug) PrintConfig();
548   fList=new TList();
549   fList->SetOwner();
550   AddQAEvents();
551   AddQACandidates();
552   if(fReadESD) MakeFilterBits();
553
554   AliFlowCommonConstants *cc = AliFlowCommonConstants::GetMaster();
555   cc->SetNbinsMult(3000); cc->SetMultMin(0);   cc->SetMultMax(30000);
556   cc->SetNbinsPt(100); cc->SetPtMin(0.0);   cc->SetPtMax(20.0);
557   cc->SetNbinsPhi(100);  cc->SetPhiMin(0.0);  cc->SetPhiMax(TMath::TwoPi());
558   cc->SetNbinsEta(100);  cc->SetEtaMin(-5.0); cc->SetEtaMax(+5.0);
559   cc->SetNbinsQ(100);    cc->SetQMin(0.0);    cc->SetQMax(3.0);
560   cc->SetNbinsMass(fMassBins);
561   cc->SetMassMin(fMinMass);
562   cc->SetMassMax(fMaxMass);
563
564   //loading pid response
565   AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
566   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
567   fPIDResponse = inputHandler->GetPIDResponse();
568
569   if(fUseFP) {
570     fTPCevent = new AliFlowEvent(100);
571     fVZEevent = new AliFlowEvent(100);
572     //array of candidates
573     fCandidates = new TObjArray(100);
574     fCandidates->SetOwner();
575   }
576   PostData(1,fList);
577   if(fUseFP) { // for connection to the flow package
578     PostData(2,fTPCevent);
579     PostData(3,fVZEevent);
580   }
581
582   gRandom->SetSeed();
583 }
584 //=======================================================================
585 void AliAnalysisTaskFlowStrange::MyUserCreateOutputObjects() {
586   TList *tList;
587   TH1D *tH1D;
588   TH2D *tH2D;
589
590   //reconstruction
591   if(fReadESD) {
592     tList=new TList(); tList->SetName("ESD_TrkAll"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
593     tList=new TList(); tList->SetName("ESD_TrkSel"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
594     tH2D = new TH2D("NPAIR", "NPAIR;NPOS;NNEG",1000,0,5000,1000,0,5000); tList->Add(tH2D);
595     tH2D = new TH2D("PtIPXY","PtIPXY;Pt;IPxy", 100,0,10,200,-10,+10); tList->Add(tH2D);
596   }
597   //aod prefilter candidates
598   tList=new TList(); tList->SetName("V0SAll"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
599   tH2D = new TH2D("V0SADC","V0S AFTER DAUGHTER CUTS;V0ALL;V0IMW",100,0,1000,100,0,1000); tList->Add(tH2D);
600   tList=new TList(); tList->SetName("AllDau"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
601   //candidates
602   tList=new TList(); tList->SetName("V0SSel"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
603   tList=new TList(); tList->SetName("SelDau"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
604   //flow
605   if(!fSkipVn) {
606     tList=new TList(); tList->SetName("V0SAllVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
607     tList=new TList(); tList->SetName("V0SSelVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
608   }
609   // IN-OUT
610   if(fQAlevel>1) {
611     tList=new TList(); tList->SetName("V0SAllIP"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
612     tList=new TList(); tList->SetName("V0SAllOP"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
613     tList=new TList(); tList->SetName("V0SSelIP"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
614     tList=new TList(); tList->SetName("V0SSelOP"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
615   }
616   //match
617   if(fReadMC) {
618     tList=new TList(); tList->SetName("STATMC"); tList->SetOwner(); fList->Add(tList);
619     tH1D = new TH1D("Events", "Events",5,0.5,5.5); tList->Add(tH1D);
620     tH1D->GetXaxis()->SetBinLabel(1,"Selected events");
621     tH1D->GetXaxis()->SetBinLabel(2,"Stack found");
622     tH1D->GetXaxis()->SetBinLabel(3,"Daughters in stack");
623     tH1D->GetXaxis()->SetBinLabel(4,"Correspond to decay");
624     tH1D->GetXaxis()->SetBinLabel(5,"Decay has mother");
625     tList=new TList(); tList->SetName("Mth"); tList->SetOwner(); AddCandidatesSpy(tList,true); fList->Add(tList);
626     tList=new TList(); tList->SetName("MthDau"); tList->SetOwner(); AddTrackSpy(tList,true); fList->Add(tList);
627     tList=new TList(); tList->SetName("MthPosPos"); tList->SetOwner(); AddCandidatesSpy(tList,true); fList->Add(tList);
628     tList=new TList(); tList->SetName("MthNegNeg"); tList->SetOwner(); AddCandidatesSpy(tList,true); fList->Add(tList);
629     tList=new TList(); tList->SetName("MthPosNeg"); tList->SetOwner(); AddCandidatesSpy(tList,true); fList->Add(tList);
630     tList=new TList(); tList->SetName("MthNegDau"); tList->SetOwner(); AddTrackSpy(tList,true); fList->Add(tList);
631     tList=new TList(); tList->SetName("MthPosDau"); tList->SetOwner(); AddTrackSpy(tList,true); fList->Add(tList);
632     tList=new TList(); tList->SetName("MthFeedDown"); tList->SetOwner(); AddCandidatesSpy(tList,true); fList->Add(tList);
633     tList=new TList(); tList->SetName("UnMth"); tList->SetOwner(); AddCandidatesSpy(tList,false); fList->Add(tList);
634     tList=new TList(); tList->SetName("UnMthDau"); tList->SetOwner(); AddTrackSpy(tList,false); fList->Add(tList);
635     if(!fSkipVn) {
636       tList=new TList(); tList->SetName("V0SMthVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
637       tList=new TList(); tList->SetName("V0SMthPosPosVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
638       tList=new TList(); tList->SetName("V0SMthNegNegVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
639       tList=new TList(); tList->SetName("V0SMthPosNegVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
640       tList=new TList(); tList->SetName("V0SUnMthVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
641     }
642   }
643 }
644 //=======================================================================
645 void AliAnalysisTaskFlowStrange::AddQAEvents() {
646   // function to add event qa
647   TH1D *tH1D;
648   TProfile *tProfile;
649   TList *tQAEvents=new TList();
650   tQAEvents->SetName("Event");
651   tQAEvents->SetOwner();
652   tH1D = new TH1D("Events","Number of Events",6,0,6); tQAEvents->Add(tH1D);
653   tH1D->GetXaxis()->SetBinLabel(1,"exec");
654   tH1D->GetXaxis()->SetBinLabel(2,"userexec");
655   tH1D->GetXaxis()->SetBinLabel(3,"reached");
656   tH1D->GetXaxis()->SetBinLabel(4,"selected");
657   tH1D->GetXaxis()->SetBinLabel(5,"rejectedByLowQw");
658   tH1D->GetXaxis()->SetBinLabel(6,"rejectedByErrorLoadVZEcal");
659   tProfile = new TProfile("Configuration","Configuration",10,0.5,10.5); tQAEvents->Add(tProfile);
660   tProfile->Fill(kSpecie,fSpecie,1);
661   tProfile->GetXaxis()->SetBinLabel(kSpecie,"fSpecie");
662   tProfile->Fill(kHarmonic,fHarmonic,1);
663   tProfile->GetXaxis()->SetBinLabel(kHarmonic,"fHarmonic");
664   tProfile->Fill(kReadMC,fReadMC,1);
665   tProfile->GetXaxis()->SetBinLabel(kReadMC,"fReadMC");
666   tProfile->Fill(kSkipSelection,fSkipSelection,1);
667   tProfile->GetXaxis()->SetBinLabel(kSkipSelection,"fSkipSelection");
668   tH1D = new TH1D("POI","POIs;multiplicity",800,0,800);         tQAEvents->Add(tH1D);
669   tH1D = new TH1D("UNTAG","UNTAG;Untagged Daughters",800,0,800);tQAEvents->Add(tH1D);
670   tH1D = new TH1D("RealTime","RealTime;LogT sec",2000,-10,+10); tQAEvents->Add(tH1D);
671   fList->Add(tQAEvents);
672   AddEventSpy("EventsRaw");
673   AddEventSpy("EventsReached");
674   AddEventSpy("EventsSelected");
675   AddEventSpy("EventsAnalyzed");
676   AddMakeQSpy();
677   AddVZEQA();
678 }
679 //=======================================================================
680 void AliAnalysisTaskFlowStrange::AddEventSpy(TString name) {
681   TH1D *tH1D;
682   TH2D *tH2D;
683   TList *tList=new TList();
684   tList->SetName(name.Data());
685   tList->SetOwner();
686   tH2D = new TH2D("VTXZ","VTXZ;PriVtxZ;SPDVtxZ",60,-25,+25,60,-25,+25); tList->Add( tH2D );
687   tH2D = new TH2D("CCCC","CCCC;V0M;TRK",60,-10,110,60,-10,110);         tList->Add( tH2D );
688   tH2D = new TH2D("HYBTPC","HYBTPC;TPC ONLY;HYBRID",100,0,3000,100,0,3000); tList->Add( tH2D );
689   tH1D = new TH1D("HYBTPCRat","HYBTPCRat;TPC/HYB",120,0.2,2.2); tList->Add( tH1D );
690   tH2D = new TH2D("SPDVZE","SPDVZE;SPD Tracklets;Total Multiplicity in VZERO",100,0,3500,100,0,25000); tList->Add( tH2D );
691   tH1D = new TH1D("SPDVZERat","SPDVZERat;TotalMultiplicityVZERO/SPDTracklets",120,2,+12); tList->Add( tH1D );
692   if(fReadMC) {
693     tH1D = new TH1D("MCEP","MCEP;MCEP",100,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
694   }
695   fList->Add(tList);
696 }
697 //=======================================================================
698 void AliAnalysisTaskFlowStrange::FillEventSpy(TString name) {
699   ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject("VTXZ"))->Fill( fPriVtxZ, fSPDVtxZ );
700   ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject("CCCC"))->Fill( fV0M, fTRK );
701   ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject("HYBTPC"))->Fill( fRefMultTPC, fRefMultHyb );
702   if(fRefMultHyb>0)
703     ((TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject("HYBTPCRat"))->Fill( double(fRefMultTPC)/double(fRefMultHyb) );
704   ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject("SPDVZE"))->Fill( fSPDtracklets, fVZETotM );
705   if(fSPDtracklets>0)
706     ((TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject("SPDVZERat"))->Fill( fVZETotM/fSPDtracklets );
707   if(fReadMC) {
708     ((TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject("MCEP"))->Fill( fMCEP );
709   }
710 }
711 //=======================================================================
712 void AliAnalysisTaskFlowStrange::AddMakeQSpy() {
713   TH1D *tH1D;
714   TH2D *tH2D;
715   TProfile *tPF1;
716   TList *tList=new TList();
717   tList->SetName("MakeQSpy");
718   tList->SetOwner();
719   fList->Add(tList);
720   tH1D = new TH1D("RFPTPC","TPC Refrence Multiplicity;multiplicity",3000,0,3000);     tList->Add( tH1D );
721   tH1D = new TH1D("RFPVZE","VZERO Reference Multiplicity;multiplicity",3000,0,30000); tList->Add( tH1D );
722   tH1D = new TH1D("QmTPC","TPC Normalized Q vector;|Q|/#sqrt{M}",360,0,7);   tList->Add( tH1D );
723   tH1D = new TH1D("QmVZEA","VZEROA Normalized Q vector;|Q|/#sqrt{W}",360,0,7); tList->Add( tH1D );
724   tH1D = new TH1D("QmVZEC","VZEROC Normalized Q vector;|Q|/#sqrt{W}",360,0,7); tList->Add( tH1D );
725   tH2D = new TH2D("TPCAllPhiEta","TPCall;Phi;Eta",180,0,TMath::TwoPi(),80,-0.9,+0.9); tList->Add( tH2D );
726   tH2D = new TH2D("VZEAllPhiEta","VZEall;Phi;Eta",20,0,TMath::TwoPi(),40,-4.0,+6.0);  tList->Add( tH2D );
727   tH1D = new TH1D("TPCPSI","TPCPSI;PSI",72,0,TMath::Pi()); tList->Add( tH1D );
728   tH1D = new TH1D("TPCPSIA","TPCPSIA;PSIA",72,0,TMath::Pi()); tList->Add( tH1D );
729   tH1D = new TH1D("TPCPSIC","TPCPSIC;PSIC",72,0,TMath::Pi()); tList->Add( tH1D );
730   tH1D = new TH1D("VZEPSI","VZEPSI;PSI",72,0,TMath::Pi()); tList->Add( tH1D );
731   tH1D = new TH1D("VZEPSIA","VZEPSIA;PSIA",72,0,TMath::Pi()); tList->Add( tH1D );
732   tH1D = new TH1D("VZEPSIC","VZEPSIC;PSIC",72,0,TMath::Pi()); tList->Add( tH1D );
733   tH2D = new TH2D("PSI_TPCAVZEC","PSI_TPCAVZEC",72,0,TMath::Pi(),72,0,TMath::Pi()); tList->Add( tH2D );
734   tH2D = new TH2D("PSI_TPCCVZEA","PSI_TPCAVZEC",72,0,TMath::Pi(),72,0,TMath::Pi()); tList->Add( tH2D );
735   tH2D = new TH2D("PSI_TPCVZE","PSI_TPCVZE",72,0,TMath::Pi(),72,0,TMath::Pi()); tList->Add( tH2D );
736   tPF1 = new TProfile("TPCQm","TPCQm",6,0.5,6.5); tList->Add( tPF1 );
737   tPF1->GetXaxis()->SetBinLabel(1,"Qcy"); tPF1->GetXaxis()->SetBinLabel(2,"Qcx");
738   tPF1->GetXaxis()->SetBinLabel(3,"Qay"); tPF1->GetXaxis()->SetBinLabel(4,"Qax");
739   tPF1->GetXaxis()->SetBinLabel(5,"Qy");  tPF1->GetXaxis()->SetBinLabel(6,"Qx");
740   tPF1 = new TProfile("VZEQm","VZEQm",6,0.5,6.5); tList->Add( tPF1 );
741   tPF1->GetXaxis()->SetBinLabel(1,"Qcy"); tPF1->GetXaxis()->SetBinLabel(2,"Qcx");
742   tPF1->GetXaxis()->SetBinLabel(3,"Qay"); tPF1->GetXaxis()->SetBinLabel(4,"Qax");
743   tPF1->GetXaxis()->SetBinLabel(5,"Qy");  tPF1->GetXaxis()->SetBinLabel(6,"Qx");
744   tPF1 = new TProfile("QmVZEAQmVZEC","QmVZEAQmVZEC",1,0.5,1.5,"s"); tList->Add( tPF1 );
745   tPF1 = new TProfile("QmVZEASQUARED","QmVZEASQUARED",1,0.5,1.5,"s"); tList->Add( tPF1 );
746   tPF1 = new TProfile("QmVZECSQUARED","QmVZECSQUARED",1,0.5,1.5,"s"); tList->Add( tPF1 );
747   tPF1 = new TProfile("QmTPCQmVZEA","QmTPCQmVZEA",1,0.5,1.5,"s"); tList->Add( tPF1 );
748   tPF1 = new TProfile("QmTPCQmVZEC","QmTPCQmVZEC",1,0.5,1.5,"s"); tList->Add( tPF1 );
749   tH1D = new TH1D("ChiSquaredVZEA","ChiSquaredVZEC",1,0.5,1.5); tList->Add( tH1D );
750   tH1D = new TH1D("ChiSquaredVZEC","ChiSquaredVZEC",1,0.5,1.5); tList->Add( tH1D );
751   if(fReadMC) {
752     tH1D = new TH1D("PSIMCDIFFTPC","PSIMCDIFFTPC;MC-TPC",72,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
753     tH1D = new TH1D("PSIMCDIFFTPCA","PSIMCDIFFTPCA;MC-TPCA",72,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
754     tH1D = new TH1D("PSIMCDIFFTPCC","PSIMCDIFFTPCC;MC-TPCC",72,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
755     tH1D = new TH1D("PSIMCDIFFVZE","PSIMCDIFFVZE;MC-VZE",72,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
756     tH1D = new TH1D("PSIMCDIFFVZEA","PSIMCDIFFVZEA;MC-VZEA",72,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
757     tH1D = new TH1D("PSIMCDIFFVZEC","PSIMCDIFFVZEC;MC-VZEC",72,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
758   }
759   tList=new TList(); tList->SetName("TPCRFPall"); tList->SetOwner(); AddTPCRFPSpy(tList); fList->Add(tList);
760   tList=new TList(); tList->SetName("TPCRFPsel"); tList->SetOwner(); AddTPCRFPSpy(tList); fList->Add(tList);
761 }
762 //=======================================================================
763 void AliAnalysisTaskFlowStrange::AddQACandidates() {
764   // function to add histogramming for candidates
765   if(fSkipSelection) return;
766   TList *tList;
767   TH1D *tH1D;
768
769   //charge particles (benchmark)
770   if(fSpecie>=90) {
771     tList=new TList(); tList->SetName("TrkAll"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
772     tList=new TList(); tList->SetName("TrkSel"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
773     if(!fSkipVn) {
774       tList=new TList(); tList->SetName("TrkAllVn"); tList->SetOwner(); AddTrackVn(tList); fList->Add(tList);
775       tList=new TList(); tList->SetName("TrkSelVn"); tList->SetOwner(); AddTrackVn(tList); fList->Add(tList);
776     }
777     //match
778     if(fReadMC) {
779       tList=new TList(); tList->SetName("STATMC"); tList->SetOwner(); fList->Add(tList);
780       tH1D = new TH1D("Events", "Events",3,0.5,3.5); tList->Add(tH1D);
781       tH1D->GetXaxis()->SetBinLabel(1,"Selected events");
782       tH1D->GetXaxis()->SetBinLabel(2,"Stack found");
783       tH1D->GetXaxis()->SetBinLabel(3,"Track in stack");
784       tList=new TList(); tList->SetName("Mth"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
785       tList=new TList(); tList->SetName("MthPos"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
786       tList=new TList(); tList->SetName("MthNeg"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
787       if(!fSkipVn) {
788         tList=new TList(); tList->SetName("MthVn"); tList->SetOwner(); AddTrackVn(tList); fList->Add(tList);
789         tList=new TList(); tList->SetName("MthPosVn"); tList->SetOwner(); AddTrackVn(tList); fList->Add(tList);
790         tList=new TList(); tList->SetName("MthNegVn"); tList->SetOwner(); AddTrackVn(tList); fList->Add(tList);
791       }
792     }
793   }
794   //stack
795   if(fReadMC) {
796     tList=new TList(); tList->SetName("MCTPionGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
797     tList=new TList(); tList->SetName("MCTKaonGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
798     tList=new TList(); tList->SetName("MCTK0sGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
799     tList=new TList(); tList->SetName("MCTProtonGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
800     tList=new TList(); tList->SetName("MCTLdaGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
801     tList=new TList(); tList->SetName("MCTPhiGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
802     tList=new TList(); tList->SetName("MCTXiGenAcc");  tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
803     tList=new TList(); tList->SetName("MCTOmegaGenAcc");  tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
804     tList=new TList(); tList->SetName("MCTPion"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
805     tList=new TList(); tList->SetName("MCTKaon"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
806     tList=new TList(); tList->SetName("MCTK0s"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
807     tList=new TList(); tList->SetName("MCTLda"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
808     tList=new TList(); tList->SetName("MCTProton"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
809   }
810   MyUserCreateOutputObjects();
811 }
812 //=======================================================================
813 void AliAnalysisTaskFlowStrange::Exec(Option_t* option) {
814   // bypassing ::exec (needed because of AMPT)
815   ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(0);
816   if(fAvoidExec) {
817     AliAnalysisTaskFlowStrange::UserExec(option);
818   } else {
819     AliAnalysisTaskSE::Exec(option);
820   }
821 }
822 //=======================================================================
823 void AliAnalysisTaskFlowStrange::UserExec(Option_t *option) {
824   // bridge
825   ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(1);
826   AliAnalysisTaskFlowStrange::MyUserExec(option);
827 }
828 //=======================================================================
829 void AliAnalysisTaskFlowStrange::MyNotifyRun() {
830   if(fVZEsave) AddVZEROResponse();
831 }
832 //=======================================================================
833 Bool_t AliAnalysisTaskFlowStrange::CalibrateEvent() {
834   if(fVZEsave) SaveVZEROResponse();
835   Bool_t okay=kTRUE;
836   if(fVZEload) {
837     LoadVZEROResponse();
838     if(!fVZEResponse) okay = kFALSE;
839   }
840   return okay;
841 }
842 //=======================================================================
843 Bool_t AliAnalysisTaskFlowStrange::AcceptAAEvent(AliESDEvent*) {
844   // ESD reading discontinued: TO BE UPDATED
845   /*
846   Double_t acceptEvent=kTRUE;
847   Double_t tTPCVtxZ = tESD->GetPrimaryVertexTPC()->GetZ();
848   if(tESD->GetPrimaryVertexTPC()->GetNContributors()<=0) return kFALSE;
849   Double_t tSPDVtxZ = tESD->GetPrimaryVertexSPD()->GetZ();
850   if(tESD->GetPrimaryVertexSPD()->GetNContributors()<=0) return kFALSE;
851   // EventCuts
852   AliCentrality *cent = tESD->GetCentrality();
853   Double_t cc1, cc2;
854   cc1 = cent->GetCentralityPercentile("V0M");
855   cc2 = cent->GetCentralityPercentile("TRK");
856   TString mycent = fCentMethod;
857   if(fCentMethod.Contains("V0MTRK")) {
858     acceptEvent = TMath::Abs(cc1-cc2)>5.0?kFALSE:acceptEvent; // a la Alex
859     mycent = "V0M";
860   }
861   fThisCent = cent->GetCentralityPercentile( mycent );
862   acceptEvent = (fThisCent<fCentPerMin||fThisCent>fCentPerMax)?kFALSE:acceptEvent;
863   acceptEvent = TMath::Abs(tTPCVtxZ-tSPDVtxZ)>0.5?kFALSE:acceptEvent;
864   acceptEvent = TMath::Abs(tTPCVtxZ)>fVertexZcut?kFALSE:acceptEvent;
865   ((TH2D*)((TList*)fList->FindObject("EventsReached"))->FindObject("VTXZ"))->Fill( tTPCVtxZ, tSPDVtxZ );
866   ((TH2D*)((TList*)fList->FindObject("EventsReached"))->FindObject("CCCC"))->Fill( cc1, cc2 );
867   // EndOfCuts
868   if(acceptEvent) {
869     ((TH2D*)((TList*)fList->FindObject("EventsSelected"))->FindObject("VTXZ"))->Fill( tTPCVtxZ, tSPDVtxZ );
870     ((TH2D*)((TList*)fList->FindObject("EventsSelected"))->FindObject("CCCC"))->Fill( cc1, cc2 );
871   }
872   return acceptEvent;
873   */
874   return kFALSE;
875 }
876 //=======================================================================
877 Bool_t AliAnalysisTaskFlowStrange::MinimumRequirementsAA(AliAODEvent *tAOD) {
878   fRunNumber = tAOD->GetRunNumber();
879   AliCentrality *cent = ((AliVAODHeader*)tAOD->GetHeader())->GetCentralityP();
880   fV0M = cent->GetCentralityPercentile("V0M");
881   fTRK = cent->GetCentralityPercentile("TRK");
882   TString mycent = fCentMethod;
883   if(fCentMethod.Contains("V0MTRK")) {
884     mycent = "V0M";
885   }
886   fThisCent = cent->GetCentralityPercentile( mycent );
887   fPriVtxZ = tAOD->GetPrimaryVertex()->GetZ();
888   fSPDVtxZ = tAOD->GetPrimaryVertexSPD()->GetZ();
889   fSPDtracklets = tAOD->GetTracklets()->GetNumberOfTracklets();
890   fVZETotM = tAOD->GetVZEROData()->GetMTotV0A() + tAOD->GetVZEROData()->GetMTotV0C();
891   int hyb_fb = 272; // for 2010h::AOD086
892   if(fRunNumber>=166529&&fRunNumber<=170593) {
893     hyb_fb = 768; // for 2011h::AOD145
894   }
895   fRefMultTPC = RefMult(tAOD,128);
896   fRefMultHyb = RefMult(tAOD,hyb_fb);
897   if(fReadMC) {
898     fMCEP = -999;
899     AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(tAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
900     if(mcHeader) {
901       fMCEP = mcHeader->GetReactionPlaneAngle();
902       if(fAddPiToMCReactionPlane) fMCEP += (gRandom->Rndm()>0.5)*TMath::Pi();
903     }
904   }
905   // centrality selection health
906   // cut in Vtx 10 & NContributors
907   if(!fSkipCentralitySelection) if(fThisCent<0||fThisCent>100) return kFALSE;
908   // vtx z position compatibility within 5 mm
909   if(TMath::Abs(fPriVtxZ-fSPDVtxZ)>0.5) return kFALSE;
910   if(fExtraEventRejection) {
911     // specific cuts for 2010h (AOD086)
912     if(fRunNumber>=136851&&fRunNumber<=139517) {
913       if(fRefMultTPC>1.118*fRefMultHyb+100) return kFALSE;
914       if(fRefMultTPC<1.118*fRefMultHyb-100) return kFALSE;
915     }
916     // specific cuts for 2011h (AOD145)
917     if(fRunNumber>=166529&&fRunNumber<=170593) {
918       if(fRefMultTPC>1.205*fRefMultHyb+100) return kFALSE;
919       if(fRefMultTPC<1.205*fRefMultHyb-100) return kFALSE;
920     }
921   }
922   return kTRUE;
923 }
924 //=======================================================================
925 Bool_t AliAnalysisTaskFlowStrange::AcceptAAEvent(AliAODEvent *tAOD) {
926   Bool_t minimum = MinimumRequirementsAA(tAOD);
927   FillEventSpy("EventsRaw");
928   if(!minimum) return kFALSE;
929
930   Double_t acceptEvent=kTRUE;
931   TString mycent = fCentMethod;
932   if(fCentMethod.Contains("V0MTRK")) {
933     acceptEvent = TMath::Abs(fV0M-fTRK)>5.0?kFALSE:acceptEvent;
934     mycent = "V0M";
935   }
936   if(!fSkipCentralitySelection) acceptEvent = (fThisCent<fCentPerMin||fThisCent>fCentPerMax)?kFALSE:acceptEvent;
937   acceptEvent = TMath::Abs(fPriVtxZ)>fVertexZcut?kFALSE:acceptEvent;
938   // HISTOGRAMMING
939   FillEventSpy("EventsReached");
940   if(acceptEvent) FillEventSpy("EventsSelected");
941   return acceptEvent;
942 }
943 //=======================================================================
944 Bool_t AliAnalysisTaskFlowStrange::AcceptPPEvent(AliAODEvent*) {
945   // PP reading discontinued: TO BE UPDATED
946   /*
947   Double_t acceptEvent=kTRUE;
948   Double_t tVtxZ = tAOD->GetPrimaryVertex()->GetZ();
949   if(tAOD->GetPrimaryVertex()->GetNContributors()<=0) return kFALSE;
950   Double_t tSPDVtxZ = tAOD->GetPrimaryVertexSPD()->GetZ();
951   if(tAOD->GetPrimaryVertexSPD()->GetNContributors()<=0) return kFALSE;
952   // EventCuts
953   AliCentrality *cent = tAOD->GetHeader()->GetCentralityP();
954   Double_t cc1, cc2;
955   cc1 = cent->GetCentralityPercentile("V0M");
956   cc2 = cent->GetCentralityPercentile("TRK");
957   fThisCent = GetReferenceMultiplicity();
958   //for pp i use fCentPerXXX to select on multiplicity
959   acceptEvent = (fThisCent<fCentPerMin||fThisCent>fCentPerMax)?kFALSE:acceptEvent;
960   acceptEvent = TMath::Abs(tVtxZ-tSPDVtxZ)>0.5?kFALSE:acceptEvent;
961   acceptEvent = TMath::Abs(tVtxZ)>fVertexZcut?kFALSE:acceptEvent;
962   ((TH2D*)((TList*)fList->FindObject("EventsReached"))->FindObject("VTXZ"))->Fill( tVtxZ, tSPDVtxZ );
963   ((TH2D*)((TList*)fList->FindObject("EventsReached"))->FindObject("CCCC"))->Fill( cc1, cc2 );
964   // EndOfCuts
965   if(acceptEvent) {
966     ((TH2D*)((TList*)fList->FindObject("EventsSelected"))->FindObject("VTXZ"))->Fill( tVtxZ, tSPDVtxZ );
967     ((TH2D*)((TList*)fList->FindObject("EventsSelected"))->FindObject("CCCC"))->Fill( cc1, cc2 );
968   }
969   return acceptEvent;
970   */
971   return kFALSE;
972 }
973 //=======================================================================
974 Int_t AliAnalysisTaskFlowStrange::GetReferenceMultiplicity() { //toberefined
975   AliAODEvent *tAOD = (AliAODEvent *) InputEvent();
976   if(!tAOD) return -1;
977   AliAODTrack *track;
978   Int_t rawN = tAOD->GetNumberOfTracks();
979   Int_t ref=0;
980   for(Int_t id=0; id!=rawN; ++id) {
981     track = tAOD->GetTrack(id);
982     if(!track->TestFilterBit(fRFPFilterBit)) continue;
983     ++ref;
984   }
985   return ref;
986 }
987 //=======================================================================
988 Bool_t AliAnalysisTaskFlowStrange::AcceptPAEvent(AliAODEvent*) {
989   // PA reading discontinued: TO BE UPDATED
990   /*
991   //if(aod->GetHeader()->GetEventNumberESDFile() == 0) return; //rejecting first chunk NOT NEEDED ANYMORE
992   Int_t bc2 = ((AliVAODHeader*)tAOD->GetHeader())->GetIRInt2ClosestInteractionMap();
993   if(bc2!=0) return kFALSE;
994   Int_t bc1 = ((AliVAODHeader*)tAOD->GetHeader())->GetIRInt1ClosestInteractionMap();
995   if(bc1!=0) return kFALSE;
996   Short_t isPileup = tAOD->IsPileupFromSPD(5);
997   if(isPileup!=0) return kFALSE;
998   if(tAOD->GetHeader()->GetRefMultiplicityComb08()<0) return kFALSE;
999
1000   const AliAODVertex* spdVtx = tAOD->GetPrimaryVertexSPD();
1001   if(!spdVtx) return kFALSE;
1002   if(spdVtx->GetNContributors()<=0) return kFALSE;
1003
1004   const AliAODVertex* tpcVtx=NULL;
1005   Int_t nVertices = tAOD->GetNumberOfVertices();
1006   for(Int_t iVertices = 0; iVertices < nVertices; iVertices++){
1007     const AliAODVertex* vertex = tAOD->GetVertex(iVertices);
1008     if (vertex->GetType() != AliAODVertex::kMainTPC) continue;
1009     tpcVtx = vertex;
1010   }
1011   if(!tpcVtx) return kFALSE;
1012   if(tpcVtx->GetNContributors()<=0) return kFALSE;
1013   Double_t tTPCVtxZ = tpcVtx->GetZ();
1014   Double_t tSPDVtxZ = spdVtx->GetZ();
1015   if (TMath::Abs(tSPDVtxZ - tTPCVtxZ)>2.0) return kFALSE;
1016   if(plpMV(tAOD)) return kFALSE;
1017
1018   Double_t acceptEvent=kTRUE;
1019   // EventCuts
1020   AliCentrality *cent = tAOD->GetHeader()->GetCentralityP();
1021   Double_t cc1, cc2;
1022   cc1 = cent->GetCentralityPercentile("V0M");
1023   cc2 = cent->GetCentralityPercentile("TRK");
1024   if(fCentMethod.Contains("V0MTRK")) fCentMethod = "V0M";
1025   fThisCent = cent->GetCentralityPercentile( fCentMethod );
1026   acceptEvent = (fThisCent<fCentPerMin||fThisCent>fCentPerMax)?kFALSE:acceptEvent;
1027   acceptEvent = TMath::Abs(tTPCVtxZ)>fVertexZcut?kFALSE:acceptEvent;
1028   // EndOfCuts
1029   ((TH2D*)((TList*)fList->FindObject("EventsReached"))->FindObject("VTXZ"))->Fill( tTPCVtxZ, tSPDVtxZ );
1030   ((TH2D*)((TList*)fList->FindObject("EventsReached"))->FindObject("CCCC"))->Fill( cc1, cc2 );
1031   if(acceptEvent) {
1032     ((TH2D*)((TList*)fList->FindObject("EventsSelected"))->FindObject("VTXZ"))->Fill( tTPCVtxZ, tSPDVtxZ );
1033     ((TH2D*)((TList*)fList->FindObject("EventsSelected"))->FindObject("CCCC"))->Fill( cc1, cc2 );    
1034   }
1035   return acceptEvent;
1036   */
1037   return kFALSE;
1038 }
1039 //=======================================================================
1040 void AliAnalysisTaskFlowStrange::MyUserExec(Option_t *) {
1041   // MAIN ROUTINE
1042   TStopwatch tTime;
1043   tTime.Start();
1044   if(fDebug) {
1045     printf("****************\n");
1046     printf("****************\n");
1047     printf("**::MyUserExec()\n");
1048   }
1049   if(fUseFP) fCandidates->SetLast(-1);
1050   AliESDEvent *tESD=dynamic_cast<AliESDEvent*>(InputEvent());
1051   AliAODEvent *tAOD=dynamic_cast<AliAODEvent*>(InputEvent());
1052   Int_t prevRun = fRunNumber;
1053   //=>check event
1054   Bool_t acceptEvent=kFALSE;
1055   if(fReadESD) {
1056     if(!tESD) {ResetContainers(); Publish(); return;}
1057     acceptEvent = fRunOnpp?kFALSE:fRunOnpA?kFALSE:AcceptAAEvent(tESD);
1058   } else {
1059     if(!tAOD) {ResetContainers(); Publish(); return;}
1060     acceptEvent = fRunOnpp?AcceptPPEvent(tAOD):fRunOnpA?AcceptPAEvent(tAOD):AcceptAAEvent(tAOD);
1061   }
1062   if(prevRun!=fRunNumber) {
1063     MyNotifyRun();
1064   }
1065   ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(2);
1066   //=>does the event clear?
1067   if(!acceptEvent) {ResetContainers(); Publish(); return;}
1068   // healthy event incomming
1069   if( !CalibrateEvent() ) { // saves/retrieves/qas VZEROCAL
1070     ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(5);
1071     ResetContainers(); Publish(); return; // issue retrieving callibration
1072   }
1073   // loads Q vectors
1074   MakeQVectors();
1075   if(fPsi2<-0.1) {
1076     ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(4);
1077     ResetContainers(); Publish(); return;
1078   }
1079   //}
1080   ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(3);
1081   //=>great, lets do our stuff!
1082   FillEventSpy("EventsAnalyzed");
1083   FillVZEQA();
1084   //=>load candidates
1085   if(!fSkipSelection) {
1086     if(fReadESD) {
1087       ReadFromESD(tESD);
1088     } else {
1089       if(fSpecie<10) ReadFromAODv0(tAOD);
1090       else ChargeParticles(tAOD);
1091     }
1092     if(fUseFP) AddCandidates();
1093     //=>flow
1094     //=>done
1095   }
1096   tTime.Stop();
1097   ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("RealTime"))->Fill( TMath::Log( tTime.RealTime() ) );
1098   Publish();
1099 }
1100 //=======================================================================
1101 void AliAnalysisTaskFlowStrange::Publish() {
1102   PostData(1,fList);
1103   if(fUseFP) {
1104     PostData(2,fTPCevent);
1105     PostData(3,fVZEevent);
1106   }
1107 }
1108 //=======================================================================
1109 void AliAnalysisTaskFlowStrange::ReadFromESD(AliESDEvent *tESD) {
1110   AliStack *stack=NULL;
1111   if(fReadMC) {
1112     AliMCEvent *mcevent=NULL;
1113     mcevent = MCEvent();
1114     if(mcevent) stack = mcevent->Stack();
1115   }
1116
1117   Int_t num = tESD->GetNumberOfTracks();
1118   AliESDtrack *myTrack;
1119   Int_t plist[3000], nlist[3000], np=0, nn=0;
1120   Double_t pd0[3000], nd0[3000];
1121   for (Int_t i=0; i!=num; ++i) {
1122     myTrack = (AliESDtrack*) tESD->GetTrack(i);
1123     if(!myTrack) continue;
1124     LoadTrack(myTrack);
1125     FillTrackSpy("ESD_TrkAll");
1126     if(!AcceptDaughter()) continue;
1127     FillTrackSpy("ESD_TrkSel");
1128     ((TH2D*)((TList*)fList->FindObject("ESD_TrkSel"))->FindObject("PtIPXY" ))->Fill( myTrack->Pt(), fDaughterImpactParameterXY );
1129     if( myTrack->Charge()>0 ) {
1130       pd0[np] = fDaughterImpactParameterXY;
1131       plist[np++] = i;
1132     } else {
1133       nd0[nn] = fDaughterImpactParameterXY;
1134       nlist[nn++] = i;
1135     }
1136   }
1137   ((TH1D*)((TList*)fList->FindObject("ESD_TrkSel"))->FindObject("NPAIR" ))->Fill( np,nn );
1138   const AliESDVertex *vtx = tESD->GetPrimaryVertex();
1139   AliESDtrack *pT, *nT;
1140   for(int p=0; p!=np; ++p) {
1141     pT = (AliESDtrack*) tESD->GetTrack( plist[p] );
1142     for(int n=0; n!=nn; ++n) {
1143       nT = (AliESDtrack*) tESD->GetTrack( nlist[n] );
1144       fDecayProductIPXY = pd0[p]*nd0[n];
1145       AliExternalTrackParam pETP(*pT), nETP(*nT);
1146       Double_t xa, xb;
1147       pETP.GetDCA(&nETP,tESD->GetMagneticField(),xa,xb);
1148       fDecayDCAdaughters = pETP.PropagateToDCA(&nETP,tESD->GetMagneticField());
1149       AliESDv0 vertex( nETP,nlist[n], pETP,plist[p] );
1150       fDecayCosinePointingAngleXY = CosThetaPointXY( &vertex, vtx );
1151       fDecayRadXY = DecayLengthXY( &vertex, vtx );
1152       fDecayPt = vertex.Pt();
1153       fDecayPhi = vertex.Phi();
1154       fDecayEta = vertex.Eta();
1155       Double_t pmx, pmy, pmz, nmx, nmy, nmz;
1156       vertex.GetNPxPyPz(nmx,nmy,nmz);
1157       vertex.GetPPxPyPz(pmx,pmy,pmz);
1158       TVector3 mom1(pmx,pmy,pmz), mom2(nmx,nmy,nmz), mom(vertex.Px(),vertex.Py(),vertex.Pz());
1159       Double_t qlpos = mom1.Dot(mom)/mom.Mag();
1160       Double_t qlneg = mom2.Dot(mom)/mom.Mag();
1161       fDecayQt = mom1.Perp(mom);
1162       fDecayAlpha = (qlpos-qlneg)/(qlpos+qlneg);
1163       Double_t mpi = 0.13957018;
1164       if(fSpecie==0) {
1165         Double_t eppi = TMath::Sqrt( mpi*mpi + pmx*pmx + pmy*pmy + pmz*pmz );
1166         Double_t enpi = TMath::Sqrt( mpi*mpi + nmx*nmx + nmy*nmy + nmz*nmz );
1167         fDecayMass = TMath::Sqrt( mpi*mpi + mpi*mpi + 2*(eppi*enpi - pmx*nmx - pmy*nmy - pmz*nmz ) );
1168         fDecayRapidity = vertex.RapK0Short();
1169       } else {
1170         Double_t mpr = 0.938272013;
1171         Double_t epi, epr;
1172         if(fDecayAlpha>0) {
1173           epr = TMath::Sqrt( mpr*mpr + pmx*pmx + pmy*pmy + pmz*pmz );
1174           epi = TMath::Sqrt( mpi*mpi + nmx*nmx + nmy*nmy + nmz*nmz );
1175         } else {
1176           epi = TMath::Sqrt( mpi*mpi + pmx*pmx + pmy*pmy + pmz*pmz );
1177           epr = TMath::Sqrt( mpr*mpr + nmx*nmx + nmy*nmy + nmz*nmz );
1178         }
1179         fDecayMass = TMath::Sqrt( mpi*mpi + mpr*mpr + 2*(epi*epr - pmx*nmx - pmy*nmy - pmz*nmz ) );
1180         fDecayRapidity = vertex.RapLambda();
1181       }
1182       Double_t energy = TMath::Sqrt( fDecayMass*fDecayMass + vertex.Px()*vertex.Px() + vertex.Py()*vertex.Py() + vertex.Pz()*vertex.Pz() );
1183       Double_t gamma = energy/fDecayMass;
1184       fDecayDecayLength = DecayLength( &vertex, vtx )/gamma;
1185       fDecayDecayLengthLab = DecayLength( &vertex, vtx );
1186       Double_t dPHI = fDecayPhi;
1187       Double_t dDPHI = dPHI - fPsi2;
1188       if( dDPHI < 0 ) dDPHI += TMath::TwoPi();
1189       if( dDPHI > TMath::Pi() ) dDPHI = TMath::TwoPi()-dDPHI;
1190       if(fQAlevel>1) {
1191         if( (dDPHI>TMath::PiOver4()) && (dDPHI<3*TMath::PiOver4()) ) FillCandidateSpy("V0SAllOP");
1192         else FillCandidateSpy("V0SAllIP");
1193       }
1194       FillCandidateSpy("V0SAll");
1195       ((TH2D*)((TList*)fList->FindObject("V0SAll"))->FindObject("D0PD0N"))->Fill( pd0[p],nd0[n] );
1196       ((TH2D*)((TList*)fList->FindObject("V0SAll"))->FindObject("XPOSXNEG"))->Fill( xa, xb );
1197       if(!AcceptCandidate()) continue;
1198       if(fDecayMass<fMinMass) continue;
1199       if(fDecayMass>fMaxMass) continue;
1200       // PID missing
1201       if(fQAlevel>1) {
1202         if( (dDPHI>TMath::PiOver4()) && (dDPHI<3*TMath::PiOver4()) ) FillCandidateSpy("V0SSelOP");
1203         else FillCandidateSpy("V0SSelIP");
1204       }
1205       FillCandidateSpy("V0SSel");
1206       ((TH2D*)((TList*)fList->FindObject("V0SSel"))->FindObject("D0PD0N"))->Fill( pd0[p],nd0[n] );
1207       ((TH2D*)((TList*)fList->FindObject("V0SSel"))->FindObject("XPOSXNEG"))->Fill( xa, xb );
1208
1209       fDecayIDneg = nT->GetID();
1210       fDecayIDpos = pT->GetID();
1211       if(fUseFP) MakeTrack();
1212       LoadTrack(pT); FillTrackSpy("SelDau");
1213       LoadTrack(nT); FillTrackSpy("SelDau");
1214
1215       //===== BEGIN OF MCMATCH
1216       if(stack) {
1217         bool matched = false;
1218         Int_t labelpos = pT->GetLabel();
1219         Int_t labelneg = nT->GetLabel();
1220         Double_t rOri=-1;
1221         if( labelpos>0 && labelneg>0 ) {
1222           TParticle *mcpos = stack->Particle( labelpos );
1223           TParticle *mcneg = stack->Particle( labelneg );
1224           Int_t pdgRecPos = mcpos->GetPdgCode();
1225           Int_t pdgRecNeg = mcneg->GetPdgCode();
1226           if( pdgRecPos==211&&pdgRecNeg==-211 ) if(mcpos->GetMother(0)>0) {
1227             if( mcpos->GetMother(0)==mcneg->GetMother(0) ) {
1228               TParticle *mcmot = stack->Particle( mcpos->GetMother(0) );
1229               rOri = TMath::Sqrt( mcmot->Vx()*mcmot->Vx() + mcmot->Vy()*mcmot->Vy() );
1230               if( TMath::Abs(mcmot->GetPdgCode())==310) {
1231                 if(mcmot->GetNDaughters()==2) matched=true;
1232               }
1233             }
1234           }
1235         }
1236         if(matched) {
1237           FillCandidateSpy("Mth");
1238           ((TH2D*)((TList*)fList->FindObject("Mth"))->FindObject("D0PD0N"))->Fill( pd0[p],nd0[n] );
1239           ((TH2D*)((TList*)fList->FindObject("Mth"))->FindObject("XPOSXNEG"))->Fill( xa, xb );
1240           ((TH1D*)((TList*)fList->FindObject("Mth"))->FindObject("MCOrigin"))->Fill( rOri );
1241           LoadTrack(pT); FillTrackSpy("MthDau");
1242           LoadTrack(nT); FillTrackSpy("MthDau");
1243         }
1244       }
1245       //===== END OF MCMATCH
1246     }
1247   }
1248 }
1249 //=======================================================================
1250 void AliAnalysisTaskFlowStrange::ReadStack(TClonesArray* mcArray) {
1251   if(!mcArray) return;
1252   AliAODMCParticle *myMCTrack;//, *iMCDau, *jMCDau;
1253   for(int i=0; i!=mcArray->GetEntriesFast(); ++i) {
1254     myMCTrack = dynamic_cast<AliAODMCParticle*>(mcArray->At( i ));
1255     if(!myMCTrack) continue;
1256     /*
1257     int tPDG=310;
1258     if(fSpecie>0) tPDG = 3122;
1259     if( TMath::Abs(myMCTrack->PdgCode())==tPDG )
1260       if( myMCTrack->GetNDaughters() == 2 ) {
1261         Int_t iDau = myMCTrack->GetDaughter(0);
1262         Int_t jDau = myMCTrack->GetDaughter(1);
1263         AliAODMCParticle *posDau=NULL;
1264         AliAODMCParticle *negDau=NULL;
1265         if(iDau>0&&jDau>0) {
1266           iMCDau = dynamic_cast<AliAODMCParticle*>(mcArray->At( iDau ));
1267           jMCDau = dynamic_cast<AliAODMCParticle*>(mcArray->At( jDau ));
1268           if(iMCDau) {
1269             if(iMCDau->Charge()>0) posDau=iMCDau;
1270             else negDau=iMCDau;
1271           }
1272           if(jMCDau) {
1273             if(jMCDau->Charge()>0) posDau=jMCDau;
1274             else negDau=jMCDau;
1275           }
1276         } //got two daughters
1277         if(posDau&&negDau) {
1278           Double_t dx = myMCTrack->Xv() - posDau->Xv();
1279           Double_t dy = myMCTrack->Yv() - posDau->Yv();
1280           Double_t dz = myMCTrack->Zv() - posDau->Zv();
1281           fDecayRadXY = TMath::Sqrt( dx*dx + dy*dy );
1282           TVector3 momPos(posDau->Px(),posDau->Py(),posDau->Pz());
1283           TVector3 momNeg(negDau->Px(),negDau->Py(),negDau->Pz());
1284           TVector3 momTot(myMCTrack->Px(),myMCTrack->Py(),myMCTrack->Pz());
1285           Double_t qlpos = momPos.Dot(momTot)/momTot.Mag();
1286           Double_t qlneg = momNeg.Dot(momTot)/momTot.Mag();
1287           fDecayQt = momPos.Perp(momTot);
1288           fDecayAlpha = 1.-2./(1.+qlpos/qlneg);
1289           fDecayMass = myMCTrack->GetCalcMass();
1290           Double_t energy = myMCTrack->E();
1291           Double_t gamma = energy/fDecayMass;
1292           fDecayDecayLength = TMath::Sqrt(dx*dx+dy*dy+dz*dz)/gamma;
1293           fDecayPt = myMCTrack->Pt();
1294           fDecayPhi = myMCTrack->Phi();
1295           fDecayEta = myMCTrack->Eta();
1296           fDecayRapidity = myMCTrack->Y();
1297           fDecayDCAdaughters = 0;
1298           fDecayCosinePointingAngleXY = 1;
1299           fDecayProductIPXY = -1;
1300           if(AcceptCandidate()) FillCandidateSpy("GenTru");
1301         }
1302       } // k0/lda with two daughters
1303     */
1304     //==== BEGIN TRACK CUTS
1305     if(myMCTrack->Eta()<-0.8) continue;
1306     if(myMCTrack->Eta()>+0.8) continue;
1307     if(myMCTrack->Y()<-0.5) continue;
1308     if(myMCTrack->Y()>+0.5) continue;
1309     //==== END TRACK CUTS
1310     switch( TMath::Abs(myMCTrack->PdgCode()) ) {
1311     case (211): //pi
1312       FillMCParticleSpy( "MCTPion", myMCTrack );
1313       if( myMCTrack->IsPrimary() )
1314         FillMCParticleSpy( "MCTPionGenAcc", myMCTrack );
1315       break;
1316     case (321): //kaon
1317       FillMCParticleSpy( "MCTKaon", myMCTrack );
1318       if( myMCTrack->IsPrimary() )
1319         FillMCParticleSpy( "MCTKaonGenAcc", myMCTrack );
1320       break;
1321     case (310): //k0s
1322       FillMCParticleSpy( "MCTK0s", myMCTrack );
1323       if( myMCTrack->IsPrimary() )
1324         FillMCParticleSpy( "MCTK0sGenAcc", myMCTrack );
1325       break;
1326     case (2212): //proton
1327       FillMCParticleSpy( "MCTProton", myMCTrack );
1328       if( myMCTrack->IsPrimary() )
1329         FillMCParticleSpy( "MCTProtonGenAcc", myMCTrack );
1330       break;
1331     case (3122): //lda
1332       FillMCParticleSpy( "MCTLda", myMCTrack );
1333       if( myMCTrack->IsPrimary() )
1334         FillMCParticleSpy( "MCTLdaGenAcc", myMCTrack );
1335       break;
1336     case (333): //phi
1337       if( myMCTrack->IsPrimary() )
1338         FillMCParticleSpy( "MCTPhiGenAcc", myMCTrack );
1339       break;
1340     case (3312): //xi
1341       if( myMCTrack->IsPrimary() )
1342         FillMCParticleSpy( "MCTXiGenAcc", myMCTrack );
1343       break;
1344     case (3334): //omega
1345       if( myMCTrack->IsPrimary() )
1346         FillMCParticleSpy( "MCTOmegaGenAcc", myMCTrack );
1347       break;
1348     }
1349   }
1350 }
1351 //=======================================================================
1352 Double_t AliAnalysisTaskFlowStrange::CosThetaPointXY(AliESDv0 *me, const AliVVertex *vtx) {
1353   TVector3 mom( me->Px(), me->Py(), 0 );
1354   TVector3 fli( me->Xv()-vtx->GetX(), me->Yv()-vtx->GetY(), 0 );
1355   Double_t ctp = mom.Dot(fli) / mom.Mag() / fli.Mag();
1356   return ctp;
1357 }
1358 //=======================================================================
1359 Double_t AliAnalysisTaskFlowStrange::CosThetaPointXY(AliAODv0 *me, const AliVVertex *vtx) {
1360   TVector3 mom( me->Px(), me->Py(), 0 );
1361   TVector3 fli( me->Xv()-vtx->GetX(), me->Yv()-vtx->GetY(), 0 );
1362   Double_t ctp = mom.Dot(fli) / mom.Mag() / fli.Mag();
1363   return ctp;
1364 }
1365 //=======================================================================
1366 Double_t AliAnalysisTaskFlowStrange::DecayLengthXY(AliESDv0 *me, const AliVVertex *vtx) {
1367   Double_t dx = me->Xv()-vtx->GetX();
1368   Double_t dy = me->Yv()-vtx->GetY();
1369   Double_t dxy = TMath::Sqrt( dx*dx + dy*dy );
1370   return dxy;
1371 }
1372 //=======================================================================
1373 Double_t AliAnalysisTaskFlowStrange::DecayLengthXY(AliAODv0 *me, const AliVVertex *vtx) {
1374   Double_t dx = me->Xv()-vtx->GetX();
1375   Double_t dy = me->Yv()-vtx->GetY();
1376   Double_t dxy = TMath::Sqrt( dx*dx + dy*dy );
1377   return dxy;
1378 }
1379 //=======================================================================
1380 Double_t AliAnalysisTaskFlowStrange::DecayLength(AliESDv0 *me, const AliVVertex *vtx) {
1381   Double_t dx = me->Xv()-vtx->GetX();
1382   Double_t dy = me->Yv()-vtx->GetY();
1383   Double_t dz = me->Zv()-vtx->GetZ();
1384   Double_t dxy = TMath::Sqrt( dx*dx + dy*dy + dz*dz );
1385   return dxy;
1386 }
1387 //=======================================================================
1388 Double_t AliAnalysisTaskFlowStrange::DecayLength(AliAODv0 *me, const AliVVertex *vtx) {
1389   Double_t dx = me->Xv()-vtx->GetX();
1390   Double_t dy = me->Yv()-vtx->GetY();
1391   Double_t dz = me->Zv()-vtx->GetZ();
1392   Double_t dxy = TMath::Sqrt( dx*dx + dy*dy + dz*dz );
1393   return dxy;
1394 }
1395 //=======================================================================
1396 void AliAnalysisTaskFlowStrange::ReadFromAODv0(AliAODEvent *tAOD) {
1397   TClonesArray* mcArray=NULL;
1398   if(fReadMC) {
1399     mcArray = dynamic_cast<TClonesArray*>(tAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1400     ReadStack(mcArray);
1401   }
1402
1403   Int_t nV0s = tAOD->GetNumberOfV0s();
1404   AliAODv0 *myV0;
1405   Int_t v0all=0, v0imw=0;
1406   for (Int_t i=0; i!=nV0s; ++i) {
1407     myV0 = (AliAODv0*) tAOD->GetV0(i);
1408     if(!myV0) continue;
1409     if(!fOnline) if(myV0->GetOnFlyStatus() ) continue;
1410     if(fOnline) if(!myV0->GetOnFlyStatus() ) continue;
1411
1412     fDecayPt = myV0->Pt();
1413     fDecayPhi = myV0->Phi();
1414     fDecayEta = myV0->Eta();
1415
1416     AliAODTrack *iT, *jT;
1417     AliAODVertex *vtx = tAOD->GetPrimaryVertex();
1418     Double_t pos[3],cov[6];
1419     vtx->GetXYZ(pos);
1420     vtx->GetCovarianceMatrix(cov);
1421     const AliESDVertex vESD(pos,cov,100.,100);
1422     // TESTING CHARGE
1423     int iPos, iNeg;
1424     iT=(AliAODTrack*) myV0->GetDaughter(0);
1425     if(iT->Charge()>0) {
1426       iPos = 0; iNeg = 1;
1427     } else {
1428       iPos = 1; iNeg = 0;
1429     }
1430     // END OF TEST
1431
1432     iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1433     AliESDtrack ieT( iT );
1434     ieT.SetTPCClusterMap( iT->GetTPCClusterMap() );
1435     ieT.SetTPCSharedMap( iT->GetTPCSharedMap() );
1436     ieT.SetTPCPointsF( iT->GetTPCNclsF() );
1437     ieT.PropagateToDCA(&vESD, tAOD->GetMagneticField(), 100);
1438     LoadTrack(&ieT,iT->Chi2perNDF());
1439     Float_t ip[2];
1440     ieT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1441     fDaughterImpactParameterXY = ip[0];
1442     fDaughterImpactParameterZ = ip[1];
1443     fDecayIPpos = fDaughterImpactParameterXY; //ieT.GetD(pos[0], pos[1], tAOD->GetMagneticField());
1444     FillTrackSpy("AllDau");
1445     if(!AcceptDaughter(fDecayPt<2.0?kTRUE:kFALSE)) continue;
1446
1447     jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1448     AliESDtrack jeT( jT );
1449     jeT.SetTPCClusterMap( jT->GetTPCClusterMap() );
1450     jeT.SetTPCSharedMap( jT->GetTPCSharedMap() );
1451     jeT.SetTPCPointsF( jT->GetTPCNclsF() );
1452     jeT.PropagateToDCA(&vESD, tAOD->GetMagneticField(), 100);
1453     LoadTrack(&jeT,jT->Chi2perNDF());
1454     jeT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1455     fDaughterImpactParameterXY = ip[0];
1456     fDaughterImpactParameterZ = ip[1];
1457     fDecayIPneg = fDaughterImpactParameterXY; //jeT.GetD(pos[0], pos[1], tAOD->GetMagneticField());
1458     FillTrackSpy("AllDau");
1459     if(!AcceptDaughter(fDecayPt<2.0?kTRUE:kFALSE)) continue;
1460
1461     if( fExcludeTPCEdges ) {
1462       if( IsAtTPCEdge(iT->Phi(),iT->Pt(),+1,tAOD->GetMagneticField()) ) continue;
1463       if( IsAtTPCEdge(jT->Phi(),jT->Pt(),-1,tAOD->GetMagneticField()) ) continue;
1464     }
1465     ieT.GetDCA(&jeT,tAOD->GetMagneticField(),fDecayXpos,fDecayXneg);
1466     /*
1467     // cutting out population close to TPC edges :: strange excess saw in 2010
1468     if( fExcludeTPCEdges ) {
1469     Double_t phimod = myV0->Phi();
1470     int sectors[6] = {5,6,9,10,11,12};
1471     for(int ii=0; ii!=6; ++ii)
1472     if( (phimod<(sectors[ii]+1)*TMath::Pi()/9.0) && (phimod>sectors[ii]*TMath::Pi()/9.0) )
1473     return 0;
1474     }
1475     */
1476     if(fSpecie==0)
1477       fDecayRapidity = myV0->RapK0Short();
1478     else
1479       fDecayRapidity = myV0->RapLambda();
1480     fDecayDCAdaughters = myV0->DcaV0Daughters();
1481     fDecayCosinePointingAngleXY = CosThetaPointXY( myV0, vtx );
1482     fDecayRadXY = DecayLengthXY( myV0, vtx );
1483     fDecayProductIPXY = fDecayIPpos*fDecayIPneg;
1484     fDecayQt = myV0->PtArmV0();
1485     fDecayAlpha = myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1486     if(myV0->ChargeProng(iPos)<0) fDecayAlpha = -fDecayAlpha; // protects for a change in convention
1487     fDecayPt = myV0->Pt();
1488     fDecayEta = myV0->Eta();
1489     if( fSpecie==0 ) {
1490       fDecayMass = myV0->MassK0Short();
1491     } else {
1492       if(fDecayAlpha>0) fDecayMass = myV0->MassLambda();
1493       else fDecayMass = myV0->MassAntiLambda();
1494     }
1495     v0all++;
1496     if(fDecayMass<fMinMass) continue;
1497     if(fDecayMass>fMaxMass) continue;
1498     v0imw++;
1499     Double_t energy = TMath::Sqrt( fDecayMass*fDecayMass + myV0->Px()*myV0->Px() + myV0->Py()*myV0->Py() + myV0->Pz()*myV0->Pz() );
1500     Double_t gamma = energy/fDecayMass;
1501     fDecayDecayLength = DecayLength( myV0, vtx )/gamma;
1502     fDecayDecayLengthLab = DecayLength( myV0, vtx );
1503     Double_t dPHI = fDecayPhi;
1504     Double_t dDPHI = dPHI - fPsi2;
1505     if( dDPHI < 0 ) dDPHI += TMath::TwoPi();
1506     if( dDPHI > TMath::Pi() ) dDPHI = TMath::TwoPi()-dDPHI;
1507     if(fQAlevel>1) {
1508       if( (dDPHI>TMath::PiOver4()) && (dDPHI<3*TMath::PiOver4()) ) FillCandidateSpy("V0SAllOP");
1509       else FillCandidateSpy("V0SAllIP");
1510     }
1511     FillCandidateSpy("V0SAll");
1512     if(!fSkipVn)
1513       FillDecayVn("V0SAllVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
1514     
1515     if(!AcceptCandidate()) continue;
1516
1517     if(fDecayPt<fDecayStopPIDAtPt) {
1518       if( fSpecie==0 ) {//PID for kzero::pion+pion
1519         if( !PassesPIDCuts(&ieT,AliPID::kPion) ) continue; //positive track
1520         if( !PassesPIDCuts(&jeT,AliPID::kPion) ) continue; //negative track
1521       } else { //PID for lambda::proton+pion
1522         if(fDecayAlpha>0) {
1523           if( !PassesPIDCuts(&ieT,AliPID::kProton) ) continue; //positive track
1524           if( !PassesPIDCuts(&jeT,AliPID::kPion) ) continue; //negative track
1525         } else {
1526           if( !PassesPIDCuts(&jeT,AliPID::kProton) ) continue; //negative track
1527           if( !PassesPIDCuts(&ieT,AliPID::kPion) ) continue; //positive track
1528         }
1529       }
1530     }
1531     if(fQAlevel>1) {
1532       if( (dDPHI>TMath::PiOver4()) && (dDPHI<3*TMath::PiOver4()) ) FillCandidateSpy("V0SSelOP");
1533       else FillCandidateSpy("V0SSelIP");
1534     }
1535     FillCandidateSpy("V0SSel");
1536     if(!fSkipVn)
1537       FillDecayVn("V0SSelVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
1538     // ============================
1539     // Posting for FlowAnalysis
1540     if(!fPostMatched) {
1541       fDecayIDneg = iT->GetID();
1542       fDecayIDpos = jT->GetID();
1543       if(fUseFP) MakeTrack();
1544     }
1545     // ============================
1546     LoadTrack(&ieT,iT->Chi2perNDF());
1547     ieT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1548     fDaughterImpactParameterXY = ip[0];
1549     fDaughterImpactParameterZ = ip[1];
1550     FillTrackSpy("SelDau");
1551     LoadTrack(&jeT,jT->Chi2perNDF()); 
1552     jeT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1553     fDaughterImpactParameterXY = ip[0];
1554     fDaughterImpactParameterZ = ip[1];
1555     FillTrackSpy("SelDau");
1556     //===== BEGIN OF MCMATCH
1557     if(fReadMC) ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 1 ); // Selected event
1558     if(mcArray) {
1559       ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 2 ); // Stack found
1560       bool matched = false;
1561       bool feeddown = false;
1562       Int_t labelpos = iT->GetLabel();
1563       Int_t labelneg = jT->GetLabel();
1564       AliAODMCParticle *mcpos = (AliAODMCParticle*) mcArray->At( TMath::Abs(labelpos) );
1565       AliAODMCParticle *mcneg = (AliAODMCParticle*) mcArray->At( TMath::Abs(labelneg) );
1566       if( mcpos && mcneg ) {
1567         ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 3 ); // Daughters in stack
1568         Int_t pdgRecPos = mcpos->GetPdgCode();
1569         Int_t pdgRecNeg = mcneg->GetPdgCode();
1570         int pospdg=211, negpdg=211;
1571         int mompdg=310, fdwpdg=333;
1572         if(fSpecie>0) {
1573           mompdg=3122;
1574           fdwpdg=3312;
1575           if(fDecayAlpha>0) {
1576             pospdg=2212; negpdg=211;
1577           } else {
1578             negpdg=2212; pospdg=211;
1579           }
1580         }
1581         if( TMath::Abs(pdgRecPos)==pospdg&&TMath::Abs(pdgRecNeg)==negpdg )
1582           if(mcpos->GetMother()>-1)
1583             if( mcpos->GetMother()==mcneg->GetMother() ) {
1584               AliAODMCParticle *mcmot = (AliAODMCParticle*) mcArray->At( mcpos->GetMother() );
1585               fDecayMatchOrigin = TMath::Sqrt( mcmot->Xv()*mcmot->Xv() + mcmot->Yv()*mcmot->Yv() );
1586               fDecayMatchPt = mcmot->Pt();
1587               fDecayMatchEta = mcmot->Eta();
1588               fDecayMatchPhi = mcmot->Phi();
1589               if( TMath::Abs(mcmot->GetPdgCode())==mompdg) {
1590                 if(mcmot->GetNDaughters()==2) {
1591                   ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 4 ); // Correspond to decay
1592                   matched=true;
1593                   Double_t dx = mcmot->Xv() - mcpos->Xv();
1594                   Double_t dy = mcmot->Yv() - mcpos->Yv();
1595                   fDecayMatchRadXY = TMath::Sqrt( dx*dx + dy*dy );
1596                 }
1597                 if(mcmot->GetMother()>-1) {
1598                   ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 5 ); // Decay has mother
1599                   AliAODMCParticle *mcfdw = (AliAODMCParticle*) mcArray->At( mcmot->GetMother() );
1600                   if( TMath::Abs(mcfdw->GetPdgCode())==fdwpdg)
1601                     feeddown=true;
1602                 } // k0/lda have mother
1603               } // mother matches k0/lda
1604             } // both have same mother
1605       }
1606       if(matched) {
1607         FillCandidateSpy("Mth",true);
1608         if(!fSkipVn)
1609           FillDecayVn("V0SMthVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
1610         if(fPostMatched>0) {
1611           fDecayIDneg = iT->GetID();
1612           fDecayIDpos = jT->GetID();
1613           if(fUseFP) MakeTrack();
1614         }
1615         if(labelpos<0&&labelneg<0) {
1616           FillCandidateSpy("MthNegNeg",true);
1617           if(!fSkipVn)
1618             FillDecayVn("V0SMthNegNegVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
1619         } else if(labelpos>0&&labelneg>0) {
1620           if(!fSkipVn)
1621             FillDecayVn("V0SMthPosPosVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
1622         } else if(labelpos*labelneg<0) {
1623           FillCandidateSpy("MthPosNeg",true);
1624           if(!fSkipVn)
1625             FillDecayVn("V0SMthPosNegVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
1626         }
1627         AliAODVertex *secvtx = myV0->GetSecondaryVtx();
1628         Double_t possec[3],covsec[6];
1629         secvtx->GetXYZ(possec);
1630         secvtx->GetCovarianceMatrix(covsec);
1631         const AliESDVertex vSecVtx(possec,covsec,100.,100);
1632         AliESDtrack trackAtSecI( iT );
1633         trackAtSecI.SetTPCClusterMap( iT->GetTPCClusterMap() );
1634         trackAtSecI.SetTPCSharedMap( iT->GetTPCSharedMap() );
1635         trackAtSecI.SetTPCPointsF( iT->GetTPCNclsF() );
1636         trackAtSecI.PropagateToDCA(&vSecVtx, tAOD->GetMagneticField(), 100);
1637         fDaughterAtSecPhi = trackAtSecI.Phi();
1638         fDaughterAtSecEta = trackAtSecI.Eta();
1639         fDaughterAtSecPt = trackAtSecI.Pt();
1640         LoadTrack(&ieT,iT->Chi2perNDF());
1641         fDaughterMatchPhi=mcpos->Phi();
1642         fDaughterMatchEta=mcpos->Eta();
1643         fDaughterMatchPt=mcpos->Pt();
1644         ieT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1645         fDaughterImpactParameterXY = ip[0];
1646         fDaughterImpactParameterZ = ip[1];
1647         FillTrackSpy("MthDau",true);
1648         if(labelpos<0||labelneg<0) FillTrackSpy("MthNegDau",true);
1649         else FillTrackSpy("MthPosDau",true);
1650         AliESDtrack trackAtSecJ( jT );
1651         trackAtSecJ.SetTPCClusterMap( jT->GetTPCClusterMap() );
1652         trackAtSecJ.SetTPCSharedMap( jT->GetTPCSharedMap() );
1653         trackAtSecJ.SetTPCPointsF( jT->GetTPCNclsF() );
1654         trackAtSecJ.PropagateToDCA(&vSecVtx, tAOD->GetMagneticField(), 100);
1655         fDaughterAtSecPhi = trackAtSecJ.Phi();
1656         fDaughterAtSecEta = trackAtSecJ.Eta();
1657         fDaughterAtSecPt = trackAtSecJ.Pt();
1658         LoadTrack(&jeT,jT->Chi2perNDF());
1659         fDaughterMatchPhi=mcneg->Phi();
1660         fDaughterMatchEta=mcneg->Eta();
1661         fDaughterMatchPt=mcneg->Pt();
1662         jeT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1663         fDaughterImpactParameterXY = ip[0];
1664         fDaughterImpactParameterZ = ip[1];
1665         FillTrackSpy("MthDau",true);
1666         if(labelpos<0||labelneg<0) FillTrackSpy("MthNegDau",true);
1667         else FillTrackSpy("MthPosDau",true);
1668       } else {
1669         FillCandidateSpy("UnMth",false);
1670         if(!fSkipVn)
1671           FillDecayVn("V0SUnMthVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
1672         if(fPostMatched<0) {
1673           fDecayIDneg = iT->GetID();
1674           fDecayIDpos = jT->GetID();
1675           if(fUseFP) MakeTrack();
1676         }
1677         LoadTrack(&ieT,iT->Chi2perNDF());
1678         ieT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1679         fDaughterImpactParameterXY = ip[0];
1680         fDaughterImpactParameterZ = ip[1];
1681         FillTrackSpy("UnMthDau",false);
1682         LoadTrack(&jeT,jT->Chi2perNDF());
1683         jeT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1684         fDaughterImpactParameterXY = ip[0];
1685         fDaughterImpactParameterZ = ip[1];
1686         FillTrackSpy("UnMthDau",false);
1687       }
1688       if(feeddown) {
1689         FillCandidateSpy("MthFeedDown",true);
1690       }
1691     }
1692     //===== END OF MCMATCH
1693   }
1694   ((TH2D*)((TList*)fList->FindObject("V0SAll"))->FindObject("V0SADC"))->Fill( v0all,v0imw );
1695   if(!fSkipVn) {
1696     QCStoreDecayVn("V0SAllVn");
1697     QCStoreDecayVn("V0SSelVn");
1698     if(fReadMC) {
1699       QCStoreDecayVn("V0SMthVn");
1700       QCStoreDecayVn("V0SMthNegNegVn");
1701       QCStoreDecayVn("V0SMthPosPosVn");
1702       QCStoreDecayVn("V0SMthPosNegVn");
1703     }
1704   }
1705   return;
1706 }
1707 //=======================================================================
1708 Bool_t AliAnalysisTaskFlowStrange::PassesPIDCuts(AliESDtrack *myTrack, AliPID::EParticleType pid) {
1709   Bool_t pass=kTRUE;
1710   if(fPIDResponse) {
1711     fDaughterNSigmaPID = fPIDResponse->NumberOfSigmasTPC(myTrack,pid);
1712     if( TMath::Abs(fDaughterNSigmaPID) > fDaughterMaxNSigmaPID )
1713       pass = kFALSE;
1714   }
1715   return pass;
1716 }
1717 //=======================================================================
1718 void AliAnalysisTaskFlowStrange::ChargeParticles(AliAODEvent *tAOD) {
1719   //benchmark purposes
1720   if(!tAOD) return;
1721   TClonesArray* mcArray=NULL;
1722   if(fReadMC) {
1723     mcArray = dynamic_cast<TClonesArray*>(tAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1724     ReadStack(mcArray);
1725   }
1726   for(int i=0; i!=tAOD->GetNumberOfTracks(); ++i) {
1727     AliAODTrack *t = tAOD->GetTrack( i );
1728     if(!t) continue;
1729     if( !t->TestFilterBit(1) ) continue;
1730     fDecayMass=0.0; // using mass as nsigmas control plot
1731     if(fPIDResponse) { // PID
1732       switch(fSpecie) { // TPC PID only
1733       case(kPION):
1734         fDecayMass = fPIDResponse->NumberOfSigmasTPC(t,AliPID::kPion);
1735         break;
1736       case(kKAON):
1737         fDecayMass = fPIDResponse->NumberOfSigmasTPC(t,AliPID::kKaon);
1738         break;
1739       case(kPROTON):
1740         fDecayMass = fPIDResponse->NumberOfSigmasTPC(t,AliPID::kProton);
1741         break;
1742       }
1743     }
1744     Bool_t pass = kTRUE;
1745     if( TMath::Abs(fDecayMass) > 3.0 ) pass=kFALSE;
1746     if( t->Eta()<-0.5 || t->Eta()>+0.5 ) pass=kFALSE;
1747     if( t->Pt()<0.2 || t->Pt()>20.0 ) pass=kFALSE;
1748     AliESDtrack et( t );
1749     et.SetTPCClusterMap( t->GetTPCClusterMap() );
1750     et.SetTPCSharedMap( t->GetTPCSharedMap() );
1751     et.SetTPCPointsF( t->GetTPCNclsF() );
1752     Float_t ip[2];
1753     LoadTrack(&et,t->Chi2perNDF()); 
1754     AliAODVertex *vtx = tAOD->GetPrimaryVertex();
1755     Double_t pos[3];
1756     vtx->GetXYZ(pos);
1757     et.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1758     fDaughterImpactParameterXY = ip[0];
1759     fDaughterImpactParameterZ = ip[1];
1760
1761     FillTrackSpy("TrkAll");
1762     if(!fSkipVn)
1763       FillTrackVn("TrkAllVn",t->Pt(),t->Phi(),t->Eta(),t->GetID());
1764     if(!pass) continue;
1765     FillTrackSpy("TrkSel");
1766     if(!fSkipVn)
1767       FillTrackVn("TrkSelVn",t->Pt(),t->Phi(),t->Eta(),t->GetID());
1768     if(fReadMC) {
1769       ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 1 ); // Selected event 
1770       if(mcArray) {
1771         ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 2 ); // Stack found
1772         bool matched = false;
1773         Int_t label = t->GetLabel();
1774         AliAODMCParticle *mcpar = (AliAODMCParticle*) mcArray->At( TMath::Abs(label) );
1775         if( mcpar ) {
1776           ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 3 ); // Particle in stack
1777           Int_t pdgmcpar = TMath::Abs(mcpar->GetPdgCode());
1778           switch(fSpecie) {
1779           case(kPION):
1780             if(pdgmcpar==211) matched = true;
1781             break;
1782           case(kKAON):
1783             if(pdgmcpar==211) matched = true;
1784             break;
1785           case(kPROTON):
1786             if(pdgmcpar==2212) matched = true;
1787             break;
1788           }
1789           if(!mcpar->IsPrimary()) matched = false;
1790         }
1791         if(matched) {
1792           FillTrackSpy("Mth");
1793           if(!fSkipVn)
1794             FillTrackVn("MthVn",t->Pt(),t->Phi(),t->Eta(),t->GetID());
1795           if(label<0) {
1796             FillTrackSpy("MthNeg");
1797             if(!fSkipVn)
1798               FillTrackVn("MthNegVn",t->Pt(),t->Phi(),t->Eta(),t->GetID());
1799           } else {
1800             FillTrackSpy("MthPos");
1801             if(!fSkipVn)
1802               FillTrackVn("MthPosVn",t->Pt(),t->Phi(),t->Eta(),t->GetID());
1803           }
1804         }
1805       }
1806     }
1807     if(fUseFP) {
1808       fDecayPt=t->Pt();
1809       fDecayPhi=t->Phi();
1810       fDecayEta=t->Eta();
1811       fDecayID=t->GetID();
1812       MakeTrack();
1813     }
1814   }
1815   if(!fSkipVn) {
1816     QCStoreTrackVn("TrkAllVn");
1817     QCStoreTrackVn("TrkSelVn");
1818     if(fReadMC) {
1819       QCStoreTrackVn("MthVn");
1820       QCStoreTrackVn("MthNegVn");
1821       QCStoreTrackVn("MthPosVn");
1822     }
1823   }
1824   return;
1825 }
1826 //=======================================================================
1827 void AliAnalysisTaskFlowStrange::ComputeChi2VZERO() {
1828   Double_t MeanQaQc = ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZEAQmVZEC"))->GetBinContent( 1 );
1829   Double_t MeanQaQa = ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZEASQUARED"))->GetBinContent( 1 );
1830   Double_t MeanQcQc = ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZECSQUARED"))->GetBinContent( 1 );
1831   Double_t MeanQaQt = ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmTPCQmVZEA"))->GetBinContent( 1 );
1832   Double_t MeanQcQt = ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmTPCQmVZEC"))->GetBinContent( 1 );
1833   if(!TMath::AreEqualAbs(MeanQaQt,0,1e-10)&&!TMath::AreEqualAbs(MeanQcQt,0,1e-10)&&!TMath::AreEqualAbs(MeanQaQc,0,1e-10)) {
1834     Double_t OneOverChiSquaredVZEA = MeanQaQa*MeanQcQt/MeanQaQc/MeanQaQt-1;
1835     Double_t OneOverChiSquaredVZEC = MeanQcQc*MeanQaQt/MeanQaQc/MeanQcQt-1;
1836     if(!TMath::AreEqualAbs(OneOverChiSquaredVZEA,0,1e-10)&&!TMath::AreEqualAbs(OneOverChiSquaredVZEC,0,1e-10)) {
1837       ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("ChiSquaredVZEA"))->SetBinContent( 1, 1/OneOverChiSquaredVZEA );
1838       ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("ChiSquaredVZEC"))->SetBinContent( 1, 1/OneOverChiSquaredVZEC );
1839     }
1840   }
1841 }
1842 //=======================================================================
1843 void AliAnalysisTaskFlowStrange::Terminate(Option_t *) {
1844   //terminate
1845   if(fSkipTerminate) return;
1846   ComputeChi2VZERO();
1847   if(fSkipSelection) return;
1848   if(fSkipVn) return;
1849   if(fSpecie<10) {
1850     ComputeDecayVn("V0SAllVn");
1851     ComputeDecayVn("V0SSelVn");
1852     if(fReadMC) {
1853       ComputeDecayVn("V0SMthVn");
1854       ComputeDecayVn("V0SMthPosPosVn");
1855       ComputeDecayVn("V0SMthNegNegVn");
1856       ComputeDecayVn("V0SMthPosNegVn");
1857       ComputeDecayVn("V0SUnMthVn");
1858     }
1859   } else {
1860     ComputeTrackVn("TrkAllVn");
1861     ComputeTrackVn("TrkSelVn");
1862     if(fReadMC) {
1863       ComputeTrackVn("MthVn");
1864       ComputeTrackVn("MthPosVn");
1865       ComputeTrackVn("MthNegVn");
1866     }
1867   }
1868 }
1869 //=======================================================================
1870 void AliAnalysisTaskFlowStrange::MakeTrack() {
1871   // create track for flow tasks
1872   if(fCandidates->GetLast()+5>fCandidates->GetSize()) {
1873     fCandidates->Expand( fCandidates->GetSize()+20 );
1874   }
1875   Bool_t overwrite = kTRUE;
1876   AliFlowCandidateTrack *oTrack = (static_cast<AliFlowCandidateTrack*> (fCandidates->At( fCandidates->GetLast()+1 )));
1877   if( !oTrack ) { // creates new
1878     oTrack = new AliFlowCandidateTrack();
1879     overwrite = kFALSE;
1880   } else { // overwrites
1881     oTrack->ClearMe();
1882   }
1883   oTrack->SetMass(fDecayMass);
1884   oTrack->SetPt(fDecayPt);
1885   oTrack->SetPhi(fDecayPhi);
1886   oTrack->SetEta(fDecayEta);
1887   if(fSpecie<10) {
1888     oTrack->AddDaughter(fDecayIDpos);
1889     oTrack->AddDaughter(fDecayIDneg);
1890   } else {
1891     oTrack->SetID( fDecayID );
1892   }
1893   oTrack->SetForPOISelection(kTRUE);
1894   oTrack->SetForRPSelection(kFALSE);
1895   if(overwrite) {
1896     fCandidates->SetLast( fCandidates->GetLast()+1 );
1897   } else {
1898     fCandidates->AddLast(oTrack);
1899   }
1900   return;
1901 }
1902 //=======================================================================
1903 void AliAnalysisTaskFlowStrange::AddCandidates() {
1904   // adds candidates to flow events (untaging if necessary)
1905   if(fDebug) printf("FlowEventTPC %d tracks | %d RFP | %d POI\n",fTPCevent->NumberOfTracks(),fTPCevent->GetNumberOfRPs(),fTPCevent->GetNumberOfPOIs());
1906   if(fDebug) printf("FlowEventVZE %d tracks | %d RFP | %d POI\n",fVZEevent->NumberOfTracks(),fVZEevent->GetNumberOfRPs(),fVZEevent->GetNumberOfPOIs());
1907   if(fDebug) printf("I received %d candidates\n",fCandidates->GetEntriesFast());
1908   Int_t untagged=0;
1909   Int_t poi=0;
1910   for(int iCand=0; iCand!=fCandidates->GetEntriesFast(); ++iCand ) {
1911     AliFlowCandidateTrack *cand = static_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
1912     if(!cand) continue;
1913     cand->SetForPOISelection(kTRUE);
1914     cand->SetForRPSelection(kFALSE);
1915     poi++;
1916     //if(fDebug) printf(" >Checking at candidate %d with %d daughters: mass %f\n",iCand,cand->GetNDaughters(),cand->Mass());
1917     if(fSpecie<10) { // DECAYS
1918       // untagging ===>
1919       if(fDaughterUnTag) {
1920         for(int iDau=0; iDau!=cand->GetNDaughters(); ++iDau) {
1921           if(fDebug) printf("  >Daughter %d with fID %d", iDau, cand->GetIDDaughter(iDau));
1922           for(int iRPs=0; iRPs!=fTPCevent->NumberOfTracks(); ++iRPs ) {
1923             AliFlowTrack *iRP = static_cast<AliFlowTrack*>(fTPCevent->GetTrack( iRPs ));
1924             if(!iRP) continue;
1925             if(!iRP->InRPSelection()) continue;
1926             if(cand->GetIDDaughter(iDau) == iRP->GetID()) {
1927               if(fDebug) printf(" was in RP set");
1928               ++untagged;
1929               iRP->SetForRPSelection(kFALSE);
1930               fTPCevent->SetNumberOfRPs( fTPCevent->GetNumberOfRPs() -1 );
1931             }
1932           }
1933           if(fDebug) printf("\n");
1934         }
1935       }
1936       // <=== untagging 
1937       fTPCevent->InsertTrack( ((AliFlowTrack*) cand) );
1938     } else {  // CHARGED
1939       // adding only new tracks and tagging accordingly ===>
1940       Bool_t found=kFALSE;
1941       for(int iRPs=0; iRPs!=fTPCevent->NumberOfTracks(); ++iRPs ) {
1942         AliFlowTrack *iRP = static_cast<AliFlowTrack*>(fTPCevent->GetTrack( iRPs ));
1943         if(!iRP) continue;
1944         if(!iRP->InRPSelection()) continue;
1945         if(cand->GetID() == iRP->GetID()) {
1946           if(fDebug) printf("  >charged track (%d) was also found in RP set (adding poi tag)\n",cand->GetID());
1947           iRP->SetMass( cand->Mass() );
1948           iRP->SetForPOISelection(kTRUE);
1949           found = kTRUE;
1950         }
1951       }
1952       if(!found) // not found adding track
1953         fTPCevent->InsertTrack( ((AliFlowTrack*) cand) );
1954     }
1955     fVZEevent->InsertTrack( ((AliFlowTrack*) cand) );
1956   } //END OF LOOP
1957   fTPCevent->SetNumberOfPOIs( poi );
1958   fVZEevent->SetNumberOfPOIs( poi );
1959   ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("POI"))->Fill( poi );
1960   ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("UNTAG"))->Fill( untagged );
1961   if(fDebug) printf("FlowEventTPC %d tracks | %d RFP | %d POI\n",fTPCevent->NumberOfTracks(),fTPCevent->GetNumberOfRPs(),fTPCevent->GetNumberOfPOIs());
1962   if(fDebug) printf("FlowEventVZE %d tracks | %d RFP | %d POI\n",fVZEevent->NumberOfTracks(),fVZEevent->GetNumberOfRPs(),fVZEevent->GetNumberOfPOIs());
1963 }
1964 //=======================================================================
1965 void AliAnalysisTaskFlowStrange::PushBackFlowTrack(AliFlowEvent *flowevent, Double_t pt, Double_t phi, Double_t eta, Double_t w, Int_t id) {
1966   AliFlowTrack rfp;
1967   rfp.SetPt(pt);
1968   rfp.SetPhi(phi);
1969   rfp.SetEta(eta);
1970   rfp.SetWeight(w);
1971   rfp.SetForRPSelection(kTRUE);
1972   rfp.SetForPOISelection(kFALSE);
1973   rfp.SetMass(-999);
1974   rfp.SetID(id);
1975   flowevent->InsertTrack( &rfp );
1976 }
1977 //=======================================================================
1978 Bool_t AliAnalysisTaskFlowStrange::IsAtTPCEdge(Double_t phi,Double_t pt,Int_t charge,Double_t b) {
1979   // Origin: Alex Dobrin
1980   // Implemented by Carlos Perez
1981   TF1 cutLo("cutLo", "-0.01/x+pi/18.0-0.015", 0, 100);
1982   TF1 cutHi("cutHi", "0.55/x/x+pi/18.0+0.03", 0, 100);
1983   Double_t phimod = phi;
1984   if(b<0) phimod = TMath::TwoPi()-phimod;  //for negatve polarity field
1985   if(charge<0) phimod = TMath::TwoPi()-phimod; //for negatve charge
1986   phimod += TMath::Pi()/18.0;
1987   phimod = fmod(phimod, TMath::Pi()/9.0);
1988   if( phimod<cutHi.Eval(pt) && phimod>cutLo.Eval(pt) )
1989     return kTRUE;
1990
1991   return kFALSE;
1992 }
1993 //=======================================================================
1994 void AliAnalysisTaskFlowStrange::MakeQVectors() {
1995   //computes event plane and updates fPsi2
1996   //if there is a problem fPsi->-1
1997   fPsi2=-1;
1998   fVZEWarning=kFALSE;
1999   //=>loading event
2000   MakeQVZE(InputEvent());
2001   MakeQTPC(InputEvent());
2002   if(fUseFP&&fReadMC) {
2003     fVZEevent->SetMCReactionPlaneAngle( fMCEP );    
2004     fTPCevent->SetMCReactionPlaneAngle( fMCEP );    
2005   }
2006   if(fDebug) {
2007     printf("**::MakeQVectors()");
2008     printf("  fQVZEACos %.16f | fQVZEASin %.16f || fQVZEA %.3f | fQVZEC %.3f \n",fQVZEACos, fQVZEASin, fQVZEA, fQVZEC);
2009     printf("  nQTPA_nTracks %d | fQTPC_nTracks %d || fQTPCA %.3f | fQTPCC %.3f \n",fQTPCA_nTracks, fQTPCC_nTracks, fQTPCA, fQTPCC);
2010     printf("  fQTPCACos %.16f | fQTPCASin %.16f || fQTPC2hCos %.16f | fQTPC2hSin %.16f \n",fQTPCACos, fQTPCASin, fQTPC2hCos, fQTPC2hSin);
2011    }
2012   FillMakeQSpy();
2013 }
2014 //=======================================================================
2015 void AliAnalysisTaskFlowStrange::FillMakeQSpy() {
2016   //=>computing psi
2017   //VZERO
2018   Double_t qvzecos,qvzesin,psivzea,psivzec,psivze,qvze, qvzea, qvzec;
2019   psivzea = ( TMath::Pi()+TMath::ATan2(-fQVZEASin,-fQVZEACos) )/fHarmonic;
2020   psivzec = ( TMath::Pi()+TMath::ATan2(-fQVZECSin,-fQVZECCos) )/fHarmonic;
2021   qvzecos = fQVZEACos + fQVZECCos;
2022   qvzesin = fQVZEASin + fQVZECSin;
2023   qvzea = fQVZEA;
2024   qvzec = fQVZEC;
2025   qvze = fQVZEA + fQVZEC;
2026   psivze = ( TMath::Pi()+TMath::ATan2(-qvzesin,-qvzecos) )/fHarmonic;
2027   //TPC
2028   Double_t qtpccos,qtpcsin,psitpca,psitpcc,psitpc,qtpc;
2029   psitpca = ( TMath::Pi()+TMath::ATan2(-fQTPCASin,-fQTPCACos) )/fHarmonic;
2030   psitpcc = ( TMath::Pi()+TMath::ATan2(-fQTPCCSin,-fQTPCCCos) )/fHarmonic;
2031   qtpccos = fQTPCACos + fQTPCCCos;
2032   qtpcsin = fQTPCASin + fQTPCCSin;
2033   qtpc = fQTPCA + fQTPCC;
2034   psitpc = ( TMath::Pi()+TMath::ATan2(-qtpcsin,-qtpccos) )/fHarmonic;
2035   //=>does the event clear?
2036   switch(fWhichPsi) {
2037   case(1): //VZERO
2038     if(fVZEWarning) return;
2039     fPsi2 = psivze;
2040     break;
2041   case(2): //TPC
2042     if(fQTPCA<2||fQTPCC<2) return;
2043     fPsi2 = psitpc;
2044     break;
2045   }
2046   //computing physical Qm vectors
2047   Double_t vzec_qmcos = fQVZECCos/fQVZEC;
2048   Double_t vzec_qmsin = fQVZECSin/fQVZEC;
2049   Double_t vzec_qmnor = TMath::Sqrt( vzec_qmcos*vzec_qmcos + vzec_qmsin*vzec_qmsin );
2050   Double_t vzea_qmcos = fQVZEACos/fQVZEA;
2051   Double_t vzea_qmsin = fQVZEASin/fQVZEA;
2052   Double_t vzea_qmnor = TMath::Sqrt( vzea_qmcos*vzea_qmcos + vzea_qmsin*vzea_qmsin );
2053   Double_t vze_qmcos = qvzecos/qvze;
2054   Double_t vze_qmsin = qvzesin/qvze;
2055   Double_t vze_qmnor = TMath::Sqrt( vze_qmcos*vze_qmcos + vze_qmsin*vze_qmsin );
2056   Double_t tpcc_qmcos = fQTPCCCos/fQTPCC;
2057   Double_t tpcc_qmsin = fQTPCCSin/fQTPCC;
2058   Double_t tpcc_qmnor = TMath::Sqrt( tpcc_qmcos*tpcc_qmcos + tpcc_qmsin*tpcc_qmsin );
2059   Double_t tpca_qmcos = fQTPCACos/fQTPCA;
2060   Double_t tpca_qmsin = fQTPCASin/fQTPCA;
2061   Double_t tpca_qmnor = TMath::Sqrt( tpca_qmcos*tpca_qmcos + tpca_qmsin*tpca_qmsin );
2062   Double_t tpc_qmcos = qtpccos/qtpc;
2063   Double_t tpc_qmsin = qtpcsin/qtpc;
2064   Double_t tpc_qmnor = TMath::Sqrt( tpc_qmcos*tpc_qmcos + tpc_qmsin*tpc_qmsin );
2065   //=>great! recording
2066   ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEPSI"))->Fill( psivze );
2067   ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEPSIA"))->Fill( psivzea );
2068   ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEPSIC"))->Fill( psivzec );
2069   ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("RFPVZE"))->Fill( qvze );
2070   ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZEA"))->Fill( vzea_qmnor*TMath::Sqrt(qvzea) );
2071   ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZEC"))->Fill( vzec_qmnor*TMath::Sqrt(qvzec) );
2072   //------
2073   ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCPSI"))->Fill( psitpc );
2074   ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCPSIA"))->Fill( psitpca );
2075   ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCPSIC"))->Fill( psitpcc );
2076   ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("RFPTPC"))->Fill( qtpc );
2077   ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmTPC"))->Fill( tpc_qmnor*TMath::Sqrt(qtpc) );
2078   //------
2079   ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSI_TPCAVZEC"))->Fill( psitpca, psivzec );
2080   ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSI_TPCCVZEA"))->Fill( psitpcc, psivzea );
2081   ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSI_TPCVZE"))->Fill( psitpc, psivze );
2082
2083   if(fReadMC) {
2084     ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSIMCDIFFTPC"))->Fill( psitpc-fMCEP );
2085     ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSIMCDIFFTPCA"))->Fill( psitpca-fMCEP );
2086     ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSIMCDIFFTPCC"))->Fill( psitpcc-fMCEP );
2087     ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSIMCDIFFVZE"))->Fill( psivze-fMCEP );
2088     ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSIMCDIFFVZEA"))->Fill( psivzea-fMCEP );
2089     ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSIMCDIFFVZEC"))->Fill( psivzec-fMCEP );
2090   }
2091   //------
2092   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQm"))->Fill( 1., tpcc_qmsin, tpcc_qmnor );
2093   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQm"))->Fill( 2., tpcc_qmcos, tpcc_qmnor );
2094   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQm"))->Fill( 3., tpca_qmsin, tpca_qmnor );
2095   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQm"))->Fill( 4., tpca_qmcos, tpca_qmnor );
2096   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQm"))->Fill( 5., tpc_qmsin, tpc_qmnor );
2097   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQm"))->Fill( 6., tpc_qmcos, tpc_qmnor );
2098   //------
2099   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQm"))->Fill( 1., vzec_qmsin, vzec_qmnor );
2100   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQm"))->Fill( 2., vzec_qmcos, vzec_qmnor );
2101   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQm"))->Fill( 3., vzea_qmsin, vzea_qmnor );
2102   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQm"))->Fill( 4., vzea_qmcos, vzea_qmnor );
2103   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQm"))->Fill( 5., vze_qmsin, vze_qmnor );
2104   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQm"))->Fill( 6., vze_qmcos, vze_qmnor );
2105   //------
2106   Double_t vzeqaqc = vzec_qmcos*vzea_qmcos + vzec_qmsin*vzea_qmsin;
2107   Double_t vzeqatpcq = vzea_qmcos*tpc_qmcos + vzea_qmsin*tpc_qmsin;
2108   Double_t vzeqctpcq = vzec_qmcos*tpc_qmcos + vzec_qmsin*tpc_qmsin;
2109   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZEAQmVZEC"))->Fill( 1., vzeqaqc, vze_qmnor );
2110   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZEASQUARED"))->Fill( 1., vzea_qmnor*vzea_qmnor, vze_qmnor );
2111   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZECSQUARED"))->Fill( 1., vzec_qmnor*vzec_qmnor, vze_qmnor );
2112   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmTPCQmVZEA"))->Fill( 1., vzeqatpcq, vze_qmnor );
2113   ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmTPCQmVZEC"))->Fill( 1., vzeqctpcq, vze_qmnor );
2114   return;
2115 }
2116 //=======================================================================
2117 void AliAnalysisTaskFlowStrange::MakeQVZE(AliVEvent *tevent) {
2118   //=>cleaning
2119   if(fUseFP) fVZEevent->ClearFast(); // flowpackage
2120   //=>computing
2121   fQVZEACos=fQVZEASin=fQVZEA=fQVZECCos=fQVZECSin=fQVZEC=0;
2122   Int_t rfp=0;
2123   Double_t eta, phi, w;
2124   //v0c -> qa
2125   for(int id=fVZECa*8;id!=8+fVZECb*8;++id) {
2126     eta = -3.45+0.5*(id/8);
2127     phi = TMath::PiOver4()*(0.5+id%8);
2128     w = tevent->GetVZEROEqMultiplicity(id);
2129     if(w<3) fVZEWarning=kTRUE;
2130     w *= fVZEextW[id];
2131     fQVZECCos += w*TMath::Cos(fHarmonic*phi);
2132     fQVZECSin += w*TMath::Sin(fHarmonic*phi);
2133     fQVZEC += w;
2134     ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEAllPhiEta"))->Fill( phi, eta, w );
2135     rfp++;
2136     if(fUseFP) PushBackFlowTrack(fVZEevent,0,phi,eta,w,0); // flowpackage
2137   }
2138   //v0a -> qb
2139   for(int id=32+fVZEAa*8;id!=40+fVZEAb*8;++id) {
2140     eta = +4.8-0.6*((id/8)-4);
2141     phi = TMath::PiOver4()*(0.5+id%8);
2142     w = tevent->GetVZEROEqMultiplicity(id);
2143     if(w<3) fVZEWarning=kTRUE;
2144     w *= fVZEextW[id];
2145     fQVZEACos += w*TMath::Cos(fHarmonic*phi);
2146     fQVZEASin += w*TMath::Sin(fHarmonic*phi);
2147     fQVZEA += w;
2148     ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEAllPhiEta"))->Fill( phi, eta, w );
2149     rfp++;
2150     if(fUseFP) PushBackFlowTrack(fVZEevent,0,phi,eta,w,0); // flowpackage
2151   }
2152   if(fUseFP) { // flowpackage
2153     fVZEevent->SetNumberOfRPs(rfp);
2154     if(fDebug>0) printf("Inserted tracks in FlowEventVZE %d ==> %.1f\n",fVZEevent->NumberOfTracks(),fQVZEC+fQVZEA);
2155   }
2156 }
2157 //=======================================================================
2158 void AliAnalysisTaskFlowStrange::AddTPCRFPSpy(TList *me) {
2159   TH1D *tH1D;
2160   tH1D = new TH1D("PT",   "PT",           50,0,5);     me->Add(tH1D);
2161   tH1D = new TH1D("PHI",  "PHI", 90,0,TMath::TwoPi()); me->Add(tH1D);
2162   tH1D = new TH1D("ETA",  "ETA",          40,-1,+1);   me->Add(tH1D);
2163   tH1D = new TH1D("TPCS", "TPC Signal",   100,0,500);  me->Add(tH1D);
2164   tH1D = new TH1D("IPXY", "IPXY",         100,-2,+2);  me->Add(tH1D);
2165   tH1D = new TH1D("IPZ",  "IPZ",          100,-2,+2);  me->Add(tH1D);
2166   // TPC
2167   tH1D = new TH1D("TPCNCLS", "NCLS", 170,-0.5,+169.5);   me->Add(tH1D);
2168   tH1D = new TH1D("TPCSHCL", "NSCLS / NCLS", 100,0,1);   me->Add(tH1D);
2169   tH1D = new TH1D("TPCFICL", "NCLS1I / NCLS",100,0,1);   me->Add(tH1D);
2170   tH1D = new TH1D("TPCXRNF", "XROW / NFCLS", 100,0,1.5); me->Add(tH1D);
2171   tH1D = new TH1D("TPCRCHI", "CHI2 / NCLS",  50,0,5);    me->Add(tH1D);
2172   // ITS
2173   tH1D = new TH1D("ITSNCLS", "NCLS",   7,-0.5,+6.5); me->Add(tH1D);
2174   tH1D = new TH1D("ITSRCHI", "CHI2 / NCLS", 50,0,5); me->Add(tH1D);
2175
2176 }
2177 //=======================================================================
2178 Bool_t AliAnalysisTaskFlowStrange::PassesRFPTPCCuts(AliESDtrack *track, Double_t aodchi2cls, Float_t aodipxy, Float_t aodipz) {
2179   if(track->GetKinkIndex(0)>0) return kFALSE;
2180   if( (track->GetStatus()&AliESDtrack::kTPCrefit)==0 ) return kFALSE;
2181   Double_t pt = track->Pt();
2182   Double_t phi = track->Phi();
2183   Double_t eta = track->Eta();
2184   Double_t tpcs = track->GetTPCsignal();
2185   Float_t ipxy, ipz;
2186   track->GetImpactParameters(ipxy,ipz);
2187   Int_t cls = track->GetTPCclusters(0);
2188   Double_t xrows, findcls, chi2;
2189   findcls = track->GetTPCNclsF();
2190   xrows = track->GetTPCCrossedRows();
2191   chi2 = track->GetTPCchi2();
2192   Double_t rchi2 = chi2/cls;
2193   if(!fReadESD) {
2194     rchi2 = aodchi2cls;
2195     ipxy = aodipxy;
2196     ipz = aodipz;
2197   }
2198   Double_t xrnfcls = xrows/findcls;
2199   Double_t scls, cls1i, itschi2;
2200   Int_t itscls;
2201   cls1i = track->GetTPCNclsIter1();
2202   scls = track->GetTPCnclsS();
2203   itscls = track->GetITSclusters(0);
2204   itschi2 = track->GetITSchi2();
2205   Double_t shcl = scls/cls;
2206   Double_t ficl = cls1i/cls;
2207   Double_t itsrchi2 = itscls/itschi2;
2208   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("PT"))->Fill( pt );
2209   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("PHI"))->Fill( phi );
2210   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("ETA"))->Fill( eta );
2211   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCS"))->Fill( tpcs );
2212   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("IPXY"))->Fill( ipxy );
2213   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("IPZ"))->Fill( ipz );
2214   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCNCLS"))->Fill( cls );
2215   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCSHCL"))->Fill( shcl );
2216   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCFICL"))->Fill( ficl );
2217   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCXRNF"))->Fill( xrnfcls );
2218   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCRCHI"))->Fill( rchi2 );
2219   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("ITSNCLS"))->Fill( itscls );
2220   ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("ITSRCHI"))->Fill( itsrchi2 );
2221   if(pt<fRFPminPt) return kFALSE; //0.2
2222   if(pt>fRFPmaxPt) return kFALSE; //5.0
2223   if(eta<fRFPCminEta||(eta>fRFPCmaxEta&&eta<fRFPAminEta)||eta>fRFPAmaxEta) return kFALSE; // -0.8 0.0 0.0 +0.8
2224   if(tpcs<fRFPTPCsignal) return kFALSE; //10.0
2225   if( TMath::Sqrt(ipxy*ipxy/fRFPmaxIPxy/fRFPmaxIPxy+ipz*ipz/fRFPmaxIPz/fRFPmaxIPz)>1 ) return kFALSE; //2.4 3.2
2226   if(cls<fRFPTPCncls) return kFALSE; //70
2227   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("PT"))->Fill( pt );
2228   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("PHI"))->Fill( phi );
2229   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("ETA"))->Fill( eta );
2230   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCS"))->Fill( tpcs );
2231   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("IPXY"))->Fill( ipxy );
2232   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("IPZ"))->Fill( ipz );
2233   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCNCLS"))->Fill( cls );
2234   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCSHCL"))->Fill( shcl );
2235   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCFICL"))->Fill( ficl );
2236   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCXRNF"))->Fill( xrnfcls );
2237   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCRCHI"))->Fill( rchi2 );
2238   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("ITSNCLS"))->Fill( itscls );
2239   ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("ITSRCHI"))->Fill( itsrchi2 );
2240   return kTRUE;
2241 }
2242 //=======================================================================
2243 void AliAnalysisTaskFlowStrange::MakeQTPC(AliVEvent *tevent) {
2244   AliESDEvent *tESD = (AliESDEvent*) (tevent);
2245   AliAODEvent *tAOD = (AliAODEvent*) (tevent);
2246   if(fReadESD) {
2247     if(!tESD) return;
2248     MakeQTPC(tESD);
2249   } else {
2250     if(!tAOD) return;
2251     MakeQTPC(tAOD);
2252   }
2253 }
2254 //=======================================================================
2255 void AliAnalysisTaskFlowStrange::MakeQTPC(AliAODEvent *tAOD) {
2256   //=>cleaning
2257   if(fUseFP) fTPCevent->ClearFast(); // flowpackage
2258   fQTPCACos=fQTPCASin=fQTPCA=fQTPC2hCos=0;
2259   fQTPCCCos=fQTPCCSin=fQTPCC=fQTPC2hSin=0;
2260   fQTPCA_nTracks = 0;
2261   fQTPCC_nTracks = 0;
2262   Int_t rfp=0;
2263   Double_t eta, phi, w;
2264   //=>aod stuff
2265   AliAODVertex *vtx = tAOD->GetPrimaryVertex();
2266   Double_t pos[3],cov[6];
2267   vtx->GetXYZ(pos);
2268   vtx->GetCovarianceMatrix(cov);
2269   const AliESDVertex vESD(pos,cov,100.,100);
2270   AliAODTrack *track;
2271   //=>looping
2272   Int_t rawN = tAOD->GetNumberOfTracks();
2273   for(Int_t id=0; id!=rawN; ++id) {
2274     track = tAOD->GetTrack(id);
2275     //=>cuts
2276     if(!track->TestFilterBit(fRFPFilterBit)) continue;
2277     if( fExcludeTPCEdges )
2278       if( IsAtTPCEdge( track->Phi(), track->Pt(), track->Charge(), tAOD->GetMagneticField() ) ) continue;
2279     AliESDtrack etrack( track );
2280     etrack.SetTPCClusterMap( track->GetTPCClusterMap() );
2281     etrack.SetTPCSharedMap( track->GetTPCSharedMap() );
2282     etrack.SetTPCPointsF( track->GetTPCNclsF() );
2283     Float_t ip[2];
2284     etrack.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
2285     if(!PassesRFPTPCCuts(&etrack,track->Chi2perNDF(),ip[0],ip[1])) continue;
2286     //=>collecting info
2287     phi = track->Phi();
2288     eta = track->Eta();
2289     w = 1;
2290     if(eta<0) {
2291       fQTPCCCos += w*TMath::Cos(fHarmonic*phi);
2292       fQTPCCSin += w*TMath::Sin(fHarmonic*phi);
2293       fQTPCC += w;
2294       fQTPCC_fID[fQTPCC_nTracks++] = track->GetID();
2295     } else {
2296       fQTPCACos += w*TMath::Cos(fHarmonic*phi);
2297       fQTPCASin += w*TMath::Sin(fHarmonic*phi);
2298       fQTPCA += w;
2299       fQTPCA_fID[fQTPCA_nTracks++] = track->GetID();
2300     }
2301     fQTPC2hCos += w*TMath::Cos(2.0*fHarmonic*phi);
2302     fQTPC2hSin += w*TMath::Sin(2.0*fHarmonic*phi);
2303     rfp++;
2304     ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCAllPhiEta"))->Fill( phi, eta, w );
2305     if(fUseFP) PushBackFlowTrack(fTPCevent,track->Pt(),phi,eta,w,track->GetID()); // flow package
2306   }
2307   if(fUseFP) {
2308     fTPCevent->SetNumberOfRPs(rfp);
2309     if(fDebug) printf("Inserted tracks in FlowEventTPC %d ==> %.1f\n",fTPCevent->NumberOfTracks(),fQTPCA+fQTPCC);
2310   }
2311 }
2312 //=======================================================================
2313 void AliAnalysisTaskFlowStrange::MakeQTPC(AliESDEvent *tESD) {
2314   //=>cleaning
2315   if(fUseFP) fTPCevent->ClearFast(); // flow package
2316   fQTPCACos=fQTPCASin=fQTPCA=0;
2317   fQTPCCCos=fQTPCCSin=fQTPCC=0;
2318   fQTPCA_nTracks = 0;
2319   fQTPCC_nTracks = 0;
2320   Int_t rfp=0;
2321   Double_t eta, phi, w;
2322   //=>looping
2323   AliESDtrack *track;
2324   Int_t rawN = tESD->GetNumberOfTracks();
2325   for(Int_t id=0; id!=rawN; ++id) {
2326     track = tESD->GetTrack(id);
2327     //=>cuts
2328     if( fExcludeTPCEdges )
2329       if( IsAtTPCEdge( track->Phi(), track->Pt(), track->Charge(), tESD->GetMagneticField() ) ) continue;
2330     if(!PassesFilterBit(track)) continue;
2331     if(!PassesRFPTPCCuts(track)) continue;
2332     //=>collecting info
2333     phi = track->Phi();
2334     eta = track->Eta();
2335     w = 1;
2336     if(eta<0) {
2337       fQTPCCCos += w*TMath::Cos(fHarmonic*phi);
2338       fQTPCCSin += w*TMath::Sin(fHarmonic*phi);
2339       fQTPCC += w;
2340       fQTPCC_fID[fQTPCC_nTracks++] = track->GetID();
2341     } else {
2342       fQTPCACos += w*TMath::Cos(fHarmonic*phi);
2343       fQTPCASin += w*TMath::Sin(fHarmonic*phi);
2344       fQTPCA += w;
2345       fQTPCA_fID[fQTPCA_nTracks++] = track->GetID();
2346     }
2347     rfp++;
2348     ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCAllPhiEta"))->Fill( phi, eta, w );
2349     if(fUseFP) PushBackFlowTrack(fTPCevent,track->Pt(),phi,eta,w,track->GetID()); // flowpackage
2350   }
2351   if(fUseFP) {
2352     fTPCevent->SetNumberOfRPs(rfp);
2353     if(fDebug) printf("Inserted tracks in FlowEventTPC %d ==> %.1f\n",fTPCevent->NumberOfTracks(),fQTPCA+fQTPCC);
2354   }
2355 }
2356 //=======================================================================
2357 void AliAnalysisTaskFlowStrange::AddMCParticleSpy(TList *me) {
2358   TH1D *tH1D;
2359   TH2D *tH2D;
2360   TProfile *tPro;
2361   tH1D = new TH1D("Pt",   "Pt",   fPtBins,fPtBinEdge);  me->Add(tH1D);
2362   tH1D = new TH1D("Phi",  "Phi",  100,0,TMath::TwoPi()); me->Add(tH1D);
2363   tH1D = new TH1D("Eta",  "Eta",  100,-1,+1);   me->Add(tH1D);
2364   tH1D = new TH1D("Y",    "Y",    100,-1,+1);   me->Add(tH1D);
2365   tH1D = new TH1D("Rad2", "Rad2", 1000,0,+100); me->Add(tH1D);
2366   tH2D = new TH2D("Dphi", "phi-MCEP;pt;dphi",fPtBins,fPtBinEdge, 72,0,TMath::Pi()); me->Add(tH2D);
2367   tPro = new TProfile("Cos2dphi","Cos2dphi",fPtBins,fPtBinEdge); me->Add(tPro);
2368   return;
2369 }
2370 //=======================================================================
2371 void AliAnalysisTaskFlowStrange::FillMCParticleSpy(TString listName, AliAODMCParticle *p) {
2372   ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Pt" ))->Fill( p->Pt() );
2373   ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Eta" ))->Fill( p->Eta() );
2374   ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Y" ))->Fill( p->Y() );
2375   ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Phi" ))->Fill( p->Phi() );
2376   ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Rad2" ))->Fill( TMath::Sqrt( p->Xv()*p->Xv() +
2377                                                                                                  p->Yv()*p->Yv() ) );
2378   ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Dphi" ))->Fill( p->Pt(), GetMCDPHI(p->Phi()) );
2379   ((TProfile*)((TList*)fList->FindObject(listName.Data()))->FindObject("Cos2dphi" ))->Fill( p->Pt(), TMath::Cos( 2*GetMCDPHI(p->Phi()) ), 1 );
2380   return;
2381 }
2382 //=======================================================================
2383 Double_t AliAnalysisTaskFlowStrange::GetMCDPHI(Double_t phi) {
2384   Double_t dDPHI = phi - fMCEP;
2385   //if( dDPHI < 0 ) dDPHI += TMath::TwoPi();
2386   //if( dDPHI > TMath::Pi() ) dDPHI = TMath::TwoPi()-dDPHI;
2387   return dDPHI;
2388 }
2389 //=======================================================================
2390 void AliAnalysisTaskFlowStrange::FillMCParticleSpy(TString listName, TParticle *p) {
2391   ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Pt" ))->Fill( p->Pt() );
2392   ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Eta" ))->Fill( p->Eta() );
2393   ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Phi" ))->Fill( p->Phi() );
2394   ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Rad2" ))->Fill( TMath::Sqrt( p->Vx()*p->Vx() +
2395                                                                                                  p->Vy()*p->Vy() ) );
2396   return;
2397 }
2398 //=======================================================================
2399 void AliAnalysisTaskFlowStrange::AddCandidatesSpy(TList *me,Bool_t res) {
2400   TH1D *tH1D;
2401   TH2D *tH2D;
2402   TProfile *tPro;
2403   TProfile2D *tPro2;
2404   tH2D = new TH2D("PhiEta",  "PhiEta;Phi;Eta", 100,0,TMath::TwoPi(),100,-1,+1);  me->Add(tH2D);
2405   tH2D = new TH2D("PtRAP",   "PtRAP;Pt;Y",     fPtBins,fPtBinEdge,100,-1,+1);  me->Add(tH2D);
2406   tH2D = new TH2D("PtDCA",   "PtDCA;Pt;DCA",   fPtBins,fPtBinEdge,100,0,1.5);  me->Add(tH2D);
2407   tH2D = new TH2D("PtCTP",   "PtCTP;Pt;CTP",   fPtBins,fPtBinEdge,100,0.95,+1);me->Add(tH2D);
2408   tH2D = new TH2D("PtD0D0",  "PtD0D0;Pt;D0D0", fPtBins,fPtBinEdge,100,-5,+5);  me->Add(tH2D);
2409   tH2D = new TH2D("PtRad2",  "PtRad2;Pt;RadXY",fPtBins,fPtBinEdge,100,0,+50);  me->Add(tH2D);
2410   tH2D = new TH2D("PtDL",    "PtDL;Pt;DL",     fPtBins,fPtBinEdge,100,0,+50);  me->Add(tH2D);
2411   tH2D = new TH2D("PtDLlab", "PtDL;Pt;DLlab",  fPtBins,fPtBinEdge,100,0,+100); me->Add(tH2D);
2412   tH2D = new TH2D("PtMASS",  "PtMASS;Pt;MASS", fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); me->Add(tH2D);
2413   tH2D = new TH2D("APPOS",   "APPOS;alphaPOS;QtPOS",100,-2,+2,100,0,0.3);     me->Add(tH2D);
2414   tH2D = new TH2D("D0PD0N",  "D0PD0N;D0P;D0N",      200,-10,+10,200,-10,+10); me->Add(tH2D);
2415   tH2D = new TH2D("XPOSXNEG","XPOSXNEG;XPOS;XNEG",  200,-50,+50,200,-50,+50); me->Add(tH2D);
2416   if(fReadMC) {
2417     if(res) {
2418       tH1D = new TH1D("MCOrigin","MCOrigin;Rad2",   1000,0,50); me->Add(tH1D);
2419       tH2D = new TH2D("PHIRes","PHIRes;PHI;MC-DAT", 72,0,TMath::TwoPi(),100,-0.12,+0.12); me->Add(tH2D);
2420       tH2D = new TH2D("ETARes","ETARes;ETA;MC-DAT", 16,-0.8,+0.8,100,-0.2,+0.2); me->Add(tH2D);
2421       tH2D = new TH2D("PTRes", "PTRes;Pt;MC-DAT",   fPtBins,fPtBinEdge,100,-0.4,+0.4); me->Add(tH2D);
2422       tH2D = new TH2D("RXYRes","RXYRes;RXY;MC-DAT", 100,0,50,100,-4.0,+4.0); me->Add(tH2D);
2423     }
2424     tH2D = new TH2D("PTDPHIMC","PtDPHIMC;Pt;PHI-MCEP",fPtBins,fPtBinEdge,72,0,TMath::Pi()); me->Add(tH2D);
2425     tPro = new TProfile("Cos2dphiMC",  "Cos2dphiMC",fPtBins,fPtBinEdge); me->Add(tPro);
2426     tPro2=new TProfile2D("C2DPHIMCMASS","C2DPHIMCMASS",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); me->Add(tPro2);
2427   }
2428   return;
2429 }
2430 //=======================================================================
2431 void AliAnalysisTaskFlowStrange::FillCandidateSpy(TString listName, Bool_t fillRes) {
2432   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PhiEta"))->Fill( fDecayPhi, fDecayEta );
2433   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtRAP" ))->Fill( fDecayPt, fDecayRapidity );
2434   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtDCA" ))->Fill( fDecayPt, fDecayDCAdaughters );
2435   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtCTP" ))->Fill( fDecayPt, fDecayCosinePointingAngleXY );
2436   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtD0D0"))->Fill( fDecayPt, fDecayProductIPXY );
2437   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtRad2"))->Fill( fDecayPt, fDecayRadXY );
2438   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtDL"  ))->Fill( fDecayPt, fDecayDecayLength );
2439   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtDLlab"))->Fill( fDecayPt, fDecayDecayLengthLab );
2440   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtMASS"))->Fill( fDecayPt, fDecayMass );
2441   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("APPOS" ))->Fill( fDecayAlpha, fDecayQt );
2442   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("D0PD0N"))->Fill( fDecayIPpos, fDecayIPneg );
2443   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("XPOSXNEG"))->Fill( fDecayXpos, fDecayXneg );
2444   if(fReadMC) {
2445     ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PTDPHIMC" ))->Fill( fDecayPt, GetMCDPHI( fDecayPhi ) );
2446     ((TProfile*)((TList*)fList->FindObject(listName.Data()))->FindObject("Cos2dphiMC" ))->Fill( fDecayPt, TMath::Cos( 2*GetMCDPHI(fDecayPhi) ), 1 );
2447     ((TProfile2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("C2DPHIMCMASS" ))->Fill(fDecayPt,fDecayMass,TMath::Cos(2*GetMCDPHI(fDecayPhi)), 1 );
2448     if(fillRes) {
2449       ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("MCOrigin"))->Fill( fDecayMatchOrigin );
2450       ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PHIRes"))->Fill( fDecayPhi, fDecayMatchPhi-fDecayPhi );
2451       ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("ETARes"))->Fill( fDecayEta, fDecayMatchEta-fDecayEta );
2452       ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PTRes"))->Fill( fDecayPt, fDecayMatchPt-fDecayPt );
2453       ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("RXYRes"))->Fill( fDecayRadXY, fDecayMatchRadXY-fDecayRadXY );
2454     }
2455   }
2456 }
2457 //=======================================================================
2458 Bool_t AliAnalysisTaskFlowStrange::AcceptCandidate() {
2459   if(fDecayEta<fDecayMinEta) return kFALSE;
2460   if(fDecayEta>fDecayMaxEta) return kFALSE;
2461   if(fDecayPt<fDecayMinPt) return kFALSE;
2462   if(fDecayProductIPXY>fDecayMaxProductIPXY) return kFALSE;
2463   if(fDecayDCAdaughters>fDecayMaxDCAdaughters) return kFALSE;
2464   if(fDecayCosinePointingAngleXY<fDecayMinCosinePointingAngleXY) return kFALSE;
2465   if(fDecayRadXY<fDecayMinRadXY) return kFALSE;
2466   if(TMath::Abs(fDecayRapidity)>fDecayMaxRapidity) return kFALSE;
2467   if(fSpecie==0) {
2468     if(fDecayAPCutPie) {
2469       if(fDecayQt/TMath::Abs(fDecayAlpha)<fDecayMinQt) return kFALSE;
2470     } else {
2471       if(fDecayQt<fDecayMinQt) return kFALSE;
2472     }
2473     if(fDecayDecayLength>fDecayMaxDecayLength*2.6842) return kFALSE;
2474   } else {
2475     if(fDecayDecayLength>fDecayMaxDecayLength*7.89) return kFALSE;
2476   }
2477   if(fSpecie==1) if(fDecayAlpha>0) return kFALSE;
2478   if(fSpecie==2) if(fDecayAlpha<0) return kFALSE;
2479   return kTRUE;
2480 }
2481 //=======================================================================
2482 void AliAnalysisTaskFlowStrange::AddTrackSpy(TList *me,Bool_t res) {
2483   TH2D *tH2D;
2484   tH2D = new TH2D("PHIETA",       "PHIETA;PHI;ETA",       100,0,TMath::TwoPi(),100,-2,2); me->Add(tH2D);
2485   tH2D = new TH2D("PTTRACKDECAY", "PTTRACKDECAY;PT;PT",   100,0,10,fPtBins,fPtBinEdge); me->Add(tH2D);
2486   tH2D = new TH2D("IPXYIPZ",      "IPXYIPZ;IPXY;IPZ",     100,-10,+10,100,-10,+10); me->Add(tH2D);
2487   tH2D = new TH2D("PTTPCNCLS",    "PTTPCNCLS;PT;NCLS",    fPtBins,fPtBinEdge,170,0,170);  me->Add(tH2D);
2488   tH2D = new TH2D("PTITSNCLS",    "PTITSNCLS;PT;NCLS",    fPtBins,fPtBinEdge,7,-0.5,6.5); me->Add(tH2D);
2489   tH2D = new TH2D("PTITSLAY",     "PTITSLAY;PT;ITSLAYER", fPtBins,fPtBinEdge,6,-0.5,+5.5);me->Add(tH2D);
2490   tH2D = new TH2D("PTITSTPCrefit","PTITSTPCrefit;PT",     fPtBins,fPtBinEdge,2,-0.5,+1.5);me->Add(tH2D);
2491   tH2D->GetYaxis()->SetBinLabel(1,"ITS refit");
2492   tH2D->GetYaxis()->SetBinLabel(2,"TPC refit");
2493   tH2D = new TH2D("POSTPCNCLCHI2","POSTPCNCLCHI2;NCLS;CHI2/NCLS", 170,0,170,100,0,8);   me->Add(tH2D);
2494   tH2D = new TH2D("POSTPCNFCLNXR","POSTPCNFCLNXR;NFCLS;NXR",      170,0,170,170,0,170); me->Add(tH2D);
2495   tH2D = new TH2D("POSTPCNCLNFCL","POSTPCNCLNFCL;NCLS;NFCLS",     170,0,170,170,0,170); me->Add(tH2D);
2496   tH2D = new TH2D("POSTPCNCLNSCL","POSTPCNCLNSCL;NCLS;NSCLS",     170,0,170,170,0,170); me->Add(tH2D);
2497   tH2D = new TH2D("NEGTPCNCLCHI2","NEGTPCNCLCHI2;NCLS;CHI2/NCLS", 170,0,170,100,0,8);   me->Add(tH2D);
2498   tH2D = new TH2D("NEGTPCNFCLNXR","NEGTPCNFCLNXR;NFCLS;NXR",      170,0,170,170,0,170); me->Add(tH2D);
2499   tH2D = new TH2D("NEGTPCNCLNFCL","NEGTPCNCLNFCL;NCLS;NFCLS",     170,0,170,170,0,170); me->Add(tH2D);
2500   tH2D = new TH2D("NEGTPCNCLNSCL","NEGTPCNCLNSCL;NCLS;NSCLS",     170,0,170,170,0,170); me->Add(tH2D);
2501   if(fReadMC) {
2502     TProfile *tPro;
2503     tPro = new TProfile("COSNDPHIMC","COSNDPHIMC",fPtBins,fPtBinEdge); me->Add(tPro);
2504   }
2505   if(res) {
2506     tH2D = new TH2D("PHIRes", "PHIRes;PHI;MC-DAT", 72,0,TMath::TwoPi(),100,-0.12,+0.12); me->Add(tH2D);
2507     tH2D = new TH2D("ETARes", "ETARes;ETA;MC-DAT", 16,-0.8,+0.8,100,-0.2,+0.2); me->Add(tH2D);
2508     tH2D = new TH2D("PTRes",  "PTRes;Pt;MC-DAT", fPtBins,fPtBinEdge,100,-0.4,+0.4); me->Add(tH2D);
2509   }
2510 }
2511 //=======================================================================
2512 void AliAnalysisTaskFlowStrange::FillTrackSpy(TString listName,Bool_t res) {
2513   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PHIETA" ))->Fill( fDaughterPhi, fDaughterEta );
2514   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTTRACKDECAY" ))->Fill( fDaughterPt, fDecayPt );
2515   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "IPXYIPZ" ))->Fill( fDaughterImpactParameterXY, fDaughterImpactParameterZ );
2516   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTTPCNCLS" ))->Fill( fDaughterPt, fDaughterNClsTPC );
2517   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSNCLS" ))->Fill( fDaughterPt, fDaughterNClsITS );
2518   if( TESTBIT(fDaughterITScm,0) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSLAY" ))->Fill( fDaughterPt, 0 );
2519   if( TESTBIT(fDaughterITScm,1) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSLAY" ))->Fill( fDaughterPt, 1 );
2520   if( TESTBIT(fDaughterITScm,2) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSLAY" ))->Fill( fDaughterPt, 2 );
2521   if( TESTBIT(fDaughterITScm,3) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSLAY" ))->Fill( fDaughterPt, 3 );
2522   if( TESTBIT(fDaughterITScm,4) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSLAY" ))->Fill( fDaughterPt, 4 );
2523   if( TESTBIT(fDaughterITScm,5) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSLAY" ))->Fill( fDaughterPt, 5 );
2524   if( (fDaughterStatus&AliESDtrack::kITSrefit) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSTPCrefit" ))->Fill( fDaughterPt, 0 );
2525   if( (fDaughterStatus&AliESDtrack::kTPCrefit) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSTPCrefit" ))->Fill( fDaughterPt, 1 );
2526   TString ch="NEG";
2527   if(fDaughterCharge>0) ch="POS";
2528   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( Form("%sTPCNCLCHI2",ch.Data()) ))->Fill( fDaughterNClsTPC, fDaughterChi2PerNClsTPC );
2529   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( Form("%sTPCNFCLNXR",ch.Data()) ))->Fill( fDaughterNFClsTPC, fDaughterXRows );
2530   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( Form("%sTPCNCLNFCL",ch.Data()) ))->Fill( fDaughterNClsTPC, fDaughterNFClsTPC );
2531   ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( Form("%sTPCNCLNSCL",ch.Data()) ))->Fill( fDaughterNClsTPC, fDaughterNSClsTPC );
2532   if(fReadMC) {
2533     Double_t cosn = TMath::Cos( fHarmonic*GetMCDPHI(fDaughterPhi) );
2534     ((TProfile*)((TList*)fList->FindObject(listName.Data()))->FindObject("COSNDPHIMC" ))->Fill( fDaughterPt, cosn, 1 );
2535   }
2536   if(res) {
2537     ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PHIRes"))->Fill( fDaughterPhi, fDaughterMatchPhi-fDaughterAtSecPhi );
2538     ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("ETARes"))->Fill( fDaughterEta, fDaughterMatchEta-fDaughterAtSecEta );
2539     ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PTRes"))->Fill( fDaughterPt, fDaughterMatchPt-fDaughterAtSecPt );
2540   }
2541 }
2542 //=======================================================================
2543 void AliAnalysisTaskFlowStrange::LoadTrack(AliESDtrack *myTrack, Double_t aodChi2NDF) {
2544   fDaughterCharge = myTrack->Charge();
2545   fDaughterXRows = myTrack->GetTPCCrossedRows();
2546   fDaughterNFClsTPC = myTrack->GetTPCNclsF();
2547   fDaughterNSClsTPC = myTrack->GetTPCnclsS();
2548   fDaughterNClsTPC = myTrack->GetTPCclusters(0);
2549   if(fReadESD) {
2550     if(fDaughterNClsTPC>0) fDaughterChi2PerNClsTPC = myTrack->GetTPCchi2()/fDaughterNClsTPC;
2551   } else {
2552     fDaughterChi2PerNClsTPC = aodChi2NDF;
2553   }
2554   myTrack->GetImpactParameters(fDaughterImpactParameterXY,fDaughterImpactParameterZ);
2555   fDaughterStatus = myTrack->GetStatus();
2556   fDaughterITScm = myTrack->GetITSClusterMap();
2557   fDaughterPhi = myTrack->Phi();
2558   fDaughterEta = myTrack->Eta();
2559   fDaughterPt = myTrack->Pt();
2560   fDaughterKinkIndex = myTrack->GetKinkIndex(0);
2561   fDaughterNClsITS=0;
2562   for(Int_t lay=0; lay!=6; ++lay)
2563     if(TESTBIT(fDaughterITScm,lay)) fDaughterNClsITS++;
2564 }
2565 //=======================================================================
2566 Bool_t AliAnalysisTaskFlowStrange::AcceptDaughter(Bool_t strongITS) {
2567   if(fDaughterKinkIndex>0) return kFALSE;
2568   if( (fDaughterStatus&AliESDtrack::kTPCrefit)==0 ) return kFALSE;
2569   if(fDaughterNFClsTPC<1) return kFALSE;
2570   if(fDaughterPt<fDaughterMinPt) return kFALSE;
2571   if(fDaughterEta<fDaughterMinEta) return kFALSE;
2572   if(fDaughterEta>fDaughterMaxEta) return kFALSE;
2573   if(fDaughterNClsTPC<fDaughterMinNClsTPC) return kFALSE;
2574   if(fDaughterXRows<fDaughterMinXRows) return kFALSE;
2575   if(fDaughterChi2PerNClsTPC>fDaughterMaxChi2PerNClsTPC) return kFALSE;
2576   if(TMath::Abs(fDaughterImpactParameterXY)<fDaughterMinImpactParameterXY) return kFALSE;
2577   if(fDaughterXRows<fDaughterMinXRowsOverNClsFTPC*fDaughterNFClsTPC) return kFALSE;
2578   if(strongITS) {
2579     if( (fDaughterITSrefit) & ((fDaughterStatus&AliESDtrack::kITSrefit)==0) ) return kFALSE;
2580     for(Int_t lay=0; lay!=6; ++lay)
2581       if(fDaughterITSConfig[lay]>-0.5) {
2582         if(fDaughterITSConfig[lay]) {
2583           if(!TESTBIT(fDaughterITScm,lay)) return kFALSE;
2584         } else {
2585           if(TESTBIT(fDaughterITScm,lay)) return kFALSE;
2586         }
2587       }
2588     if(fDaughterNClsITS<fDaughterMinNClsITS) return kFALSE;
2589     if(fDaughterSPDRequireAny) {
2590       if( !TESTBIT(fDaughterITScm,0)&&!TESTBIT(fDaughterITScm,1)) return kFALSE;
2591     }
2592   }
2593   return kTRUE;
2594 }
2595 //=======================================================================
2596 Double_t AliAnalysisTaskFlowStrange::GetWDist(const AliVVertex* v0, const AliVVertex* v1) {
2597   // calculate sqrt of weighted distance to other vertex
2598   if (!v0 || !v1) {
2599     printf("One of vertices is not valid\n");
2600     return 0;
2601   }
2602   static TMatrixDSym vVb(3);
2603   double dist = -1;
2604   double dx = v0->GetX()-v1->GetX();
2605   double dy = v0->GetY()-v1->GetY();
2606   double dz = v0->GetZ()-v1->GetZ();
2607   double cov0[6],cov1[6];
2608   v0->GetCovarianceMatrix(cov0);
2609   v1->GetCovarianceMatrix(cov1);
2610   vVb(0,0) = cov0[0]+cov1[0];
2611   vVb(1,1) = cov0[2]+cov1[2];
2612   vVb(2,2) = cov0[5]+cov1[5];
2613   vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
2614   vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
2615   vVb.InvertFast();
2616   if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
2617   dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
2618     +    2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
2619   return dist>0 ? TMath::Sqrt(dist) : -1;
2620 }
2621 //=======================================================================
2622 Bool_t AliAnalysisTaskFlowStrange::plpMV(const AliVEvent *event) {
2623   // check for multi-vertexer pile-up
2624   const AliAODEvent *aod = (const AliAODEvent*)event;
2625   const AliESDEvent *esd = (const AliESDEvent*)event;
2626   //
2627   const int    kMinPlpContrib = 5;
2628   const double kMaxPlpChi2 = 5.0;
2629   const double kMinWDist = 15;
2630   //
2631   if (!aod && !esd) {
2632     printf("Event is neither of AOD nor ESD\n");
2633     exit(1);
2634   }
2635   //
2636   const AliVVertex* vtPrm = 0;
2637   const AliVVertex* vtPlp = 0;
2638   int nPlp = 0;
2639   //
2640   if (aod) {
2641     if ( !(nPlp=aod->GetNumberOfPileupVerticesTracks()) ) return kFALSE;
2642     vtPrm = aod->GetPrimaryVertex();
2643     if (vtPrm == aod->GetPrimaryVertexSPD()) return kTRUE; // there are pile-up vertices but no primary
2644   }
2645   else {
2646     if ( !(nPlp=esd->GetNumberOfPileupVerticesTracks())) return kFALSE;
2647     vtPrm = esd->GetPrimaryVertexTracks();
2648     if (((AliESDVertex*)vtPrm)->GetStatus()!=1) return kTRUE; // there are pile-up vertices but no primary
2649   }
2650   //int bcPrim = vtPrm->GetBC();
2651   //
2652   for (int ipl=0;ipl<nPlp;ipl++) {
2653     vtPlp = aod ? (const AliVVertex*)aod->GetPileupVertexTracks(ipl) : (const AliVVertex*)esd->GetPileupVertexTracks(ipl);
2654     //
2655     if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
2656     if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
2657     //  int bcPlp = vtPlp->GetBC();
2658     //  if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other BC
2659     //
2660     double wDst = GetWDist(vtPrm,vtPlp);
2661     if (wDst<kMinWDist) continue;
2662     //
2663     return kTRUE; // pile-up: well separated vertices
2664   }
2665   //
2666   return kFALSE;
2667   //
2668 }
2669 //=======================================================================
2670 void AliAnalysisTaskFlowStrange::MakeFilterBits() {
2671   //FilterBit 1
2672   fFB1 = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
2673   //FilterBit1024
2674   fFB1024 = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE,1);
2675   fFB1024->SetMinNCrossedRowsTPC(120);
2676   fFB1024->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
2677   fFB1024->SetMaxChi2PerClusterITS(36);
2678   fFB1024->SetMaxFractionSharedTPCClusters(0.4);
2679   fFB1024->SetMaxChi2TPCConstrainedGlobal(36);
2680   fFB1024->SetEtaRange(-0.9,0.9);
2681   fFB1024->SetPtRange(0.15, 1e10);
2682 }
2683 //=======================================================================
2684 Bool_t AliAnalysisTaskFlowStrange::PassesFilterBit(AliESDtrack *track) {
2685   Bool_t ret=kFALSE;
2686   switch(fRFPFilterBit) {
2687     case(1024):
2688       ret = fFB1024->AcceptTrack(track);
2689       break;
2690     default:
2691       ret = fFB1->AcceptTrack(track);
2692   }
2693   return ret;
2694 }
2695 //=======================================================================
2696 void AliAnalysisTaskFlowStrange::LoadVZEROResponse() {
2697   if(fVZEResponse) {
2698     TString run = fVZEResponse->GetTitle();
2699     if( run.Atoi() == fRunNumber ) return;
2700     fVZEResponse = NULL;
2701   }
2702   //==>loading
2703   fVZEResponse = dynamic_cast<TH2D*> (fVZEload->FindObject( Form("%d",fRunNumber) ));
2704   if(fVZEResponse) {
2705     printf("New VZE calibration: run %d || %s -> Entries %.0f\n",fRunNumber, fVZEResponse->GetTitle(),fVZEResponse->GetEntries());
2706     //=>external weights
2707     for(int i=0;i!=64;++i) fVZEextW[i]=1;
2708     if(!fVZEsave) {
2709       Double_t minC = fCentPerMin, maxC = fCentPerMax;
2710       if(fVZEmb) {
2711         minC = 0;
2712         maxC = 80;
2713       }
2714       Int_t ybinmin = fVZEResponse->GetYaxis()->FindBin(minC+1e-6);
2715       Int_t ybinmax = fVZEResponse->GetYaxis()->FindBin(maxC-1e-6);
2716       if(fSkipCentralitySelection) {
2717         ybinmin=-1;
2718         ybinmax=-1;
2719       }
2720       for(int i=0;i!=64;++i) fVZEextW[i] = fVZEResponse->Integral(i+1,i+1,ybinmin,ybinmax)/(maxC-minC);
2721       //ring-wise normalization
2722       Double_t ring[8];
2723       for(int j=0; j!=8; ++j) {
2724         ring[j]=0;
2725         for(int i=0;i!=8;++i) ring[j] += fVZEextW[j*8+i]/8;
2726       }
2727       //disk-wise normalization
2728       Double_t disk[2];
2729       int xbinmin, xbinmax;
2730       xbinmin = 1+8*fVZECa;
2731       xbinmax = 8+8*fVZECb;
2732       disk[0] = fVZEResponse->Integral(xbinmin,xbinmax,ybinmin,ybinmax)/(maxC-minC)/(xbinmax-xbinmin+1);
2733       xbinmin = 33+8*fVZEAa;
2734       xbinmax = 40+8*fVZEAb;
2735       disk[1] = fVZEResponse->Integral(xbinmin,xbinmax,ybinmin,ybinmax)/(maxC-minC)/(xbinmax-xbinmin+1);
2736       //for(int i=0;i!=64;++i) printf("CELL %d -> W = %f ||",i,fVZEextW[i]);
2737       if(fVZEByDisk) {
2738         for(int i=0;i!=64;++i) fVZEextW[i] = disk[i/32]/fVZEextW[i];
2739       } else {
2740         for(int i=0;i!=64;++i) fVZEextW[i] = ring[i/8]/fVZEextW[i];
2741       }
2742       //for(int i=0;i!=64;++i) printf(" W = %f \n",fVZEextW[i]);
2743     }
2744   } else {
2745     printf("VZE calibration: requested but not found!!!\n");
2746   }
2747 }
2748 //=======================================================================
2749 void AliAnalysisTaskFlowStrange::AddVZEQA() {
2750   TProfile2D *prof;
2751   TH2D *tH2D;
2752   TList *tList = new TList();
2753   tList->SetName( "VZEQA" );
2754   tList->SetOwner();
2755   fList->Add( tList );
2756   tH2D = new TH2D("EQU","EQU;VZEeqmult-VZEmult;cell",100,-5,+5,64,0,64); tList->Add( tH2D );
2757   prof = new TProfile2D("LINbefCAL","LINbef;VZEcell;VZEeqmult;SPDtrkl", 64,0,64,350,0,700,0,10000); tList->Add( prof );
2758   prof = new TProfile2D("LINaftCAL","LINaft;VZEcell;VZEeqmult;SPDtrkl", 64,0,64,350,0,700,0,10000); tList->Add( prof );
2759 }
2760 //=======================================================================
2761 void AliAnalysisTaskFlowStrange::FillVZEQA() {
2762   AliESDEvent *tESD = (AliESDEvent*) (InputEvent());
2763   AliAODEvent *tAOD = (AliAODEvent*) (InputEvent());
2764   if(fReadESD) {
2765     if(!tESD) return;
2766     //FillVZEQA(tESD);
2767   } else {
2768     if(!tAOD) return;
2769     FillVZEQA(tAOD);
2770   }
2771 }
2772 //=======================================================================
2773 void AliAnalysisTaskFlowStrange::FillVZEQA(AliAODEvent *tAOD) {
2774   AliVVZERO *vzero = tAOD->GetVZEROData();
2775   AliAODTracklets *tracklets = tAOD->GetTracklets();
2776   if(!vzero) return;
2777   if(!tracklets) return;
2778   Double_t mult, eqmult;
2779   Int_t trkl = tracklets->GetNumberOfTracklets();
2780   for(int id=0; id!=64; ++id) {
2781     mult = vzero->GetMultiplicity(id);
2782     eqmult = tAOD->GetVZEROEqMultiplicity(id);
2783     ((TH2D*) ((TList*) fList->FindObject("VZEQA"))->FindObject("EQU"))->Fill(eqmult-mult,id);
2784     ((TProfile2D*) ((TList*) fList->FindObject("VZEQA"))->FindObject( "LINbefCAL" ))->Fill(id,eqmult,trkl,1);
2785     ((TProfile2D*) ((TList*) fList->FindObject("VZEQA"))->FindObject( "LINaftCAL" ))->Fill(id,eqmult*fVZEextW[id],trkl,1);
2786   }
2787 }
2788 //=======================================================================
2789 void AliAnalysisTaskFlowStrange::AddVZEROResponse() {
2790   fVZEResponse = NULL;
2791   AliVEvent *event = InputEvent();
2792   if(!event) return;
2793   Int_t thisrun = event->GetRunNumber();
2794   fVZEResponse = new TH2D( Form("%d",thisrun), Form("%d;cell;CC",thisrun), 64,0,64, 110, -10, 100);
2795   fList->Add(fVZEResponse);
2796 }
2797 //=======================================================================
2798 void AliAnalysisTaskFlowStrange::SaveVZEROResponse() {
2799   if(!fVZEResponse) return;
2800   AliVEvent *event = InputEvent();
2801   if(!event) return;
2802   Double_t w;
2803   // reject event with low ocupancy in VZERO
2804   // for centralities below 60% this should not happen anyway
2805   Double_t rejectEvent = kFALSE;
2806   for(int id=0; id!=64; ++id) {
2807     w = event->GetVZEROEqMultiplicity(id);
2808     if(w<3) rejectEvent = kTRUE;
2809   }
2810   if(rejectEvent) return;
2811   // saves weights
2812   for(int id=0; id!=64; ++id) {
2813     w = event->GetVZEROEqMultiplicity(id);
2814     fVZEResponse->Fill(id,fThisCent,w);
2815   }
2816 }
2817 //=======================================================================
2818 Int_t AliAnalysisTaskFlowStrange::RefMult(AliAODEvent *tAOD, Int_t fb) {
2819   Int_t found = 0;
2820   for(int i=0; i!=tAOD->GetNumberOfTracks(); ++i) {
2821     AliAODTrack *t = tAOD->GetTrack( i );
2822     if(!t) continue;
2823     if( !t->TestFilterBit(fb) ) continue;
2824     if( t->Eta()<-0.8 || t->Eta()>+0.8 ) continue;
2825     if( t->Pt()<0.2 || t->Pt()>5.0 ) continue;
2826     if( t->GetTPCNcls()<70 ) continue;
2827     //if( t->GetTPCsignal()<10.0 ) continue;
2828     if( t->Chi2perNDF()<0.2 ) continue;
2829     ++found;
2830   }
2831   return found;
2832 }
2833 //=======================================================================
2834 Int_t AliAnalysisTaskFlowStrange::RefMultTPC() {
2835   AliAODEvent *ev = dynamic_cast<AliAODEvent*> (InputEvent());
2836   if(!ev) return -1;
2837   Int_t found = 0;
2838   for(int i=0; i!=ev->GetNumberOfTracks(); ++i) {
2839     AliAODTrack *t = ev->GetTrack( i );
2840     if(!t) continue;
2841     if( !t->TestFilterBit(1) ) continue;
2842     if( t->Eta()<-0.8 || t->Eta()>+0.8 ) continue;
2843     if( t->Pt()<0.2 || t->Pt()>5.0 ) continue;
2844     if( t->GetTPCNcls()<70 ) continue;
2845     if( t->GetTPCsignal()<10.0 ) continue;
2846     if( t->Chi2perNDF()<0.2 ) continue;    
2847     ++found;
2848   }
2849   return found;
2850 }
2851 //=======================================================================
2852 Int_t AliAnalysisTaskFlowStrange::RefMultGlobal() {
2853   AliAODEvent *ev = dynamic_cast<AliAODEvent*> (InputEvent());
2854   if(!ev) return -1;
2855   Int_t found = 0;
2856   for(int i=0; i!=ev->GetNumberOfTracks(); ++i) {
2857     AliAODTrack *t = ev->GetTrack( i );
2858     if(!t) continue;
2859     if( !t->TestFilterBit(16) ) continue;
2860     if( t->Eta()<-0.8 || t->Eta()>+0.8 ) continue;
2861     if( t->Pt()<0.2 || t->Pt()>5.0 ) continue;
2862     if( t->GetTPCNcls()<70 ) continue;
2863     if( t->GetTPCsignal()<10.0 ) continue;
2864     if( t->Chi2perNDF()<0.1 ) continue;    
2865     Double_t b[3], bcov[3];
2866     if( !t->PropagateToDCA(ev->GetPrimaryVertex(),ev->GetMagneticField(),100,b,bcov) ) continue;
2867     if( b[0]>+0.3 || b[0]<-0.3 || b[1]>+0.3 || b[1]<-0.3) continue;
2868     ++found;
2869   }
2870   return found;
2871 }
2872 //=======================================================================
2873 void AliAnalysisTaskFlowStrange::ResetContainers() {
2874   if(fUseFP) {
2875     fTPCevent->ClearFast();
2876     fVZEevent->ClearFast();
2877   }
2878 }
2879 //=======================================================================
2880 void AliAnalysisTaskFlowStrange::AddTrackVn(TList *tList) {
2881   TProfile *tProfile;
2882   TH1D *tH1D;
2883   // vze
2884   tProfile = new TProfile("SP_uVZEA","u x Q_{VZEA}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2885   tProfile = new TProfile("SP_uVZEC","u x Q_{VZEC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2886   tProfile = new TProfile("SP_VZEAVZEC","Q_{VZEA} x Q_{VZEC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2887   tProfile = new TProfile("SP_VZEATPC","Q_{VZEA} x Q_{TPC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2888   tProfile = new TProfile("SP_VZECTPC","Q_{VZEC} x Q_{TPC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2889   // error
2890   tProfile = new TProfile("SP_uVZEAuVZEC","u x Q_{VZEA} . u x Q_{VZEC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2891   tProfile = new TProfile("SP_uVZEAVZEAVZEC","u x Q_{VZEA} . Q_{VZEA} x Q_{VZEC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2892   tProfile = new TProfile("SP_uVZECVZEAVZEC","u x Q_{VZEC} . Q_{VZEA} x Q_{VZEC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2893   // tpc
2894   tProfile = new TProfile("SP_uTPCA","u x Q_{TPCA}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2895   tProfile = new TProfile("SP_uTPCC","u x Q_{TPCC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2896   tProfile = new TProfile("SP_TPCATPCC","Q_{TPCA} x Q_{TPCC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2897   // error
2898   tProfile = new TProfile("SP_uTPCATPCATPCC","u x Q_{TPCA} . Q_{TPCA} x Q_{TPCC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2899   tProfile = new TProfile("SP_uTPCCTPCATPCC","u x Q_{TPCC} . Q_{TPCA} x Q_{TPCC}",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2900   // control
2901   tH1D = new TH1D("QC_HistPt_P","HistPt_P",fPtBins,fPtBinEdge); tList->Add( tH1D );
2902   tH1D = new TH1D("QC_HistPt_Q","HistPt_Q",fPtBins,fPtBinEdge); tList->Add( tH1D );
2903   // qc
2904   tProfile = new TProfile("QC_C2","QC_C2",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2905   tProfile = new TProfile("QC_C4","QC_C4",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2906   tProfile = new TProfile("QC_DC2","QC_DC2",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2907   tProfile = new TProfile("QC_DC4","QC_DC4",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2908   tProfile = new TProfile("QC_C2C4","QC_C2C4",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2909   tProfile = new TProfile("QC_C2DC2","QC_C2DC2",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2910   tProfile = new TProfile("QC_C2DC4","QC_C2DC4",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2911   tProfile = new TProfile("QC_C4DC2","QC_C4DC2",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2912   tProfile = new TProfile("QC_C4DC4","QC_C4DC4",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2913   tProfile = new TProfile("QC_DC2DC4","QC_DC2DC4",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2914   // qc transient
2915   tProfile = new TProfile("QC_pCos","QC_pCos",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2916   tProfile = new TProfile("QC_pSin","QC_pSin",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2917   tProfile = new TProfile("QC_qCos","QC_qCos",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2918   tProfile = new TProfile("QC_qSin","QC_qSin",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2919   tProfile = new TProfile("QC_q2hCos","QC_q2hCos",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2920   tProfile = new TProfile("QC_q2hSin","QC_q2hSin",fPtBins,fPtBinEdge,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2921   // measurements
2922   tH1D = new TH1D("SP_vnVZEA","SP_vnVZEA",fPtBins,fPtBinEdge); tList->Add( tH1D );
2923   tH1D = new TH1D("SP_vnVZEC","SP_vnVZEC",fPtBins,fPtBinEdge); tList->Add( tH1D );
2924   tH1D = new TH1D("SP_vnTPCA","SP_vnTPCA",fPtBins,fPtBinEdge); tList->Add( tH1D );
2925   tH1D = new TH1D("SP_vnTPCC","SP_vnTPCC",fPtBins,fPtBinEdge); tList->Add( tH1D );
2926   tH1D = new TH1D("QC_Cum2","QC_Cum2",fPtBins,fPtBinEdge); tList->Add( tH1D );
2927   tH1D = new TH1D("QC_Cum4","QC_Cum4",fPtBins,fPtBinEdge); tList->Add( tH1D );
2928   tH1D = new TH1D("QC_DCum2","QC_DCum2",fPtBins,fPtBinEdge); tList->Add( tH1D );
2929   tH1D = new TH1D("QC_DCum4","QC_DCum4",fPtBins,fPtBinEdge); tList->Add( tH1D );
2930   tH1D = new TH1D("SP_vnVZEGA","SP_vnVZEGA",fPtBins,fPtBinEdge); tList->Add( tH1D );
2931   tH1D = new TH1D("SP_vnVZEWA","SP_vnVZEWA",fPtBins,fPtBinEdge); tList->Add( tH1D );
2932   tH1D = new TH1D("SP_vnTPCAA","SP_vnTPCAA",fPtBins,fPtBinEdge); tList->Add( tH1D );
2933   tH1D = new TH1D("QC_vn2","QC_vn2",fPtBins,fPtBinEdge); tList->Add( tH1D );
2934   tH1D = new TH1D("QC_vn4","QC_vn4",fPtBins,fPtBinEdge); tList->Add( tH1D );
2935   if(fReadMC) {
2936     TH2D *tH2D;
2937     tProfile = new TProfile("MC_COSNDPHI","MC_COSNDPHI",fPtBins,fPtBinEdge,-3,+3); tList->Add( tProfile );
2938     tH2D = new TH2D("MC_COSNDPHI_uQVZEA","MC_COSNDPHI_uQVZEA",100,-1,+1,100,-0.3,+0.3); tList->Add( tH2D );
2939     tH2D = new TH2D("MC_COSNDPHI_uQVZEC","MC_COSNDPHI_uQVZEC",100,-1,+1,100,-0.3,+0.3); tList->Add( tH2D );
2940     tH2D = new TH2D("MC_COSNDPHI_uQTPCA","MC_COSNDPHI_uQTPCA",100,-1,+1,100,-0.3,+0.3); tList->Add( tH2D );
2941     tH2D = new TH2D("MC_COSNDPHI_uQTPCC","MC_COSNDPHI_uQTPCC",100,-1,+1,100,-0.3,+0.3); tList->Add( tH2D );
2942   }
2943 }
2944 //=======================================================================
2945 void AliAnalysisTaskFlowStrange::AddDecayVn(TList *tList) {
2946   TProfile2D *tProfile;
2947   TH2D *tH2D;
2948   // decay
2949   tH2D = new TH2D("DecayYield_PtMass","Decay_PtMass",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2950   tProfile = new TProfile2D("DecayAvgPt_PtMass","Decay_PtMass",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tProfile ); tProfile->Sumw2();
2951   // vze
2952   tProfile = new TProfile2D("SP_uVZEA","u x Q_{VZEA}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2953   tProfile = new TProfile2D("SP_uVZEC","u x Q_{VZEC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2954   tProfile = new TProfile2D("SP_VZEAVZEC","Q_{VZEA} x Q_{VZEC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2955   tProfile = new TProfile2D("SP_VZEATPC","Q_{VZEA} x Q_{TPC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2956   tProfile = new TProfile2D("SP_VZECTPC","Q_{VZEC} x Q_{TPC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2957   // error
2958   tProfile = new TProfile2D("SP_uVZEAuVZEC","u x Q_{VZEA} . u x Q_{VZEC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2959   tProfile = new TProfile2D("SP_uVZEAVZEAVZEC","u x Q_{VZEA} . Q_{VZEA} x Q_{VZEC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2960   tProfile = new TProfile2D("SP_uVZECVZEAVZEC","u x Q_{VZEC} . Q_{VZEA} x Q_{VZEC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2961   // tpc
2962   tProfile = new TProfile2D("SP_uTPCA","u x Q_{TPCA}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2963   tProfile = new TProfile2D("SP_uTPCC","u x Q_{TPCC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2964   tProfile = new TProfile2D("SP_TPCATPCC","Q_{TPCA} x Q_{TPCC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2965   // error
2966   tProfile = new TProfile2D("SP_uTPCATPCATPCC","u x Q_{TPCA} . Q_{TPCA} x Q_{TPCC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2967   tProfile = new TProfile2D("SP_uTPCCTPCATPCC","u x Q_{TPCC} . Q_{TPCA} x Q_{TPCC}",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2968   // control
2969   tH2D = new TH2D("QC_HistPt_P","HistPt_P",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2970   tH2D = new TH2D("QC_HistPt_Q","HistPt_Q",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2971   // qc
2972   tProfile = new TProfile2D("QC_C2","QC_C2",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2973   tProfile = new TProfile2D("QC_C4","QC_C4",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2974   tProfile = new TProfile2D("QC_DC2","QC_DC2",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2975   tProfile = new TProfile2D("QC_DC4","QC_DC4",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2976   tProfile = new TProfile2D("QC_C2C4","QC_C2C4",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2977   tProfile = new TProfile2D("QC_C2DC2","QC_C2DC2",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2978   tProfile = new TProfile2D("QC_C2DC4","QC_C2DC4",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2979   tProfile = new TProfile2D("QC_C4DC2","QC_C4DC2",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2980   tProfile = new TProfile2D("QC_C4DC4","QC_C4DC4",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2981   tProfile = new TProfile2D("QC_DC2DC4","QC_DC2DC4",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2982   // qc transient
2983   tProfile = new TProfile2D("QC_pCos","QC_pCos",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2984   tProfile = new TProfile2D("QC_pSin","QC_pSin",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2985   tProfile = new TProfile2D("QC_qCos","QC_qCos",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2986   tProfile = new TProfile2D("QC_qSin","QC_qSin",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2987   tProfile = new TProfile2D("QC_q2hCos","QC_q2hCos",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2988   tProfile = new TProfile2D("QC_q2hSin","QC_q2hSin",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2989   // measurements
2990   tH2D = new TH2D("SP_vnVZEA","SP_vnVZEA",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2991   tH2D = new TH2D("SP_vnVZEC","SP_vnVZEC",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2992   tH2D = new TH2D("SP_vnTPCA","SP_vnTPCA",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2993   tH2D = new TH2D("SP_vnTPCC","SP_vnTPCC",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2994   tH2D = new TH2D("QC_Cum2","QC_Cum2",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2995   tH2D = new TH2D("QC_Cum4","QC_Cum4",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2996   tH2D = new TH2D("QC_DCum2","QC_DCum2",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2997   tH2D = new TH2D("QC_DCum4","QC_DCum4",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2998   tH2D = new TH2D("SP_vnVZEGA","SP_vnVZEGA",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
2999   tH2D = new TH2D("SP_vnVZEWA","SP_vnVZEWA",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
3000   tH2D = new TH2D("SP_vnTPCAA","SP_vnTPCAA",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
3001   tH2D = new TH2D("QC_vn2","QC_vn2",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
3002   tH2D = new TH2D("QC_vn4","QC_vn4",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tH2D );
3003   if(fReadMC) {
3004     tProfile = new TProfile2D("MC_COSNDPHI","MC_COSNDPHI",fPtBins,fPtBinEdge,fMassBins,fMinMass,fMaxMass); tList->Add( tProfile );
3005   }
3006 }
3007 //=======================================================================
3008 TList* AliAnalysisTaskFlowStrange::RebinDecayVn(Int_t nbins, Int_t *bins) {
3009   TList *out = new TList();
3010   out->SetOwner();
3011   out->SetName( fList->GetName() );
3012
3013   TList *ori, *end;
3014
3015   ori = (TList*) fList->FindObject("Event");
3016   out->Add( ori );
3017
3018   ori = (TList*) fList->FindObject("MakeQSpy");
3019   out->Add( ori );
3020
3021   ori = (TList*) fList->FindObject("V0SAllVn");
3022   end = RebinDecayVn(ori,nbins,bins);
3023   end->SetName( ori->GetName() ); out->Add( end );
3024
3025   ori = (TList*) fList->FindObject("V0SSelVn");
3026   end = RebinDecayVn(ori,nbins,bins);
3027   end->SetName( ori->GetName() ); out->Add( end );
3028
3029   if(fReadMC) {
3030     ori = (TList*) fList->FindObject("V0SMthVn");
3031     end = RebinDecayVn(ori,nbins,bins);
3032     end->SetName( ori->GetName() ); out->Add( end );
3033
3034     ori = (TList*) fList->FindObject("V0SMthPosPosVn");
3035     end = RebinDecayVn(ori,nbins,bins);
3036     end->SetName( ori->GetName() ); out->Add( end );
3037
3038     ori = (TList*) fList->FindObject("V0SMthNegNegVn");
3039     end = RebinDecayVn(ori,nbins,bins);
3040     end->SetName( ori->GetName() ); out->Add( end );
3041
3042     ori = (TList*) fList->FindObject("V0SMthPosNegVn");
3043     end = RebinDecayVn(ori,nbins,bins);
3044     end->SetName( ori->GetName() ); out->Add( end );
3045
3046     ori = (TList*) fList->FindObject("V0SUnMthVn");
3047     end = RebinDecayVn(ori,nbins,bins);
3048     end->SetName( ori->GetName() ); out->Add( end );
3049   }
3050
3051   return out;
3052 }
3053 //=======================================================================
3054 TList* AliAnalysisTaskFlowStrange::RebinDecayVn(TList *tList,Int_t nbins, Int_t *bins) {
3055   // getting expected number of mass bins
3056   int sum=0;
3057   for(int i=0; i!=nbins; ++i) sum += bins[i];
3058
3059   TList *list = new TList();
3060   list->SetOwner();
3061
3062   Int_t npt;
3063   Double_t pt[200], mass[200];
3064   TH2D *tH2D;
3065
3066   tH2D = ((TH2D*)tList->FindObject( "DecayYield_PtMass" ));
3067   list->Add(tH2D); //keeping it as it is
3068   int nmassbins = tH2D->GetNbinsY();
3069   //consistency check
3070   if( nmassbins!=sum  ) {
3071     printf("Error: incompatible binning %d vs %d\nBYE\n",nmassbins,sum);
3072     return NULL;
3073   }
3074   //reading pts
3075   npt = tH2D->GetNbinsX();
3076   for(int i=0; i!=npt+1; ++i) pt[i] = tH2D->GetXaxis()->GetBinLowEdge(i+1);
3077   //making mass bins
3078   for(int i=0,j=0; i!=nbins+1; ++i) {
3079     mass[i] = tH2D->GetYaxis()->GetBinLowEdge(j+1);
3080     j += bins[i];
3081   }
3082   //TProfile2D migrating info
3083   TProfile2D *tProfileOld, *tProfileNew;
3084   TString tprofiles[31] = {"DecayAvgPt_PtMass","SP_uVZEA","SP_uVZEC","SP_VZEAVZEC","SP_VZEATPC",
3085                            "SP_VZECTPC","SP_uVZEAuVZEC","SP_uVZEAVZEAVZEC","SP_uVZECVZEAVZEC","SP_uTPCA",
3086                            "SP_uTPCC","SP_TPCATPCC","SP_uTPCATPCATPCC","SP_uTPCCTPCATPCC","QC_C2",
3087                            "QC_C4","QC_DC2","QC_DC4","QC_C2C4","QC_C2DC2",
3088                            "QC_C2DC4","QC_C4DC2","QC_C4DC4","QC_DC2DC4","QC_pCos",
3089                            "QC_pSin","QC_qCos","QC_qSin","QC_q2hCos","QC_q2hSin",
3090                            "MC_COSNDPHI"};
3091   for(int i=0; i!=31; ++i) {
3092     tProfileOld = (TProfile2D*) tList->FindObject( tprofiles[i].Data() );
3093     if(!tProfileOld) continue;
3094     TArrayD *oldsw2 = tProfileOld->GetSumw2();
3095     TArrayD *oldbsw2 = tProfileOld->GetBinSumw2();
3096     tProfileNew = new TProfile2D( tprofiles[i].Data(), tprofiles[i].Data(),
3097                                   npt,pt,nbins,mass,"s");
3098     tProfileNew->Sumw2();
3099     list->Add( tProfileNew );
3100     TArrayD *newsw2 = tProfileNew->GetSumw2();
3101     TArrayD *newbsw2 = tProfileNew->GetBinSumw2();
3102     for(int x=1; x!=tProfileOld->GetNbinsX()+1; ++x) { //pt
3103       for(int y=1; y!=tProfileNew->GetNbinsY()+1; ++y) { //mass
3104         Double_t sContent=0;
3105         Double_t sEntries=0;
3106         Double_t sSqWeigh=0;
3107         Double_t sSqBWeig=0;
3108         Int_t binnew = tProfileNew->GetBin(x,y);
3109         Double_t minmass = tProfileNew->GetYaxis()->GetBinLowEdge( y );
3110         Double_t maxmass = tProfileNew->GetYaxis()->GetBinLowEdge( y+1 );
3111         Int_t minbin = tProfileOld->GetYaxis()->FindBin( minmass+1e-10 );
3112         Int_t maxbin = tProfileOld->GetYaxis()->FindBin( maxmass-1e-10 );
3113         for(int k=minbin; k!=maxbin+1; ++k) {
3114           Int_t binold = tProfileOld->GetBin(x,k);
3115           Double_t wk = tProfileOld->GetBinEntries( binold );
3116           Double_t yk = tProfileOld->GetBinContent( binold );
3117           sEntries += wk;
3118           sContent += wk*yk;
3119           sSqWeigh += oldsw2->At(binold);
3120           if(oldbsw2->GetSize()) sSqBWeig += oldbsw2->At(binold);
3121         }
3122         tProfileNew->SetBinEntries( binnew, sEntries );
3123         tProfileNew->SetBinContent( binnew, sContent );
3124         newsw2->SetAt(sSqWeigh, binnew);
3125         if(oldbsw2->GetSize()) newbsw2->SetAt(sSqBWeig, binnew);
3126       }
3127     }
3128     tProfileOld = NULL;
3129   }
3130   //TH2D all new blank
3131   TString th2ds[15] = {"QC_HisPt_P","QC_HistPt_Q","SP_vnVZEA","SP_vnVZEC","SP_vnTPCA",
3132                        "SP_vnTPCC","QC_Cum2","QC_Cum4","QC_DCum2","QC_DCum4",
3133                        "SP_vnVZEGA","SP_vnVZEWA","SP_vnTPCAA","QC_vn2","QC_vn4"};
3134   for(int i=0; i!=15; ++i) {
3135     tH2D = new TH2D( th2ds[i].Data(),th2ds[i].Data(),
3136                      npt,pt,nbins,mass);
3137     list->Add(tH2D);
3138   }
3139   return list;
3140
3141   //int nmass = tH2Dold->GetNbinsY();
3142
3143   /*
3144   // decay
3145   TString nam[38] = {
3146   for(int i=0; i!=38; ++i) {
3147     tProfile = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( nam[i].Data() ));
3148     if(i==o) { // binning
3149       nxs = tProfile->GetXaxis()->GetNbins();
3150       nys = tProfile->GetYaxis()->GetNbins();
3151       if(pt) {
3152         nxs = nbins;
3153         *nx = *bins;
3154         for(int y=0; y!=nys; ++y)
3155           ny[y] = tProfile->GetYaxis()->GetBinLowEdge(y+1);
3156       } else {
3157         nys = nbins;
3158         *ny = *bins;
3159         for(int x=0; x!=nxs; ++x)
3160           nx[x] = tProfile->GetXaxis()->GetBinLowEdge(x+1);
3161       }
3162     }
3163     tProfileNew = new TProfile2D(tProfile->GetName(),tProfile->GetTitle(),nxs,nx,nys,ny,"s"); list->Add( tProfileNew );
3164     if(pt) {
3165       for(int y=0; y!=nys; ++y) {
3166         
3167         for(int x=0; x!=nxs; ++x) {
3168         }
3169       }
3170     } else {
3171
3172     }
3173   }
3174   */
3175 }
3176 //=======================================================================
3177 void AliAnalysisTaskFlowStrange::QCStoreTrackVn(TString name) {
3178   // getting transients
3179   TProfile *pCos = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pCos" ));
3180   TProfile *pSin = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pSin" ));
3181   TProfile *qCos = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qCos" ));
3182   TProfile *qSin = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qSin" ));
3183   TProfile *q2hCos = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hCos" ));
3184   TProfile *q2hSin = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hSin" ));
3185   if(fDebug) {
3186     printf("**QCStoreTrackVn( %s )\n",name.Data());
3187     printf(" Read %.0f entries in p set and %.0f entries in q set\n",pCos->GetEntries(),qCos->GetEntries());
3188     printf(" pCos[5] %.16f | pSin[5] %.16f \n", pCos->GetBinContent(5), pSin->GetBinContent(5));
3189     printf(" qCos[5] %.16f | qSin[5] %.16f \n", qCos->GetBinContent(5), qSin->GetBinContent(5));
3190 &n