changes for Vertex and Tracks classes
[u/mrichter/AliRoot.git] / PWGCF / EBYE / LRC / AliAnalysisTaskLRC.cxx
1 /**************************************************************************
2  * Author: Andrey Ivanov.                                           *
3  * Contributors are mentioned in the code where appropriate.              *
4  *                                                                        *
5  * Permission to use, copy, modify and distribute this software and its   *
6  * documentation strictly for non-commercial purposes is hereby granted   *
7  * without fee, provided that the above copyright notice appears in all   *
8  * copies and that both the copyright notice and this permission notice   *
9  * appear in the supporting documentation. The authors make no claims     *
10  * about the suitability of this software for any purpose. It is          *
11  * provided "as is" without express or implied warranty.                  *
12  **************************************************************************/
13 // Analysis task for Long Range Correlation (LRC) analysis using TPC data
14 // This includes a TList of AliLRCBase objects that are processing LRC analysis
15 // for a given Eta window 
16 // This task is worcking with ESD data only
17
18 // Authors : Andrey Ivanov , Igor Altsybeev, St.Peterburg State University
19 // Email: Andrey.Ivanov@cern.ch, Igor.Altsybeev@cern.ch
20
21 #include <AliAnalysisManager.h>
22 #include <AliESDInputHandler.h>
23 #include "AliAnalysisTaskLRC.h"
24 #include <AliLRCBase.h>
25 //#include <AliLRCProcess.h>
26 #include <AliVEvent.h>
27 #include <AliMCEvent.h>
28 #include <AliESDEvent.h>     
29 #include "AliAODEvent.h"
30 #include "AliAODTrack.h"
31 #include "AliAODInputHandler.h"
32
33 #include "AliGenEventHeader.h"
34 #include "AliGenHijingEventHeader.h"
35
36 #include <AliPID.h>
37 #include <AliPIDCombined.h>
38 #include <AliPIDResponse.h>
39 #include <AliMultiplicity.h>
40 #include <AliESDv0.h>
41
42 #include <AliRunLoader.h>
43 #include <AliRun.h>
44 #include <AliStack.h>
45 #include <AliESDtrackCuts.h>
46 #include <AliKalmanTrack.h>
47 #include <AliPhysicsSelection.h>
48 #include <AliCentrality.h>
49 #include <AliESDZDC.h>
50 #include <AliESDFMD.h>
51
52 #include "AliEventplane.h"
53
54
55 #include <TParticle.h>
56 #include <TH1.h>
57 #include <TH2.h>
58 #include "TH1I.h"
59 #include "TRandom3.h"
60 #include <TComplex.h>
61
62 #include <TF1.h>
63
64 //#include "AliSimpleEvent.h"
65
66 //flag use my Event Tree
67 //#define fSetIncludeEventTreeInOutput 1
68
69
70 #include "TStopwatch.h"
71
72 using std::endl;
73 using std::cout;
74
75 ClassImp(AliAnalysisTaskLRC)
76
77 //________________________________________________________________________
78 AliAnalysisTaskLRC::AliAnalysisTaskLRC( const char *name, Bool_t runKine)
79     : AliAnalysisTaskSE(name)
80     ,fNumberOfPhiSectors(1)
81     ,fFlagSuppressAddingSomeHistos(kTRUE)
82     ,fAnalysisLevel("ESD")
83     ,fEsdTrackCuts(0)
84     ,fAODtrackCutBit(128)
85     ,fArrTrackCuts(0x0)
86     ,fHistCutsNamesBins(0)
87     //    ,fNumberOfCutsToRemember(0)
88     ,fSwitchToListingCuts(0)
89     
90     //    ,fAnalysisType(en_AnalysisType_ESD)
91     ,fMinNumberOfSPDtracklets(0)
92     ,fMaxPtLimit(100.0)
93     ,fMinPtLimit(0.0)
94     ,fKineLowPtCut(0.0)
95     ,fMinAcceptedTracksCut(0)
96     ,fMaxAcceptedTracksCut(0)
97     ,fCheckForkVtx(kTRUE)
98     ,fCheckForVtxPosition(kFALSE)
99     ,fVxMax(0)
100     ,fVyMax(0)
101     ,fVzMax(0)
102     ,fLRCproc(0)
103     ,fOutList(0)
104     ,fRunKine(0)
105     ,fAnalyseMCESD(0)
106     ,fShowEventStats(kFALSE)
107     ,fShowPerTrackStats(kFALSE)
108     ,fEtaRegionForTests(0.8)
109     ,fMultCutInEtaRegion(0)
110     ,fHistEventCutStats(0)
111     ,fHistTrackCutStats(0)
112     ,fHistAODTrackStats(0)
113     ,fHistVx(0)
114     ,fHistVy(0)
115     ,fHistVz(0)
116     ,fHistVxMCrecoDiff(0)
117     ,fHistVyMCrecoDiff(0)
118     ,fHistVzMCrecoDiff(0)
119     ,fHistVertexNconributors(0)
120     ,fHistNumberOfPileupVerticesTracks(0)
121     ,fHistNumberOfPileupVerticesSPD(0)
122     ,fHistEventPlane(0)
123     ,fHistPt(0)
124     ,fHistEta(0)
125     ,fHistEtaAODpure(0)
126     ,fHistPhi(0)
127     ,fHistEtaPhi(0)
128     ,fHistPhiLRCrotationsCheck(0)
129     ,fHistPhiArtificialProfilingCheck(0)
130     ,fHistPhiArtificialProfilingCheckWrtEvPlane(0)
131     ,fHistPhiArtificialEvPlane(0)
132     ,fHistEtaVsZvCoverage(0)
133     ,fHistEtaVsZvCoverageAccepted(0)
134     ,fHistMultBeforeCuts(0)
135     ,fHistAcceptedMult(0)
136     ,fHistAcceptedTracks(0)
137     ,fHistMultiplicityInEtaRegion(0)
138     ,fHistMultiplicityInEtaRegionAfterPtCuts(0)
139     ,fHist2DMultiplicityMCESDInEtaRegion(0)
140     ,fHistAcceptedTracksAfterPtCuts(0)
141     ,fHistAcceptedTPCtracks(0)
142     ,fHistClustersTPC(0)
143     ,fHistClustersTPCafterCuts(0)
144     ,fHistCrossedRowsTPC(0)
145     ,fHistCrossedRowsTPCafterCuts(0)
146     ,fHistClustersITS(0)
147     ,fHistTrackletsITS(0)
148     ,fHist2DClustersTPCvsPt(0)
149     ,fHist2DClustersTPCvsEta(0)
150     ,fHist2DAcceptedTracksPtvsEta(0)
151     ,fHistMClabels(0)
152     ,fHistRejectedTracksCharge(0)
153     ,fHistTracksCharge(0)
154     ,fHistPtPlus(0)
155     ,fHistPtMinus(0)
156     ,fHistNetChargeVsPt(0)
157     ,fHistChargePlusVsPtTmp(0)
158     ,fHistChargeMinusVsPtTmp(0)
159     ,fHist2DNetChargeVsPt(0)
160     ,fHist2DNetChargeVsPtCorrectedOnEventMean(0)
161     ,fHist2DNetChargeVsPtCorrectedOnEventMeanNormOnNch(0)
162     ,fHistProbabilitiesPID(0)
163     ,fHistESDtrackMass(0)
164     ,fHistProbabilityPion(0)
165     ,fHistProbabilityKaon(0)
166     ,fHistProbabilityProton(0)
167     ,fHistParticlesDistr(0)
168     ,fHistParticlesDistrBeforeCuts(0)
169     //    ,fNumberOfSectors(1)
170     ,fHistCentralityPercentile(0)
171     ,fHistCentralityClass10(0)
172     ,fHistCentralityClass5(0)
173     //,fHistZDCenergy(0x0)
174     ,fHistZDCparticipants(0)
175     ,fHistV0multiplicity(0)
176     ,fHistV0Amultiplicity(0)
177     ,fHistV0Cmultiplicity(0)
178     ,fHist2DV0ACmultiplicity(0)
179     ,fHist2DTracksAcceptedVsV0multiplicity(0)
180     //,fHistV0spectra(0)
181     
182     ,fHistV0cells    (0)
183     ,fHistV0Acells   (0)
184     ,fHistV0Ccells   (0)
185     ,fHist2DV0ACcells(0)
186     
187     ,fThresholdOnV0mult(0)
188     
189     
190     ,fMinCentralityClass(-0.01)
191     ,fMaxCentralityClass(98)
192     ,fIsIonsAnalysis(kFALSE)
193     ,fEtInsteadOfPt(kFALSE)
194     ,fUsePhiShufflingByHand(kFALSE)
195     ,fUseToyEvents(kFALSE)
196     ,fNumberOfToyEvents(1)
197     ,fTmpCounter(0)
198     ,fPIDResponse(0x0)
199     ,fPIDCombined(0x0)
200     
201     ,fHistPidMaxProbability(0)
202     ,fHistPidPureMaxProbability(0)
203     ,fStrPIDforFwd("any")
204     ,fStrPIDforBwd("any")
205
206     ,fPIDsensingFlag(kFALSE)
207     ,fArtificialInefficiency(-1.)
208     ,fHistNumberOfDroppedByHandTracks(0)
209     ,fRand(0)
210     
211     
212     ,fPhiArtificialGapBegin(0)
213     ,fPhiArtificialGapEnd(0)
214     
215     ,fFlagWatchZDC(0)
216     ,fFlagWatchV0 (0)
217     ,fFlagWatchFMD(0)
218     
219     ,fAnalysisTimer(0)
220     
221     //MC qa and studies
222     //    ,fHistMCvertexRdeltaFromParent(0)
223     //    ,fHistMCparentsStat(0)
224     //,fHistMCparentsEta(0)
225     //,fHistMCchildsEta(0)
226     //,fHistMCdeltaEtaChildParent(0)
227     //,fProbabilitiesPID(0)
228     
229     
230     //    ,fSimpleEvent(0)
231     //    ,fNsimpleEvents(0)
232     //    ,fEventTree(0)
233     //    ,fSetIncludeEventTreeInOutput(0)
234 {
235     fAnalysisTimer = new TStopwatch;
236     //Init
237     fRunKine = runKine;
238     for ( int i = 0; i < 5; i++ )
239     {
240         fHistZDCenergy[i] = 0x0;
241     }
242     for ( int i = 0; i < 4; i++ )
243     {
244         fHistV0AmultiplicityRing[i] = 0x0;
245         fHistV0CmultiplicityRing[i] = 0x0;
246         fHist2DV0ACmultiplicityRing[i] = 0x0;
247         fHist2DTracksAcceptedVsV0AmultiplicityRing[i] = 0x0;
248         fHist2DTracksAcceptedVsV0CmultiplicityRing[i] = 0x0;
249     }
250     //fProbabilitiesPID = new Double_t[AliPID::kSPECIES];
251     // Output slot #1 writes into a TList container for common task data and QA
252     DefineOutput(1, TList::Class());
253     
254     //Defining output slots for each LRC processor (required to avoid TList of TLists on merging)
255     for(Int_t i=0; i < fLRCproc.GetEntries(); i++)
256     {
257         DefineOutput(Proc(i)->GetOutputSlotNumber(),TList::Class());
258     }
259 }
260
261
262 // ---------------------------------------  Setters ------------------
263
264 void AliAnalysisTaskLRC::SetMaxPtLimit(Double_t MaxPtLimit)
265 {
266     //Sets  Max Pt filter
267     fMaxPtLimit = MaxPtLimit;
268 }
269 void AliAnalysisTaskLRC::SetMinPtLimit(Double_t MinPtLimit)
270 {
271     //Sets  Min Pt filter
272     fMinPtLimit = MinPtLimit;
273 }
274 void AliAnalysisTaskLRC::SetKineLowPtCut(Double_t kineLowPtCut)
275 {
276     //Sets  Min Pt filter for Kine tracks
277     fKineLowPtCut = kineLowPtCut;
278 }
279 AliLRCBase*  AliAnalysisTaskLRC::Proc(Int_t index)
280 {
281     // Get Processor i
282     return (dynamic_cast<AliLRCBase*> (fLRCproc.At(index)));
283 }
284
285 void AliAnalysisTaskLRC::SetMinNumberOfSPDtracklets( Int_t MinSPDtracklets )
286 {
287     //Sets  Min SPD tracklets number
288     fMinNumberOfSPDtracklets = MinSPDtracklets;
289 }
290
291 //________________________________________________________________________
292 void AliAnalysisTaskLRC::UserCreateOutputObjects()
293 {
294     // --------- Output list
295     fOutList = new TList();
296     
297     //### added 13.12.2011
298     fOutList->SetOwner();  // IMPORTANT!
299     
300     fAnalysisTimer->Start();
301     
302     
303     //##### PID stuff
304     // ------- setup PIDCombined
305     fPIDCombined=new AliPIDCombined;
306     fPIDCombined->SetDefaultTPCPriors();
307     fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
308     
309     // no light nuclei
310     fPIDCombined->SetSelectedSpecies(AliPID::kSPECIES);
311     //Int_t lNumberOfPidSpecies = (Int_t)AliPID::kSPECIES;
312     for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec)
313     {
314         fProbTPCTOF[ispec]=new TH2D(Form("prob%s_mom_TPCTOF",AliPID::ParticleName(ispec)),
315                                     Form("%s probability vs. momentum;momentum;probability",AliPID::ParticleName(ispec)),
316                                     100,0.,20.,50,0.,1.);
317         fProbAllDets[ispec]=new TH2D(Form("prob%s_mom_AllDets",AliPID::ParticleName(ispec)),
318                                      Form("%s probability vs. momentum;momentum;probability",AliPID::ParticleName(ispec)),
319                                      100,0.,20.,50,0.,1.);
320         if ( !fFlagSuppressAddingSomeHistos )
321         {
322             fOutList->Add(fProbTPCTOF[ispec]);
323             fOutList->Add(fProbAllDets[ispec]);
324         }
325         
326         
327         // basic priors
328         fPriors[ispec] = new TH1F(Form("%s_priors",AliPID::ParticleName(ispec)),
329                                   Form("%s priors vs momentum",AliPID::ParticleName(ispec)),
330                                   100,0.,20.);
331         if ( !fFlagSuppressAddingSomeHistos )
332             fOutList->Add(fPriors[ispec]);
333         switch (ispec) {
334         case AliPID::kElectron:
335             for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.02);
336             break;
337         case AliPID::kMuon:
338             for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.02);
339             break;
340         case AliPID::kPion:
341             for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.56);
342             break;
343         case AliPID::kKaon:
344             for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.20);
345             break;
346         case AliPID::kProton:
347             for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.20);
348             break;
349         default:
350             break;
351         }
352         fPIDCombined->SetPriorDistribution((AliPID::EParticleType)ispec,fPriors[ispec]);
353         
354         // priors used
355         fPriorsUsed[ispec] = new TH2D(Form("%s_priorsUsed",AliPID::ParticleName(ispec)),
356                                       Form("%s priors vs transverse momentum;p_{t} (GeV/c);priors",AliPID::ParticleName(ispec)),
357                                       100,0.,20.,101,0,1.01);
358         if ( !fFlagSuppressAddingSomeHistos )
359             fOutList->Add(fPriorsUsed[ispec]);
360     }
361     
362     fHistPidMaxProbability = new TH1D("fHistPidMaxProbability","HistPidMaxProbability;Probability;Entries",100,0,1);
363     fOutList->Add(fHistPidMaxProbability);
364     
365     fHistPidPureMaxProbability = new TH1D("fHistPidPureMaxProbability","HistPidMaxProbability;Probability;Entries",100,0,1);
366     fOutList->Add(fHistPidPureMaxProbability);
367     
368     // ##### end PID part
369     Printf("UserCreateOutputObjects.........");
370     //Disabling "replacing existing TH2D ...etc" warning
371     
372     Bool_t lTH1oldStatus = TH1::AddDirectoryStatus();
373     TH1::AddDirectory(kFALSE);
374     
375     Bool_t lTH2oldStatus = TH2::AddDirectoryStatus();
376     TH2::AddDirectory(kFALSE);
377     
378     
379     
380     //LRC processors init
381     Int_t lLrcNum = fLRCproc.GetEntries();
382     for(Int_t i = 0; i < lLrcNum; i++)
383     {
384         AliLRCBase *lrcBase = dynamic_cast<AliLRCBase*> (fLRCproc.At(i));
385         if(lrcBase)
386             lrcBase->InitDataMembers();
387         else continue;
388         //remember pointer to be used in the analysis
389         //fLRCprocArrayPointers[i] = lrcBase;
390     }
391     
392     
393     //fOutList->Add( fEsdTrackCuts );
394     
395     
396     //create here array with bin names
397     //15.03.2013 - list with cuts to apply on tracks and remember decisions
398     if ( fSwitchToListingCuts )
399     {
400         int lNumberOfCutsToRemember = fArrTrackCuts.GetEntries();
401         fHistCutsNamesBins = new TH1I("fHistCutsNamesBins","Tracks in different cut conditions;;N_{tracks}", lNumberOfCutsToRemember,-0.5,lNumberOfCutsToRemember-0.5);
402         for(Int_t i = 1; i <= lNumberOfCutsToRemember; i++)
403         {
404             //            TString cutsName = (TString*)fArrCutsNames[i-1];
405             //            fHistCutsNamesBins->GetXaxis()->SetBinLabel(i,((TString*)fArrCutsNames[i-1])->Data());
406             fHistCutsNamesBins->GetXaxis()->SetBinLabel(i,fArrCutsNames[i-1].Data());
407             
408         }
409         fOutList->Add(fHistCutsNamesBins);
410     }
411     
412     
413     //Event level QA
414     const Int_t nEventStatBins = 10;
415     fHistEventCutStats = new TH1I("fHistEventCutStats","Event statistics;;N_{events}", nEventStatBins,-0.5,nEventStatBins-0.5);
416     TString gEventCutBinNames[nEventStatBins] = {"Total","No trigger","Wrong centrality","No reco vertex","Bad vertex params", "Bad vertex position","Few SPD tracklets","HighMult cut","LowMult cut","Analyzed"};
417     for(Int_t i = 1; i <= nEventStatBins; i++)fHistEventCutStats->GetXaxis()->SetBinLabel(i,gEventCutBinNames[i-1].Data());
418     fOutList->Add(fHistEventCutStats);
419     
420     //Track level QA
421     const Int_t nTracksStatBins = 6;
422     fHistTrackCutStats = new TH1I("fHistTrackCutStats","Track statistics;;N_{tracks}", nTracksStatBins,-0.5,nTracksStatBins-0.5);
423     TString gTrackCutBinNames[nTracksStatBins] = {"Total","QA track cuts","HighPt cut","LowPt cut","Good", "Not phys primary"};
424     for(Int_t i = 1; i <= nTracksStatBins; i++)fHistTrackCutStats->GetXaxis()->SetBinLabel(i,gTrackCutBinNames[i-1].Data());
425     fOutList->Add(fHistTrackCutStats);
426     
427     //AOD track cut bits stats
428     const int nAODtrackStats = 16;
429     fHistAODTrackStats = new TH1I("fHistAODTrackStats","AOD tracks statistics;TrackFilterBit;N_{tracks}",nAODtrackStats,-0.5,nAODtrackStats-0.5);
430     if ( fAnalysisLevel == "AOD" )
431         fOutList->Add(fHistAODTrackStats);
432     
433     
434     //Vertex distributions
435     fHistVx = new TH1D("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",200,-0.4,0.4);
436     fOutList->Add(fHistVx);
437     fHistVy = new TH1D("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",200,-0.4,0.4);
438     fOutList->Add(fHistVy);
439     fHistVz = new TH1D("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",200,-20.,20.);
440     fOutList->Add(fHistVz);
441     
442     if ( fRunKine )
443     {
444         fHistVxMCrecoDiff = new TH1D("fHistVxMCrecoDiff","Primary vertex mc-reco diff - x coordinate;V_{x} (cm);Entries",200,-0.2,0.2);
445         fOutList->Add(fHistVxMCrecoDiff);
446         fHistVyMCrecoDiff = new TH1D("fHistVyMCrecoDiff","Primary vertex mc-reco diff - y coordinate;V_{y} (cm);Entries",200,-0.2,0.2);
447         fOutList->Add(fHistVyMCrecoDiff);
448         fHistVzMCrecoDiff = new TH1D("fHistVMCrecoDiffz","Primary vertex mc-reco diff - z coordinate;V_{z} (cm);Entries",200,-10.,10.);
449         fOutList->Add(fHistVzMCrecoDiff);
450     }
451     
452     fHistVertexNconributors = new TH1I("fHistVertexNconributors","Primary vertex n contributors;N contributors;Entries",101,-0.5,100.5);
453     fOutList->Add(fHistVertexNconributors);
454     
455     fHistNumberOfPileupVerticesTracks = new TH1I("fHistNumberOfPileupVerticesTracks","Number of pilup verteces (by tracks);N verteces;Entries",11,-0.5,10.5);
456     fOutList->Add(fHistNumberOfPileupVerticesTracks);
457     
458     fHistNumberOfPileupVerticesSPD = new TH1I("fHistNumberOfPileupVerticesSPD","Number of pilup verteces (by SPD);N verteces;Entries",11,-0.5,10.5);
459     fOutList->Add(fHistNumberOfPileupVerticesSPD);
460     
461     
462     if ( fIsIonsAnalysis )
463     {
464         //Event plane
465         fHistEventPlane = new TH2F("fHistEventPlane",";#Psi, rad.;Centrality percentile;Counts",80,0,TMath::Pi()/*0,360.*/,101,-0.5,100.5);
466         fOutList->Add(fHistEventPlane);
467     }
468     //pt, eta, phi checkplots
469     fHistPt = new TH1F("fHistPt", "p_{T} distribution", 120, 0.0, 6.0);
470     fHistPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
471     fHistPt->GetYaxis()->SetTitle("dN/dp_{T} (c/GeV)");
472     fHistPt->SetMarkerStyle(kFullCircle);
473     fOutList->Add(fHistPt);
474     
475     fHistEta = new TH1F("fHistEta", "#eta distribution", 200, -4, 4);
476     fHistEta->GetXaxis()->SetTitle("#eta");
477     fHistEta->GetYaxis()->SetTitle("dN/#eta");
478     fHistEta->SetMarkerStyle(kFullCircle);
479     fOutList->Add(fHistEta);
480     
481     fHistEtaAODpure = new TH1F("fHistEtaAODpure", "#eta distribution AOD before cuts;#eta;dN/d#eta", 200, -4, 4);
482     fOutList->Add(fHistEtaAODpure);
483     
484     fHistPhi = new TH1F("fHistPhi", "#phi distribution", 200, 0, 2*TMath::Pi() );
485     fHistPhi->GetXaxis()->SetTitle("#phi");
486     fHistPhi->GetYaxis()->SetTitle("dN/#phi");
487     fHistPhi->SetMarkerStyle(kFullCircle);
488     fOutList->Add(fHistPhi);
489     
490     fHistEtaPhi = new TH2D("fHistEtaPhi","N tracks in (#eta, #phi);#eta;#phi",25,-0.8,0.8,25,0,2*TMath::Pi());
491     fOutList->Add(fHistEtaPhi);
492     
493     fHistPhiLRCrotationsCheck = new TH1F("fHistPhiLRCrotationsCheck", "#phi distribution in LRC rotations;#phi;dN/#phi", 200, -4*TMath::Pi(), 4*TMath::Pi() );
494     fOutList->Add(fHistPhiLRCrotationsCheck);
495     
496     //some histos for tracks manipulations by hand
497     fHistPhiArtificialProfilingCheck = new TH1F("fHistPhiArtificialProfilingCheck", "#phi distribution tracks distr by hand;#phi;dN/#phi", 200, -4*TMath::Pi(), 4*TMath::Pi() );
498     fHistPhiArtificialProfilingCheckWrtEvPlane = new TH1F("fHistPhiArtificialProfilingCheckWrtEvPlane", "#phi distribution by hand wrt event plane;#phi;dN/#phi", 200, -4*TMath::Pi(), 4*TMath::Pi() );
499     fHistPhiArtificialEvPlane = new TH1F("fHistPhiArtificialEvPlane", "#phi distribution of event plane;#phi;dN/#phi", 200, -4*TMath::Pi(), 4*TMath::Pi() );
500     if ( fUsePhiShufflingByHand )
501     {
502         fOutList->Add(fHistPhiArtificialProfilingCheck);
503         fOutList->Add(fHistPhiArtificialProfilingCheckWrtEvPlane);
504         fOutList->Add(fHistPhiArtificialEvPlane);
505     }
506     
507     //Eta vs Zv
508     fHistEtaVsZvCoverage = new TH2D("fHistEtaVsZvCoverage","TPC tracks #eta vs Zv;V_{z} (cm);#eta",100,-20,20,50,-2,2);
509     fHistEtaVsZvCoverageAccepted = new TH2D("fHistEtaVsZvCoverageAccepted","Accepted TPC tracks #eta vs Zv;V_{z} (cm);#eta",100,-20,20,50,-2,2);
510     fOutList->Add(fHistEtaVsZvCoverage);
511     fOutList->Add(fHistEtaVsZvCoverageAccepted);
512     
513     int lMaxAcceptedTracksInHist = fIsIonsAnalysis ? 4001 : 101;
514     fHistMultBeforeCuts = new TH1D("fHistMultBeforeCuts","N_{ch} - tracks;N_{ch} tracks;Entries"
515                                    ,lMaxAcceptedTracksInHist,-0.5,lMaxAcceptedTracksInHist-0.5);
516     fOutList->Add(fHistMultBeforeCuts);
517     
518     fHistAcceptedMult = new TH1D("fHistAcceptedMult","N_{ch} - accepted tracks;N_{ch} accepted;Entries"
519                                  ,lMaxAcceptedTracksInHist,-0.5,lMaxAcceptedTracksInHist-0.5);
520     fOutList->Add(fHistAcceptedMult);
521     
522     fHistAcceptedTracks = new TH1D("fHistAcceptedTracks","N_{ch} - accepted tracks for LRC;N_{ch} accepted;Entries"
523                                    ,lMaxAcceptedTracksInHist,-0.5,lMaxAcceptedTracksInHist-0.5);
524     fOutList->Add(fHistAcceptedTracks);
525     
526     fHistMultiplicityInEtaRegion = new TH1D("fHistMultiplicityInEtaRegion","N_{ch} in |#eta| region;N_{ch};Entries"
527                                             ,lMaxAcceptedTracksInHist,-0.5,lMaxAcceptedTracksInHist-0.5);
528     fOutList->Add(fHistMultiplicityInEtaRegion);
529     
530     fHistMultiplicityInEtaRegionAfterPtCuts = new TH1D("fHistMultiplicityInEtaRegionAfterPtCuts","N_{ch} in |#eta| region after p_{T} cuts;N_{ch};Entries"
531                                                        ,lMaxAcceptedTracksInHist,-0.5,lMaxAcceptedTracksInHist-0.5);
532     fOutList->Add(fHistMultiplicityInEtaRegionAfterPtCuts);
533     
534     if ( fAnalyseMCESD )
535     {
536         fHist2DMultiplicityMCESDInEtaRegion = new TH2D("fHist2DMultiplicityMCESDInEtaRegion","N_{ch} in |#eta| region;N_{ch} MC;N_{ch} ESD;Entries"
537                                                        ,101,-0.5,100.5,101,-0.5,100.5);
538         fOutList->Add(fHist2DMultiplicityMCESDInEtaRegion);
539     }
540     
541     
542     fHistAcceptedTracksAfterPtCuts = new TH1D("fHistAcceptedTracksAfterPtCuts","N_{ch} - accepted tracks for LRC after Pt cuts;N_{ch} accepted;Entries"
543                                               ,lMaxAcceptedTracksInHist,-0.5,lMaxAcceptedTracksInHist-0.5);
544     fOutList->Add(fHistAcceptedTracksAfterPtCuts);
545     
546     fHistAcceptedTPCtracks = new TH1D("fHistAcceptedTPCtracks","N_{ch} - accepted tracks with TPC inner param;N_{ch} accepted;Entries"
547                                       ,lMaxAcceptedTracksInHist,-0.5,lMaxAcceptedTracksInHist-0.5);
548     fOutList->Add(fHistAcceptedTPCtracks);
549     
550     fHistClustersTPC = new TH1D("fHistClustersTPC","N Clusters TPC;N_{TPC clusters};Entries",161,-0.5,160.5);
551     fOutList->Add(fHistClustersTPC);
552     
553     fHistClustersTPCafterCuts = new TH1D("fHistClustersTPCafterCuts","N Clusters TPC after cuts;N_TPC_clusters_after_cuts;Entries",161,-0.5,160.5);
554     fOutList->Add(fHistClustersTPCafterCuts);
555     
556     
557     fHistCrossedRowsTPC = new TH1D("fHistCrossedRowsTPC","N Crossed Rows TPC;N_{TPC CrossedRows};Entries",161,-0.5,160.5);
558     fOutList->Add(fHistCrossedRowsTPC);
559     
560     fHistCrossedRowsTPCafterCuts = new TH1D("fHistCrossedRowsTPCafterCuts","N CrossedRows TPC after cuts;N_TPC_clusters_after_cuts;Entries",161,-0.5,160.5);
561     fOutList->Add(fHistCrossedRowsTPCafterCuts);
562     
563     
564     
565     fHistTrackletsITS = new TH1D("fHistTrackletsITS","N Tracklets ITS;N_ITS_tracklets;Entries",101,-0.5,100.5);
566     fOutList->Add(fHistTrackletsITS);
567     
568     
569     fHistClustersITS = new TH1D("fHistClustersITS","N Clusters ITS;N_ITS_clusters;Entries",11,-0.5,10.5);
570     fOutList->Add(fHistClustersITS);
571     
572     
573     fHist2DClustersTPCvsPt = new TH2D("fHist2DClustersTPCvsPt","Num TPC clusters vs Pt;P_t;N_clusters",50,0,2,161,-0.5,160.5);
574     fOutList->Add(fHist2DClustersTPCvsPt);
575     
576     fHist2DClustersTPCvsEta = new TH2D("fHist2DClustersTPCvsEta","Num TPC clusters vs Eta;#eta;N_clusters",50,-1.5,1.5,161,-0.5,160.5);
577     fOutList->Add(fHist2DClustersTPCvsEta);
578     
579     fHist2DAcceptedTracksPtvsEta = new TH2D("fHist2DAcceptedTracksPtvsEta","Accepted Tracks Pt vs Eta;pt;#eta",50,0,2,50,-1.5,1.5);
580     fOutList->Add(fHist2DAcceptedTracksPtvsEta);
581     
582     fHistProbabilitiesPID = new TH1D("fHistProbabilitiesPID","PID Probabilities;pid;Entries",10,-0.5,9.5 );//AliPID::kSPECIES,-0.5,AliPID::kSPECIES-0.5);
583     fOutList->Add(fHistProbabilitiesPID);
584     
585     fHistESDtrackMass = new TH1D("fHistESDtrackMass","ESD Mass;Mass;Entries",100,0,1);
586     fOutList->Add(fHistESDtrackMass);
587     
588     fHistProbabilityPion = new TH1D("fHistProbabilityPion","Probability Pion;Probability;Entries",100,0,1);
589     fOutList->Add(fHistProbabilityPion);
590     
591     fHistProbabilityKaon = new TH1D("fHistProbabilityKaon","Probability Kaon;Probability;Entries",100,0,1);
592     fOutList->Add(fHistProbabilityKaon);
593     
594     fHistProbabilityProton = new TH1D("fHistProbabilityProton","Probability Proton;Probability;Entries",100,0,1);
595     fOutList->Add(fHistProbabilityProton);
596     
597     
598     
599     fHistParticlesDistrBeforeCuts = new TH1D("fHistParticlesDistrBeforeCuts","Particles Distr;Particle;Entries",5,-0.5,4.5);
600     TString gBinParticleNames[5] = {/*"Other",*/"Electron","Muon","Pion","Kaon", "Proton"};
601     for(Int_t i = 1; i <= 5; i++)fHistParticlesDistrBeforeCuts->GetXaxis()->SetBinLabel(i,gBinParticleNames[i-1].Data());
602     fOutList->Add(fHistParticlesDistrBeforeCuts);
603     
604     fHistParticlesDistr = new TH1D("fHistParticlesDistr","Particles Distr;Particle;Entries",5,-0.5,4.5);
605     //TString gBinParticleNames[6] = {"Undefined","Electron","Muon","Pion","Kaon", "Proton"};
606     for(Int_t i = 1; i <= 5; i++)fHistParticlesDistr->GetXaxis()->SetBinLabel(i,gBinParticleNames[i-1].Data());
607     fOutList->Add(fHistParticlesDistr);
608     
609     fHistCentralityPercentile = new TH1D("fHistCentralityPercentile","Centrality Percentile;Centrality;Entries",101,-1.01,100.01);
610     if ( fIsIonsAnalysis ) fOutList->Add(fHistCentralityPercentile);
611     
612     fHistCentralityClass10 = new TH1D("fHistCentralityClass10","Centrality Class 10;Centrality;Entries",10,-0.5,9.5);
613     if ( fIsIonsAnalysis ) fOutList->Add(fHistCentralityClass10);
614     
615     fHistCentralityClass5 = new TH1D("fHistCentralityClass5","Centrality Class 5;Centrality;Entries",20,-0.5,19.5);
616     if ( fIsIonsAnalysis ) fOutList->Add(fHistCentralityClass5);
617     
618     
619     fRand = new TRandom3;
620     //ZDC additional study
621     
622     fMultForZDCstudy[0] = 0;
623     fMultForZDCstudy[1] = 200;
624     fMultForZDCstudy[2] = 220;
625     fMultForZDCstudy[3] = 240;
626     fMultForZDCstudy[4] = 250;
627     
628     TString strZDCenergyName;
629     TString strZDCenergyTitle;
630     for ( int i = 0; i < 5; i++ )
631     {
632         strZDCenergyName = Form( "fHistZDCenergy%d", fMultForZDCstudy[i] );
633         strZDCenergyTitle = Form(  "ZDC Energy, multThreshold %d;Energy;Entries", fMultForZDCstudy[i] );
634         fHistZDCenergy[i] = new TH1D(strZDCenergyName, strZDCenergyTitle, 200, 0, 10000);
635         if ( fFlagWatchZDC && fIsIonsAnalysis )
636             fOutList->Add(fHistZDCenergy[i]);
637     }
638     
639     
640     fHistZDCparticipants = new TH1D("fHistZDCparticipants","ZDC Participants;Participants;Entries",800,0,800);
641     if ( fFlagWatchZDC && fIsIonsAnalysis )
642         fOutList->Add(fHistZDCparticipants);
643     
644     if ( fFlagWatchV0 )//fIsIonsAnalysis )
645     {
646         const int nBinsV0multForNonIons = 700;
647         const int nBinsV0multForIons = 220;
648         const int nMaxV0multForIons = 22000;
649         fHistV0multiplicity = new TH1D("fHistV0multiplicity","V0 Multiplicity;Multiplicity;Entries",
650                                        fIsIonsAnalysis ? nBinsV0multForIons : nBinsV0multForNonIons, 0, fIsIonsAnalysis ? nMaxV0multForIons : nBinsV0multForNonIons);
651         fHistV0Amultiplicity = new TH1D("fHistV0Amultiplicity","V0-A Multiplicity;Multiplicity;Entries",
652                                         fIsIonsAnalysis ? nBinsV0multForIons : nBinsV0multForNonIons, 0, fIsIonsAnalysis ? nMaxV0multForIons : nBinsV0multForNonIons);
653         fHistV0Cmultiplicity = new TH1D("fHistV0Cmultiplicity","V0-C Multiplicity;Multiplicity;Entries",
654                                         fIsIonsAnalysis ? nBinsV0multForIons : nBinsV0multForNonIons, 0, fIsIonsAnalysis ? nMaxV0multForIons : nBinsV0multForNonIons);
655         fHist2DV0ACmultiplicity = new TH2D("fHist2DV0ACmultiplicity","V0 A-C Multiplicity;Multiplicity A;Multiplicity C;Entries"
656                                            ,50, 0, fIsIonsAnalysis ? nMaxV0multForIons/2 : nBinsV0multForNonIons/2
657                                                                      ,50, 0, fIsIonsAnalysis ? nMaxV0multForIons/2 : nBinsV0multForNonIons/2);
658         
659         fHist2DTracksAcceptedVsV0multiplicity = new TH2D("fHist2DTracksAcceptedVsV0multiplicity","V0 Multiplicity vs Accepted;N Accepted tracks;Multiplicity V0;Entries"
660                                                          ,70, 0, fIsIonsAnalysis ? 3000 : 70
661                                                                                    ,50, 0, fIsIonsAnalysis ? nMaxV0multForIons/2 : 100 );
662         
663         //number of fired cells
664         const int nV0cells = 32;
665         //fHistV0cells = new TH1D("fHistV0cells","V0 cells;N cells;Entries", nV0cells+1, -0.5, nV0cells + 0.5 );
666         fHistV0Acells = new TH1D("fHistV0Acells","V0-A cells;N cells;Entries", nV0cells+1, -0.5, nV0cells + 0.5 );
667         fHistV0Ccells = new TH1D("fHistV0Ccells","V0-C cells;N cells;Entries", nV0cells+1, -0.5, nV0cells + 0.5 );
668         fHist2DV0ACcells = new TH2D("fHist2DV0ACcells","V0 A-C cells;N cells A;N cells C;Entries", nV0cells+1, -0.5, nV0cells + 0.5, nV0cells+1, -0.5, nV0cells + 0.5 );
669         
670         
671         fOutList->Add(fHistV0multiplicity);
672         fOutList->Add(fHistV0Amultiplicity);
673         fOutList->Add(fHistV0Cmultiplicity);
674         fOutList->Add(fHist2DV0ACmultiplicity);
675         fOutList->Add(fHist2DTracksAcceptedVsV0multiplicity);
676         //cells
677         //fOutList->Add(fHistV0cells    );
678         fOutList->Add(fHistV0Acells   );
679         fOutList->Add(fHistV0Ccells   );
680         fOutList->Add(fHist2DV0ACcells);
681         
682         //mult in rings
683         //        char strV0ringName[200];
684         //        char strV0ringTitle[200];
685         TString strV0ringName;
686         TString strV0ringTitle;
687         for ( int i = 0; i < 4; i++ )
688         {
689             strV0ringName =  Form( "fHistV0AmultiplicityRing%d", i );
690             strV0ringTitle = Form( strV0ringTitle, "V0-A Multiplicity Ring %d;Multiplicity;Entries", i );
691             fHistV0AmultiplicityRing[i] = new TH1D( strV0ringName, strV0ringTitle,
692                                                     fIsIonsAnalysis ? nBinsV0multForIons : nBinsV0multForNonIons, 0, fIsIonsAnalysis ? nMaxV0multForIons : nBinsV0multForNonIons);
693             strV0ringName =  Form(  "fHistV0CmultiplicityRing%d", i );
694             strV0ringTitle = Form( strV0ringTitle, "V0-C Multiplicity Ring %d;Multiplicity;Entries", i );
695             fHistV0CmultiplicityRing[i] = new TH1D( strV0ringName, strV0ringTitle,
696                                                     fIsIonsAnalysis ? nBinsV0multForIons : nBinsV0multForNonIons, 0, fIsIonsAnalysis ? nMaxV0multForIons : nBinsV0multForNonIons);
697             strV0ringName =  Form( "fHist2DV0ACmultiplicityRing%d", i );
698             strV0ringTitle = Form(  "V0-AC Multiplicity Ring %d;Multiplicity A;Multiplicity C;Entries", i );
699             fHist2DV0ACmultiplicityRing[i] = new TH2D( strV0ringName, strV0ringTitle
700                                                        , 100, 0, fIsIonsAnalysis ? nMaxV0multForIons/2 : 100
701                                                                                    , 100, 0, fIsIonsAnalysis ? nMaxV0multForIons/2 : 100 );
702             fOutList->Add(fHistV0AmultiplicityRing[i]);
703             fOutList->Add(fHistV0CmultiplicityRing[i]);
704             fOutList->Add(fHist2DV0ACmultiplicityRing[i]);
705             
706             //mult in barrel vs V0 rings
707             strV0ringName =  Form( "fHist2DTracksAcceptedVsV0AmultiplicityRing%d", i );
708             strV0ringTitle = Form(  "Accepted tracks vs V0-A Multiplicity in Ring %d;N Accepted tracks;Multiplicity V0A;Entries", i );
709             fHist2DTracksAcceptedVsV0AmultiplicityRing[i] = new TH2D( strV0ringName, strV0ringTitle
710                                                                       , 100, 0, fIsIonsAnalysis ? nMaxV0multForIons/2 : 100
711                                                                                                   , 100, 0, fIsIonsAnalysis ? nMaxV0multForIons/2 : 100 );
712             strV0ringName =  Form( "fHist2DTracksAcceptedVsV0CmultiplicityRing%d", i );
713             strV0ringTitle = Form(  "Accepted tracks vs V0-C Multiplicity in Ring %d;N Accepted tracks;Multiplicity V0C;Entries", i );
714             fHist2DTracksAcceptedVsV0CmultiplicityRing[i] = new TH2D( strV0ringName, strV0ringTitle
715                                                                       , 100, 0, fIsIonsAnalysis ? nMaxV0multForIons/2 : 100
716                                                                                                   , 100, 0, fIsIonsAnalysis ? nMaxV0multForIons/2 : 100 );
717             
718             fOutList->Add(fHist2DTracksAcceptedVsV0AmultiplicityRing[i]);
719             fOutList->Add(fHist2DTracksAcceptedVsV0CmultiplicityRing[i]);
720         }
721         
722     }
723     
724     
725     
726     
727     //    fHistV0spectra = new TH1D("fHistV0spectra","V0 spectra;Mass, GeV;Entries",500,0,5000);
728     //    fOutList->Add(fHistV0spectra);
729     
730     
731     fHistMClabels = new TH1D("fHistMClabels","MC label;label;Entries",102,-100.5,100.5);
732     fOutList->Add(fHistMClabels);
733     
734     fHistRejectedTracksCharge = new TH1D("fHistRejectedTracksCharge","Rejected tracks charge;charge;Entries",3,-1.5,1.5);
735     fOutList->Add(fHistRejectedTracksCharge);
736     
737     fHistTracksCharge = new TH1D("fHistTracksCharge","Accepted tracks charge;charge;Entries",3,-1.5,1.5);
738     fOutList->Add(fHistTracksCharge);
739
740     // ##### net charge vs pt study
741     const int kPtNetChargePtBins = 200;
742     const double kPtNetChargePtMin = 0.1;
743     const double kPtNetChargePtMax = 2.1;
744
745     //pt plus
746     fHistPtPlus = new TH1D("fHistPtPlus","p_{T} +;p_{T};dN/dpT"
747                                   , kPtNetChargePtBins,  kPtNetChargePtMin, kPtNetChargePtMax
748                                   );
749     fOutList->Add(fHistPtPlus);
750
751     //pt minus
752     fHistPtMinus = new TH1D("fHistPtMinus","p_{T} -;p_{T};dN/dpT"
753                                   , kPtNetChargePtBins,  kPtNetChargePtMin, kPtNetChargePtMax
754                                   );
755     fOutList->Add(fHistPtMinus);
756
757     //net charge vs pT
758     fHistNetChargeVsPt = new TH1D("fHistNetChargeVsPt","charge vs p_{T};p_{T};Q"
759                                   , kPtNetChargePtBins,  kPtNetChargePtMin, kPtNetChargePtMax
760                                   );
761     fOutList->Add(fHistNetChargeVsPt);
762
763     //plus
764     fHistChargePlusVsPtTmp = new TH1D( "fHistChargePlusVsPtTmp","charge plus vs p_{T};p_{T};q plus"
765                                        , kPtNetChargePtBins,  kPtNetChargePtMin, kPtNetChargePtMax
766                                        );
767     //    fOutList->Add(fHistChargePlusVsPtTmp);
768
769     //minus
770     fHistChargeMinusVsPtTmp = new TH1D( "fHistChargeMinusVsPtTmp","charge minus vs p_{T};p_{T};q minus"
771                                         , kPtNetChargePtBins,  kPtNetChargePtMin, kPtNetChargePtMax
772                                         );
773     //    fOutList->Add(fHistChargeMinusVsPtTmp);
774
775     //NetChargeVsPt
776     fHist2DNetChargeVsPt = new TH2D( "fHist2DNetChargeVsPt","Net charge vs p_{T};p_{T};Q"
777                                      , kPtNetChargePtBins,  kPtNetChargePtMin, kPtNetChargePtMax
778                                      , 40,  -20, 20
779                                      );
780     fOutList->Add(fHist2DNetChargeVsPt);
781
782     //NetChargeVsPt CorrectedOnEventMean
783     fHist2DNetChargeVsPtCorrectedOnEventMean = new TH2D( "fHist2DNetChargeVsPtCorrectedOnEventMean","Net charge corrected on mean e-by-e vs p_{T};p_{T};Q"
784                                                          , kPtNetChargePtBins,  kPtNetChargePtMin, kPtNetChargePtMax
785                                                          , 40,  -20, 20
786                                                          );
787     fOutList->Add(fHist2DNetChargeVsPtCorrectedOnEventMean);
788
789     //NetChargeVsPt CorrectedOnEventMean normalized on Nch
790     fHist2DNetChargeVsPtCorrectedOnEventMeanNormOnNch = new TH2D( "fHist2DNetChargeVsPtCorrectedOnEventMeanNormOnNch","Net charge vs p_{T} corrected on mean e-by-e normalized on nCh;p_{T};Q"
791                                                                   , kPtNetChargePtBins,  kPtNetChargePtMin, kPtNetChargePtMax
792                                                                   , 40,  -20, 20
793                                                                   );
794     fOutList->Add(fHist2DNetChargeVsPtCorrectedOnEventMeanNormOnNch);
795
796
797
798     if ( fArtificialInefficiency >= 0 ) //i.e. have this kind of analysis
799     {
800         fHistNumberOfDroppedByHandTracks = new TH2D("fHistNumberOfDroppedByHandTracks","Accepted tracks vs Dropped artificially;N_{ch} accepted;N_{ch} dropped;Entries", 71, -0.5, 70.5,   71, -0.5, 70.5);
801         fOutList->Add(fHistNumberOfDroppedByHandTracks);
802     }
803     
804     
805     // Returning TH1,TH2 AddDirectory status
806     TH1::AddDirectory(lTH1oldStatus);
807     TH2::AddDirectory(lTH2oldStatus);
808     
809     //fHistProbabilitiesPID->Fill(1,fStrPIDforFwd);
810     //fHistProbabilitiesPID->Fill(2,fStrPIDforBwd);
811     Printf("UserCreateOutputObjects done.");
812     
813     // NEW HISTO added to fOutput here
814     PostData(1, fOutList); // Post data for ALL output slots >0 here, to get at least an empty histogram
815     
816     for( Int_t i = 0; i < fLRCproc.GetEntries(); i++)
817     {
818         PostData( Proc(i)->GetOutputSlotNumber(),Proc(i)->CreateOutput() );
819     }
820     
821     
822     //if ( fSetIncludeEventTreeInOutput )
823     //    PostData(2, fEventTree);
824     //int a;
825     //cin >> a;
826     
827     
828 }
829
830 //________________________________________________________________________
831 //void AliAnalysisTaskLRC::UserExec(Option_t *)
832 //{
833
834 //    //cout << "TEST in Task" << endl;
835 //    // ###### tuning phi
836
837 //    //rotate phi when N_phi_sectors > 1 and call UserExecLoop:
838 //    double phiStep = 2 * TMath::Pi();
839 //    if ( fNumberOfSectors > 1 )
840 //        phiStep /= fNumberOfSectors;
841
842 //    double lPhiRotatedExtra = 0; //additional phi rotation of all tracks
843 //    for ( Int_t sectorId = 0; sectorId < fNumberOfSectors; sectorId++ )
844 //    {
845 //        //cout << "loop " << sectorId << endl;
846 //        UserExecLoop( lPhiRotatedExtra );
847 //        lPhiRotatedExtra += phiStep; //rotate track
848 //    }
849
850
851 //}
852
853 //________________________________________________________________________
854 void AliAnalysisTaskLRC::UserExec(Option_t *)   //UserExecLoop( Double_t phiAdditional )//Option_t *)
855 {
856     // ########### if use toy events
857     if ( fUseToyEvents )
858     {
859         //generate all events here and fill LRC processors
860         UseToyEvents();
861
862         //just return
863         return;
864     }
865     
866     // Main loop
867     // Called for each event
868     //printf( "starting UserExec...\n" );
869     //Pid setting
870     if ( fPIDsensingFlag )
871     {
872         SetParticleTypeToProcessors( 0 /*fwd*/, fStrPIDforFwd );
873         SetParticleTypeToProcessors( 1 /*backward*/, fStrPIDforBwd );
874         fPIDsensingFlag = kFALSE;
875     }
876     
877     //printf( "starting event...\n" );
878     //The common PID object can then be retrieved from the input handler:
879     AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
880     AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
881     fPIDResponse = inputHandler->GetPIDResponse();
882     //if (!fPIDResponse) !!!!!!!!!!!!
883     //  AliFatal("This Task needs the PID response attached to the inputHandler");
884     
885     
886     AliVEvent *event = InputEvent();
887     AliESDEvent *lESD = 0x0;
888     AliAODEvent *lAOD = 0x0;
889     
890     AliStack *stack = 0x0;
891     AliMCEvent *eventMC = 0x0;
892     
893     
894     if ( fAnalysisLevel == "ESD" ) //fAnalysisType == en_AnalysisType_AOD )
895     {
896         lESD = dynamic_cast<AliESDEvent*> (event) ;
897         //printf( "cast to ESD is Ok\n." );
898     }
899     else if ( fAnalysisLevel == "AOD" ) //fAnalysisType == en_AnalysisType_AOD )
900     {
901         //cout << "test AOD analysis... " << endl;
902         lAOD = dynamic_cast<AliAODEvent*>(event); // from TaskSE
903         //Int_t bMagneticFieldSign = (lAOD->GetMagneticField() > 0) ? 1 : -1;
904         //printf( "Number of AOD tracks = %d, magnetic field = %f\n", lAOD->GetNumberOfTracks(), lAOD->GetMagneticField() );
905         //return;
906         
907     }
908     
909     if( fRunKine ) //|| fAnalyseMCESD )
910     {
911         eventMC = MCEvent();
912         stack = eventMC->Stack();
913         //Printf("Number of primaries: %d",stack->GetNprimary());
914     }
915     
916     if (!event)
917     {
918         //Printf("ERROR: Could not retrieve event");
919         return;
920     }
921     //cout << "test start Event!" << endl;
922     
923     if( (lESD ) && ( !fEsdTrackCuts ) )
924     {
925         AliDebug(AliLog::kError, "No ESD track cuts avalible");
926         return;
927     }
928     
929     //make combined flag for esd selection during MC kine analysis
930     Bool_t lRunKine = ( fRunKine && !fAnalyseMCESD );
931     
932     // Processing event selection
933     fHistEventCutStats->Fill( "Total", 1 );
934     
935     
936     
937     //cout << "test: nOfTracks = " << event->GetNumberOfTracks() << endl;
938     
939     //Trigger
940     Bool_t lTrigger = kTRUE;
941     if( lESD && !lRunKine ) //just for tests; usually it is done before UserExecLoop in PhysSel task
942     {
943         if(fShowEventStats)
944             Printf("Trigger classes: %s:", lESD->GetFiredTriggerClasses().Data());
945         
946         lTrigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
947         
948         //8 TeV runs: TEST triggers! (29.12.2012)
949         //        Bool_t lTrigger1 = (lESD->IsTriggerClassFired("CINT7-B-NOPF-ALLNOTRD")) ? 1 : 0;
950         //        Bool_t lTrigger2 = (lESD->IsTriggerClassFired("CINT7WU-B-NOPF-ALL")) ? 1 : 0;
951         //        lTrigger = lTrigger1 && lTrigger2;
952         
953         if( !lTrigger )
954         {
955             if( fShowEventStats )
956                 Printf("Rejected!");
957             fHistEventCutStats->Fill( "No trigger", 1 );
958             PostData(1, fOutList);
959             return;
960         }
961     }
962     //Centrality      //https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies
963     Double_t lCentralityPercentile      = 0.;
964     Int_t lCentralityClass10            = 0;
965     Int_t lCentralityClass5             = 0;
966     Float_t lReactionPlane       = -1.;
967     
968     
969     if( !lRunKine && fIsIonsAnalysis )
970     {
971         AliCentrality *centrality = 0x0;
972         
973         if ( lESD )
974             centrality = lESD->GetCentrality();
975         else if ( lAOD )
976             centrality = lAOD->GetCentrality();
977         
978         if ( centrality )
979         {
980             lCentralityPercentile       = centrality->GetCentralityPercentile("V0M");   // returns the centrality percentile, a float from 0 to 100 (or to the trigger efficiency)
981             lCentralityClass10          = centrality->GetCentralityClass10("V0M");//"FMD");             // returns centrality class for 10% width (a integer from 0 to 10)
982             lCentralityClass5           = centrality->GetCentralityClass5("V0M");//"TKL");              // returns centrality class for 5% width (a integer from 0 to 20)
983             Bool_t lCentralityInClass = centrality->IsEventInCentralityClass(fMinCentralityClass,fMaxCentralityClass,"V0M"); // returns kTRUE if the centrality of the event is between a and b, otherwise kFALSE
984             //cout << lCentralityPercentile << " "
985             //<< fMinCentralityClass << " "
986             //<< fMaxCentralityClass << " "
987             //<< lCentralityInClass << endl;
988             
989             if ( !lCentralityInClass )
990             {
991                 //cout << "outside of centrality class!" << endl;
992                 fHistEventCutStats->Fill("Wrong centrality", 1);
993                 PostData(1, fOutList);
994                 return;
995             }
996             
997             fHistCentralityPercentile->Fill( lCentralityPercentile );
998             fHistCentralityClass10->Fill( lCentralityClass10 );
999             fHistCentralityClass5->Fill( lCentralityClass5 );
1000             
1001             // get the reaction plane
1002             lReactionPlane = GetEventPlane( event );
1003             fHistEventPlane->Fill( lReactionPlane, lCentralityPercentile );
1004         }
1005     }
1006     
1007     //number of verteces in ESD (pileup)
1008     if ( lESD )
1009     {
1010         int lNumberOfPileUpVerteces = 0;
1011         TClonesArray *lPileupVertecesTracks = lESD->GetPileupVerticesTracks();
1012         lNumberOfPileUpVerteces = lPileupVertecesTracks->GetSize();
1013         fHistNumberOfPileupVerticesTracks->Fill( lNumberOfPileUpVerteces );
1014         
1015         TClonesArray *lPileupVertecesSPD = lESD->GetPileupVerticesSPD();
1016         lNumberOfPileUpVerteces = lPileupVertecesSPD->GetSize();
1017         fHistNumberOfPileupVerticesSPD->Fill( lNumberOfPileUpVerteces );
1018     }
1019     else if ( lAOD ) //number of verteces in AOD (pileup)
1020     {
1021         int lNumberOfPileUpVerteces = 0;
1022         lNumberOfPileUpVerteces = lAOD->GetNumberOfPileupVerticesTracks();
1023         fHistNumberOfPileupVerticesTracks->Fill( lNumberOfPileUpVerteces );
1024         lNumberOfPileUpVerteces = lAOD->GetNumberOfPileupVerticesSPD();
1025         fHistNumberOfPileupVerticesSPD->Fill( lNumberOfPileUpVerteces );
1026     }
1027     
1028     
1029     // Vertex present
1030     //    const AliESDVertex *vertex = lESD->GetPrimaryVertex();
1031     Double_t lVertexX = -1000;
1032     Double_t lVertexY = -1000;
1033     Double_t lVertexZ = -1000;
1034     Double_t lMCVertexX = -1000;
1035     Double_t lMCVertexY = -1000;
1036     Double_t lMCVertexZ = -1000;
1037     
1038     Bool_t lReconstructedVertexPresent = kFALSE;
1039     Bool_t lVertexAcceptable = kFALSE;
1040     
1041     if ( fAnalysisLevel == "ESD" ||  fAnalysisLevel == "AOD" )
1042     {
1043         const AliVVertex *vertex = event->GetPrimaryVertex();
1044         fHistVertexNconributors->Fill( vertex->GetNContributors() );
1045         
1046         Double32_t lCov[6];
1047         vertex->GetCovarianceMatrix(lCov);
1048         
1049         //check nContributors and z*z>0
1050         lReconstructedVertexPresent = ( ( vertex->GetNContributors() > 0 ) && ( lCov[5] != 0 ) );
1051         //        if ( (vertex) && (vertex->GetNContributors() > 0 ) && ( vertex->GetZRes() == 0 ) )
1052         //        {
1053         //            int aa;
1054         //            cout << "stop with nContributors = " << vertex->GetNContributors() << endl;
1055         //            cin >> aa;
1056         //        }
1057         lVertexX = vertex->GetX();
1058         lVertexY = vertex->GetY();
1059         lVertexZ = vertex->GetZ();
1060     }
1061     
1062     if ( fRunKine ) //take MC vertex params and rewrite position
1063     {
1064         AliGenEventHeader *header = eventMC->GenEventHeader();
1065         if(!header) return;
1066         
1067         TArrayF gVertexArray;
1068         header->PrimaryVertex(gVertexArray);
1069         lMCVertexX = gVertexArray.At(0);
1070         lMCVertexY = gVertexArray.At(1);
1071         lMCVertexZ = gVertexArray.At(2);
1072         
1073         //fill mc-reco diff
1074         if ( lReconstructedVertexPresent )
1075         {
1076             fHistVxMCrecoDiff->Fill( lMCVertexX - lVertexX );
1077             fHistVyMCrecoDiff->Fill( lMCVertexY - lVertexY );
1078             fHistVzMCrecoDiff->Fill( lMCVertexZ - lVertexZ );
1079         }
1080         
1081     }
1082     
1083     //cut on reco vertex params
1084     if ( !lRunKine )
1085     {
1086         if( ( !lReconstructedVertexPresent ) && fCheckForkVtx )
1087         {
1088             if( fShowEventStats )
1089                 Printf("Bad vertex params");
1090             fHistEventCutStats->Fill("Bad vertex params",1);
1091             PostData(1, fOutList);
1092             return;
1093         }
1094     }
1095     else
1096     {
1097         //rewrite vert x,y,z into mc values
1098         lVertexX = lMCVertexX;
1099         lVertexY = lMCVertexY;
1100         lVertexZ = lMCVertexZ;
1101     }
1102     
1103     
1104     // Vertex in range
1105     lVertexAcceptable = (TMath::Abs(lVertexX) < fVxMax) && (TMath::Abs(lVertexY) < fVyMax);
1106     if( lVertexAcceptable )
1107     {
1108         if( fVzMax > 0 )   //   fVzMax < 0 -> select Zv outside selected range
1109             lVertexAcceptable = (TMath::Abs(lVertexZ) < fVzMax);
1110         else
1111             lVertexAcceptable = (TMath::Abs(lVertexZ) > -fVzMax);
1112     }
1113     
1114     if( (!lVertexAcceptable) && fCheckForVtxPosition && !lRunKine )
1115     {
1116         if(fShowEventStats)
1117             Printf("Vertex out of range");
1118         fHistEventCutStats->Fill("Bad vertex position",1);
1119         PostData(1, fOutList);
1120         return;
1121     }
1122     
1123     //fill bin if no reco vertex in MC analysis
1124     if ( lRunKine && !lReconstructedVertexPresent )
1125         fHistEventCutStats->Fill( "No reco vertex", 1 );
1126     
1127     fHistVx->Fill(lVertexX);
1128     fHistVy->Fill(lVertexY);
1129     fHistVz->Fill(lVertexZ);
1130     //end with vertex
1131     
1132     
1133     
1134     
1135     //cut on number of SPD tracklets (25.03.2012)
1136     if( lESD && !lRunKine )
1137     {
1138         //How to get tracklets
1139         const AliMultiplicity *tracklets = lESD->GetMultiplicity();
1140         Int_t multSPD = tracklets->GetNumberOfTracklets();
1141         //Int_t nITStracklets = kalmanTrack->GetNumberOfTracklets();
1142         fHistTrackletsITS->Fill( multSPD );
1143         if ( multSPD < fMinNumberOfSPDtracklets )
1144         {
1145             if(fShowEventStats)
1146                 Printf("Too few SPD tracklets");
1147             fHistEventCutStats->Fill( "Few SPD tracklets", 1 );
1148             PostData(1, fOutList);
1149             return;
1150         }
1151     }
1152     
1153     
1154     //tmp eta ranges (for ZDC respond study)
1155     double etaBmin = -0.8;
1156     double etaBmax = -0.6;
1157     double etaFmin = 0.6;
1158     double etaFmax = 0.8;
1159     int countTracksEtaB = 0;
1160     int countTracksEtaF = 0;
1161     
1162     //n accepted tracks
1163     int lNchTrigger = 0;
1164     int lTPCtracks  = 0; // added 25.03.2012
1165     fHistMultBeforeCuts->Fill( event->GetNumberOfTracks() );
1166     
1167     // Pre event loop
1168     if( !lRunKine )
1169     {
1170         for ( Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++)
1171         {
1172             if ( fAnalysisLevel == "ESD" )
1173             {
1174                 AliESDtrack *track = lESD->GetTrack(iTracks);
1175                 if ( fEsdTrackCuts->AcceptTrack(track) )
1176                 {
1177                     lNchTrigger++;
1178                     //todo: try to implement pt cuts here?....
1179                     //How to get TPC standalone tracks
1180                     AliExternalTrackParam *tpcTrack
1181                             = (AliExternalTrackParam *)track->GetTPCInnerParam();
1182                     if ( tpcTrack )
1183                         lTPCtracks++;
1184                 }
1185                 else
1186                 {
1187                     fHistRejectedTracksCharge->Fill( track->Charge() );
1188                 }
1189             }
1190             else if ( fAnalysisLevel == "AOD" )
1191             {
1192                 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1193                 if( aodTrack->TestFilterBit( fAODtrackCutBit ) )
1194                     lNchTrigger++;
1195             }
1196         }
1197     }
1198     fHistAcceptedTPCtracks->Fill( lTPCtracks );
1199     
1200     
1201     
1202     
1203     //if(lRunKine)lNchTrigger=eventMC->GetNumberOfPrimaries();
1204     
1205     // Nch bins cut
1206     if( ( lNchTrigger > fMaxAcceptedTracksCut ) && ( fMaxAcceptedTracksCut != 0 ) )
1207     {
1208         fHistEventCutStats->Fill( "HighMult cut", 1 );
1209         PostData(1, fOutList);
1210         return;
1211     }
1212     
1213     if( lNchTrigger < fMinAcceptedTracksCut )
1214     {
1215         fHistEventCutStats->Fill( "LowMult cut", 1 );
1216         PostData(1, fOutList);
1217         return;
1218     }
1219     fHistAcceptedMult->Fill(lNchTrigger);
1220     fHistEventCutStats->Fill("Analyzed",1);
1221     
1222     if ( lESD && fFlagWatchV0 ) // cut on V0 multiplicity "radius" in 2D-hist for both A and C sides
1223     {
1224         const AliESDVZERO* vzrData = lESD->GetVZEROData(); //aod the same
1225         
1226         //        const int lThresholdMultV0RingId = 3; //which ring is considered
1227         //        Double_t lThisEventV0MultSumRing3 = vzrData->GetMRingV0A(2) + vzrData->GetMRingV0C(2);
1228         //        Double_t lThisEventV0MultSumRing4 = vzrData->GetMRingV0A(3) + vzrData->GetMRingV0C(3);
1229         //        Double_t lThisEventV0MultSum = lThisEventV0MultSumRing3 + lThisEventV0MultSumRing4;
1230         double sumV0Amult = vzrData->GetMTotV0A();
1231         double sumV0Cmult = vzrData->GetMTotV0C();
1232
1233         Double_t lThisEventV0MultSum = sumV0Amult + sumV0Cmult;
1234
1235         //sqrt ( pow ( vzrData->GetMRingV0A(lThresholdMultV0RingId), 2 ) + pow ( vzrData->GetMRingV0C(lThresholdMultV0RingId), 2 ) );
1236         if ( lThisEventV0MultSum < fThresholdOnV0mult
1237              || sumV0Amult < fThresholdOnV0mult / 5 //additionally, to "exclude" diffraction
1238              || sumV0Cmult < fThresholdOnV0mult / 5 //additionally, to "exclude" diffraction
1239              )
1240         {
1241             PostData(1, fOutList);
1242             return;
1243         }
1244     }
1245     
1246     //    if ( fAnalysisLevel == "AOD" )
1247     //        return;
1248     
1249
1250     //reset tmp histos
1251     fHistChargePlusVsPtTmp->Reset();
1252     fHistChargeMinusVsPtTmp->Reset();
1253
1254     //Track selection counters
1255     int lNaccept=0;
1256     int lNacceptEtaInRegion=0;
1257     int lNacceptAfterPtCuts=0;
1258     int lNacceptEtaInRegionAfterPtCuts=0;
1259     //    int lNacceptPlusEtaInRegion=0;
1260     //    int lNacceptMinusEtaInRegion=0;
1261     int lPtOver=0;
1262     int lPtUnder=0;
1263     int lNoCharge=0;
1264     int lCutfail=0;
1265     int lDecay=0;
1266     //AliLRCBase::LRCparticleType lParticleType = AliLRCBase::kLRCany;
1267     
1268     int lNumberOfDropByHandTracks = 0; //number of artificially dropped tracks
1269     
1270     //Double_t probTPC[AliPID::kSPECIES]={0.};
1271     Double_t probTOF[AliPID::kSPECIES]={0.};
1272     Double_t probTPCTOF[AliPID::kSPECIES]={0.};
1273     
1274     Int_t nTracksInEvent = ( !lRunKine )
1275             ? event->GetNumberOfTracks()
1276             : eventMC->GetNumberOfPrimaries();//GetNumberOfTracks();
1277     
1278     //        if ( fSetIncludeEventTreeInOutput )
1279     //        {
1280     //            fSimpleEvent->SetHeader( fNsimpleEvents, -1, -1, lCentralityPercentile, lReactionPlane );
1281     //            fSimpleEvent->SetVertexPos( vertex->GetX(), vertex->GetY(), vertex->GetZ() );
1282     //            fNsimpleEvents++;
1283     //            fSimpleEvent->StartEventFilling();
1284     //        }
1285     //if ( !lRunKine )
1286     int lNumberOfAcceptedTracksForLRC = 0;
1287     // Track loop -----------------------------------------------------------------
1288     for (Int_t iTracks = 0; iTracks < nTracksInEvent/*event->GetNumberOfTracks()*/; iTracks++)
1289     {
1290         
1291         //Track variables
1292         double lPt = 0;   // Temp Pt
1293         double lEta = -100;       // Temp ETA
1294         double lPhi = -100;    // Temp Phi
1295         double lMass = 0;    // Temp Mass
1296         double lEt = 0;
1297         Short_t lCharge = 0;
1298         
1299         AliVParticle* track = ( !lRunKine )
1300                 ? event->GetTrack(iTracks)
1301                 : eventMC->GetTrack(iTracks);
1302         
1303         
1304         if ( !track ) {
1305             Printf("ERROR: Could not receive track %d", iTracks);
1306             continue;
1307         }
1308         
1309         fHistTrackCutStats->Fill("Total",1);
1310         fHistEtaVsZvCoverage->Fill(lVertexZ,track->Eta());
1311         
1312         lPt = track->Pt();
1313         //if ( iTracks == 0 ) cout << "pt of 1st track is " << lPt << endl;
1314         lEta = track->Eta();
1315         lPhi = track->Phi();
1316         
1317         if ( !lRunKine && fAnalysisLevel == "AOD" )
1318         {
1319             AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1320             if ( aodTrack->TestFilterBit( fAODtrackCutBit ) )
1321                 fHistEtaAODpure->Fill( lEta );
1322         }
1323         
1324         //1.10.2012 - "inefficiency" in phi-acceptance "by hand"
1325         if ( lPhi > fPhiArtificialGapBegin &&
1326              lPhi < fPhiArtificialGapEnd )
1327             continue; // drop track in this "ineffective" area
1328         
1329         
1330         //lPhi += phiAdditional; // add here extra phi angle! (when applying tracks rotation by hand)
1331         
1332         
1333         //        if ( lPhi > 2 * TMath::Pi() )
1334         //            lPhi -= 2 * TMath::Pi();
1335         //cout << "track pt = " <<  lPt << endl;
1336         lCharge = track->Charge();
1337         
1338         Int_t lTriggerMask = -1; //for manual event tree extraction
1339         
1340         // ESD or AOD track cuts
1341         if( !lRunKine )
1342         {
1343             if ( fAnalysisLevel == "ESD" )
1344             {
1345                 AliESDtrack* lESDtrack = lESD->GetTrack(iTracks);
1346                 if( fShowPerTrackStats )
1347                     Printf("ESD Track N%d , Eta=%f, Pt=%f , TPC Nclusters =%d Sigma=%f ",  iTracks , lEta , lPt, lESDtrack-> GetTPCNcls(),fEsdTrackCuts->GetSigmaToVertex( lESDtrack) );
1348                 //cluster histograms (for watching)
1349                 fHistClustersTPC->Fill( lESDtrack->GetTPCNcls() );
1350                 fHistCrossedRowsTPC->Fill( lESDtrack->GetTPCCrossedRows() );
1351                 Int_t nITSclusters = lESDtrack->GetITSclusters(0);
1352                 fHistClustersITS->Fill( nITSclusters );//kalmanTrack->GetNumberOfClusters()  );
1353                 //Printf( " after \n" );
1354                 
1355                 fHist2DClustersTPCvsPt->Fill( lPt, lESDtrack->GetTPCNcls() );
1356                 fHist2DClustersTPCvsEta->Fill( lEta, lESDtrack->GetTPCNcls() );
1357                 
1358                 //now check track cuts: either look at fEsdTrackCuts or take cuts from array and look
1359                 if ( !fSwitchToListingCuts )
1360                 {
1361                     if( !fEsdTrackCuts->AcceptTrack(lESDtrack) )
1362                     {
1363                         lCutfail++;
1364                         if( fShowPerTrackStats )
1365                             Printf("Rejected by cuts");
1366                         fHistTrackCutStats->Fill( "QA track cuts", 1 );
1367                         continue;
1368                     }
1369                     else
1370                     {
1371                         if(fShowPerTrackStats)
1372                             Printf("OK");
1373                     }
1374                 }
1375                 else
1376                 {
1377                     lTriggerMask = 0;
1378                     //                    cout << "cuts: " << endl;
1379                     for( Int_t cutsId = 0; cutsId < fArrTrackCuts.GetEntries()/*fNumberOfCutsToRemember*/; cutsId++ )
1380                     {
1381                         //                        if ( !fArrTrackCuts[cutsId] )
1382                         //                            continue;
1383                         int cutsPassed = ( ((AliESDtrackCuts*)fArrTrackCuts[cutsId])->AcceptTrack( lESDtrack ) );
1384                         if ( cutsPassed )
1385                         {
1386                             fHistCutsNamesBins->Fill( cutsId );
1387                             Int_t cutMask = 1;
1388                             for( int iPow = 0; iPow < cutsId; iPow++ )
1389                                 cutMask *= 2;
1390                             lTriggerMask ^=  cutMask;//(cutsId+1);
1391                         }
1392                         //                        cout << cutsPassed;
1393                     }
1394                 }
1395                 
1396                 //fHist2DRejectedTracksPtvsEta->Fill( lPt, lEta );
1397                 fHistClustersTPCafterCuts->Fill( lESDtrack->GetTPCNcls() );
1398                 fHistCrossedRowsTPCafterCuts->Fill( lESDtrack->GetTPCCrossedRows() );
1399             }
1400             else if ( fAnalysisLevel == "AOD" )
1401             {
1402                 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1403                 if (!aodTrack)
1404                 {
1405                     AliError(Form("Could not receive track %d", iTracks));
1406                     continue;
1407                 }
1408                 
1409                 for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
1410                     fHistAODTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
1411                     //                    cout << aodTrack->TestFilterBit(1<<iTrackBit) << " ";
1412                 }
1413                 //                cout << endl;
1414                 if( !aodTrack->TestFilterBit( fAODtrackCutBit ) )
1415                 {
1416                     //                    cout << "fAODtrackCutBit=" << fAODtrackCutBit << endl;
1417                     fHistTrackCutStats->Fill( "QA track cuts", 1 );
1418                     continue;
1419                 }
1420                 //                int a;
1421                 //                cin >> a;
1422             }
1423             
1424             fHistTracksCharge->Fill(lCharge);
1425
1426         }
1427         //        lEta = fRand->Uniform(-1,1); //!!!!!!!!!!!!!!!!!!!!!!! TEST!!!!!!
1428         // end of ESD or AOD track cuts
1429         
1430         
1431         //MC truth
1432         //Int_t label = -1000;
1433         TParticle * part = 0x0;
1434         if( lRunKine )   // Dropping undetectable particles in Kine
1435         {
1436             //            if ( fAnalyseMCESD ) //run MC+ESD
1437             //            {
1438             //                label = track->GetLabel();//TMath::Abs(track->GetLabel()); //abs!!!!
1439             //                fHistMClabels->Fill( label );
1440             //                if ( label < 0 ) //by Peter Hristov - it's ghost tracks
1441             //                    continue;
1442             //                part = stack->Particle(label);
1443             //            }
1444             //            else
1445             {
1446                 AliMCParticle* trackMC = dynamic_cast<AliMCParticle *>(eventMC->GetTrack(iTracks));
1447                 if(!trackMC)
1448                     continue;
1449                 part = trackMC->Particle();
1450                 if( !part ) continue;
1451             }
1452             
1453             Int_t lNoD = part->GetNDaughters();
1454             if(fShowPerTrackStats)
1455             {
1456                 printf("Track %d, Mother %d, track Pt=%f, MC Pt=%f, PDG=%d  Nof Dothers=%d  ETA=%f   ## ",  iTracks , part->GetFirstMother() ,lPt,      part->Pt() ,part->GetPdgCode(),lNoD,lEta);
1457                 printf("%s", part->GetName());
1458                 printf(" ");
1459                 printf("%s", part->GetPDG()->ParticleClass());
1460             }
1461             
1462             if( !stack->IsPhysicalPrimary(iTracks) )//label))
1463             {
1464                 fHistTrackCutStats->Fill( "Not phys primary", 1 );
1465                 continue;
1466             }
1467             
1468             //charge in TParticlePDG is in units of |e|/3
1469             fHistTracksCharge->Fill( part->GetPDG()->Charge() / 3 ); // 0-charge bin fills only for MC truth events
1470             //cout << " charge = " << part->GetPDG()->Charge();
1471             if ( part->GetPDG()->Charge() == 0 )
1472             {
1473                 if(fShowPerTrackStats)
1474                     Printf(" ChargeReject");
1475                 lNoCharge++;
1476                 continue;
1477             }
1478             
1479             //cut on low-pt kine particles
1480             if ( lPt < fKineLowPtCut )
1481                 continue;
1482         }    //End of  Kine particle filter
1483         
1484         //now decide that we have a good track
1485         lNaccept++;
1486         
1487         if( fabs(lEta) < fEtaRegionForTests ) //look at eta region
1488         {
1489             lNacceptEtaInRegion++;
1490             //net charge fillings
1491             fHistNetChargeVsPt->Fill( lPt, lCharge );
1492             if ( lCharge > 0 )
1493             {
1494                 fHistChargePlusVsPtTmp->Fill( lPt );
1495                 fHistPtPlus->Fill( lPt );
1496             }
1497             else if ( lCharge < 0 )
1498             {
1499                 fHistChargeMinusVsPtTmp->Fill( lPt );
1500                 fHistPtMinus->Fill( lPt );
1501             }
1502         }
1503         
1504         if( lPt > fMaxPtLimit )
1505         {
1506             lPtOver++;
1507             fHistTrackCutStats->Fill(2);
1508             continue;
1509         } // Dropping tracks with hi Pt
1510         
1511         if( lPt < fMinPtLimit )
1512         {
1513             lPtUnder++;
1514             fHistTrackCutStats->Fill(3);
1515             continue;
1516         } // Dropping tracks with low Pt
1517         lNacceptAfterPtCuts++;
1518         fHistTrackCutStats->Fill(4);
1519         fHistPt->Fill( lPt );
1520         fHistEta->Fill( lEta );
1521         fHistPhi->Fill( lPhi );
1522         
1523         if( fabs(lEta) < fEtaRegionForTests ) //look at eta region
1524             lNacceptEtaInRegionAfterPtCuts++;
1525         
1526         //fill counters for ZDC-comparing with a 'shelve'
1527         if ( lEta > etaBmin && lEta < etaBmax )
1528             countTracksEtaB++;
1529         if ( lEta > etaFmin && lEta < etaFmax )
1530             countTracksEtaF++;
1531         
1532         if ( fabs( lEta ) < 0.8 ) //eta-phi for tracks in "ALICE barrel acceptance"
1533             fHistEtaPhi->Fill(lEta, lPhi);
1534         
1535         
1536         // ###### New PID
1537         Int_t lMostProbablePIDPure = -1;
1538         Int_t lMostProbablePIDdirty = -1;
1539         if ( (lESD || lAOD) && !lRunKine && fPIDResponse )
1540         {
1541             AliVTrack *trackV = (AliVTrack*)event->GetTrack(iTracks);
1542             //##### from Pid Combined Task
1543             // compute priors for TPC+TOF, even if we ask just TOF for PID
1544             fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
1545             UInt_t detUsed = fPIDCombined->ComputeProbabilities(trackV, fPIDResponse, probTOF);
1546             Double_t priors[5];         // check priors used for TOF
1547             fPIDCombined->GetPriors(trackV,priors,fPIDResponse,detUsed);
1548             for(Int_t ispec=0;ispec<5;ispec++) fPriorsUsed[ispec]->Fill(TMath::Abs(trackV->Pt()),priors[ispec]);
1549             
1550             fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
1551             detUsed = fPIDCombined->ComputeProbabilities(trackV, fPIDResponse, probTPCTOF);
1552             fHistProbabilityPion->Fill( probTPCTOF[AliPID::kPion] );////probDensity[iParticle] );
1553             fHistProbabilityKaon->Fill( probTPCTOF[AliPID::kKaon] );
1554             fHistProbabilityProton->Fill( probTPCTOF[AliPID::kProton] );////probDensity[iParticle] );
1555             //cout << "chto-to s detId: " << detUsed << endl;
1556 //            for ( Int_t ispec=0; ispec < (int)(AliPID::kSPECIES); ++ispec ) {
1557 //                //cout << probTPCTOF[ispec] << "  " << endl;
1558 //            }
1559             //cout << endl;
1560             if ( detUsed == (UInt_t)fPIDCombined->GetDetectorMask() )
1561             {
1562                 Double_t lMaxProb = 0;
1563                 for ( Int_t ispec=0; ispec < (int)(AliPID::kSPECIES); ++ispec )
1564                 {
1565                     //Double_t nSigmaTPC = fPIDResponse->NumberOfSigmasTPC(trackV,(AliPID::EParticleType)ispec);
1566                     fProbTPCTOF[ispec]->Fill( lPt, probTPCTOF[ispec] );
1567                     //fProbTPCTOFnSigmaTPC[ispec]->Fill(nSigmaTPC,probTPCTOF[ispec]);
1568                     //fProbTPCTOFnSigTPCMom[ibin][ispec]->Fill(nSigmaTPC,probTPCTOF[ispec]);
1569                     
1570                     if ( probTPCTOF[ispec] > lMaxProb )
1571                     {
1572                         lMaxProb = probTPCTOF[ispec];
1573                         lMostProbablePIDPure = ispec;
1574                     }
1575                 }
1576                 //fill for watching at max prob-s in pid species arrays
1577                 fHistPidPureMaxProbability->Fill( lMaxProb );
1578             }
1579             
1580             const bool useDirtyPID = true;
1581             if ( useDirtyPID )
1582             {
1583                 //Int_t lMostProbablePIDdirty = -1;
1584                 
1585                 Double_t lMaxProb = 0;
1586                 for ( Int_t ispec=0; ispec < (int)(AliPID::kSPECIES); ++ispec )
1587                 {
1588                     fProbAllDets[ispec]->Fill( lPt, probTPCTOF[ispec] );
1589                     if ( probTPCTOF[ispec] > lMaxProb )
1590                     {
1591                         lMaxProb = probTPCTOF[ispec];
1592                         lMostProbablePIDdirty = ispec;  // define most probable particle!!!
1593                     }
1594                 }
1595                 //fill for watching at max prob-s in pid species arrays
1596                 fHistPidMaxProbability->Fill( lMaxProb );
1597                 fHistParticlesDistr->Fill( lMostProbablePIDdirty );
1598             }
1599             
1600             
1601         }
1602         // end of new PID filling
1603         
1604         fHist2DAcceptedTracksPtvsEta->Fill( lPt, lEta );
1605         
1606         //get particle mass and Et
1607         if ( lMostProbablePIDPure >= 0 )
1608         {
1609             lMass = AliPID::ParticleMass( lMostProbablePIDPure );
1610             fHistESDtrackMass->Fill( lMass );
1611             lEt = sqrt( lPt*lPt + lMass*lMass );
1612             //cout << lEt << " " << lMass << " " << lPt << endl;
1613         }
1614         
1615         //eta vs vertex z coverage
1616         fHistEtaVsZvCoverageAccepted->Fill(lVertexZ,track->Eta());
1617         
1618         //artificial inefficiency (27.09.12) - check, what if we drop some fraction of tracks randomly
1619         Bool_t lDropDecision = ( ( 1.0 - fArtificialInefficiency ) < fRand->Uniform(0,1) );
1620         
1621         if ( lDropDecision )
1622         {
1623             lNumberOfDropByHandTracks++;
1624             continue;
1625         }
1626         
1627         if ( lNumberOfAcceptedTracksForLRC >= kMaxParticlesNumber )
1628         {
1629             AliDebug(AliLog::kError, "lNumberOfAcceptedTracksForLRC is too large...");
1630             break;
1631         }
1632         
1633         //fill arrays with track data for LRC
1634         fArrayTracksPt[lNumberOfAcceptedTracksForLRC]       = (fEtInsteadOfPt ? lEt : lPt);
1635         fArrayTracksEta[lNumberOfAcceptedTracksForLRC]      = lEta;
1636         fArrayTracksPhi[lNumberOfAcceptedTracksForLRC]      = lPhi;
1637         fArrayTracksCharge[lNumberOfAcceptedTracksForLRC]   = lCharge;
1638         fArrayTracksPID[lNumberOfAcceptedTracksForLRC]      = lMostProbablePIDPure;
1639         
1640         lNumberOfAcceptedTracksForLRC++;
1641         
1642     } //end of track loop
1643     
1644     //net charge vs pt
1645 //    if ( lNacceptEtaInRegion > 20 )
1646 //    {
1647 //        fHistChargePlusVsPtTmp->DrawClone();
1648 //        fHistChargeMinusVsPtTmp->SetLineColor( kBlue );
1649 //        fHistChargeMinusVsPtTmp->DrawClone("same");
1650 //    }
1651
1652
1653     //final actions for net-charge
1654     double integralPlus     = fHistChargePlusVsPtTmp->ComputeIntegral();
1655     double integralMinus    = fHistChargeMinusVsPtTmp->ComputeIntegral();
1656     double xAxisLowEdge     = fHistChargePlusVsPtTmp->GetBinLowEdge(1);
1657     double xAxisHighEdge    = fHistChargePlusVsPtTmp->GetBinLowEdge(fHistChargePlusVsPtTmp->GetNbinsX()) + fHistChargePlusVsPtTmp->GetBinWidth(fHistChargePlusVsPtTmp->GetNbinsX());
1658     double ptRange =  xAxisHighEdge - xAxisLowEdge;
1659 //    cout << "integralPlus=" << integralPlus << ", integralMinus=" << integralMinus << ", ptRange=" << ptRange << endl;
1660
1661     for ( int bin = 0; bin < fHistChargePlusVsPtTmp->GetNbinsX(); bin++)
1662     {
1663         double binCenter = fHistChargePlusVsPtTmp->GetBinCenter(bin+1);
1664         double binContentPlus = fHistChargePlusVsPtTmp->GetBinContent(bin+1);
1665         double binContentMinus = fHistChargeMinusVsPtTmp->GetBinContent(bin+1);
1666
1667
1668
1669         double lNetCharge = binContentPlus-binContentMinus;
1670         fHist2DNetChargeVsPt->Fill( binCenter, lNetCharge );
1671         if ( ptRange>0 )
1672             fHist2DNetChargeVsPtCorrectedOnEventMean->Fill( binCenter, lNetCharge - (integralPlus-integralMinus)/ptRange );
1673         //        fHist2DNetChargeVsPtCorrectedOnEventMeanNormOnNch->Fill( );
1674
1675     }
1676
1677     //############ MC ESD QA
1678     if ( fRunKine && fAnalyseMCESD ) //run MC+ESD
1679     {
1680         Double_t lEta;
1681         Double_t lPt;
1682 //        Double_t lY;
1683         Int_t lMCtracksAcceptedEtaInRegionAfterPtCuts = 0;
1684         for (Int_t iMCtracks = 0; iMCtracks < eventMC->GetNumberOfPrimaries(); iMCtracks++)
1685         {
1686             AliMCParticle* track = dynamic_cast<AliMCParticle *>(eventMC->GetTrack(iMCtracks));
1687             if (!track) {
1688                 Printf("ERROR: Could not receive particle %d", iMCtracks);
1689                 continue;
1690             }
1691             
1692             //exclude non stable particles
1693             if(!eventMC->IsPhysicalPrimary(iMCtracks)) continue;
1694             
1695             if ( track->Particle()->GetPDG()->Charge() == 0 )
1696                 continue;
1697             
1698             lEta    = track->Eta();
1699             lPt     = track->Pt();
1700 //            lY      = track->Y();
1701             
1702             if( fabs(lEta) > fEtaRegionForTests ) //look at eta region
1703                 continue;
1704             
1705             if( lPt > fMaxPtLimit || lPt < fMinPtLimit )
1706                 continue;
1707             
1708             
1709             lMCtracksAcceptedEtaInRegionAfterPtCuts++;
1710         }
1711         fHist2DMultiplicityMCESDInEtaRegion->Fill( lMCtracksAcceptedEtaInRegionAfterPtCuts, lNacceptEtaInRegionAfterPtCuts );
1712     }
1713     
1714     
1715     
1716     
1717     //profiling tracks in phi BY HAND if requested
1718     //    if ( fUsePhiShufflingByHand )
1719     //    {
1720     //        double lFictionEventPlanePhi = fRand->Uniform(0,2*TMath::Pi());
1721     //        fHistPhiArtificialEvPlane->Fill( lFictionEventPlanePhi );
1722     //        TF1 lFictionEventPlaneFunc( "fictionEventPlane", "0.75 + 0.25*TMath::Cos(x)", 0, 2*TMath::Pi() );
1723     //        for ( Int_t trackId = 0; trackId < lNumberOfAcceptedTracksForLRC; trackId++ )
1724     //        {
1725     //            //check phi wrt event plane
1726     //            Double_t lProfiledPhiWrtEvPlane = lFictionEventPlaneFunc.GetRandom();
1727     //            //FixAngleInTwoPi( lProfiledPhiWrtEvPlane );
1728     //            fHistPhiArtificialProfilingCheckWrtEvPlane->Fill( lProfiledPhiWrtEvPlane );
1729     
1730     //            //  double lOppositePhi = ( ( fRand->Uniform(0,1) > 0.5 ) ? -TMath::Pi() : 0 );
1731     //            //  double lRandomConeAngleEta = fRand->Gaus(0,0.5);
1732     //            //  double lRandomConeAnglePhi = fRand->Gaus(0,TMath::PiOver2());
1733     
1734     //            //double lFictionPhi = fRand->Uniform(0,2*TMath::Pi());
1735     //            //double lPhiSign = ( ( fRand->Uniform(0,1) > 0.5 ) ? TMath::Pi() : 0 );
1736     
1737     //            //add event plane phi angle
1738     //            Double_t lProfiledPhi = lProfiledPhiWrtEvPlane + lFictionEventPlanePhi;
1739     //            FixAngleInTwoPi( lProfiledPhi );
1740     //            fHistPhiArtificialProfilingCheck->Fill( lProfiledPhi );
1741     
1742     //            fArrayTracksPhi[trackId] = lProfiledPhi; //lPhi + lRandomConeAnglePhi; //lPhi;
1743     
1744     //            //        fArrayTracksPhi[lNumberOfAcceptedTracksForLRC]      = lPhi - lRandomConeAnglePhi; //lPhi;
1745     //            //        if ( fArrayTracksPhi[lNumberOfAcceptedTracksForLRC] > 2 * TMath::Pi() )
1746     //            //            fArrayTracksPhi[lNumberOfAcceptedTracksForLRC] -= 2 * TMath::Pi();
1747     //        }
1748     //    }
1749     
1750     
1751     // ######### fill some QA plots
1752     fHistAcceptedTracks->Fill(lNaccept);
1753     fHistMultiplicityInEtaRegion->Fill( lNacceptEtaInRegion );
1754     fHistMultiplicityInEtaRegionAfterPtCuts->Fill( lNacceptEtaInRegionAfterPtCuts );
1755     fHistAcceptedTracksAfterPtCuts->Fill(lNacceptAfterPtCuts);
1756     if ( fArtificialInefficiency >= 0 ) //i.e. have ineff analysis ON
1757         fHistNumberOfDroppedByHandTracks->Fill( lNacceptAfterPtCuts, lNumberOfDropByHandTracks );
1758     
1759     // ########### ProfilePhiByHand
1760     if ( fUsePhiShufflingByHand )
1761         ProfilePhiByHand( lNumberOfAcceptedTracksForLRC);
1762     
1763     // ########### LRC processors filling
1764     Bool_t lFlagAcceptEventForLRC = ( lNacceptEtaInRegion >= fMultCutInEtaRegion );
1765     if ( lFlagAcceptEventForLRC )
1766         FillLRCProcessors( lNumberOfAcceptedTracksForLRC, lCentralityPercentile );
1767     
1768     // ##### ZDC, VZERO, etc
1769     if ( fFlagWatchZDC && fIsIonsAnalysis )
1770     {
1771         AliESDZDC *zdcData = lESD->GetESDZDC();
1772         //fHistZDCenergy->Fill( zdcData->GetZDCEMEnergy(0) );
1773         fHistZDCparticipants->Fill( zdcData->GetZDCParticipants ());
1774         for ( int i = 0; i < 5; i++ )
1775         {
1776             //if multiplicity in eta windows above threshold - fill ZDC energy hist
1777             if ( countTracksEtaB >= fMultForZDCstudy[i] && countTracksEtaF >= fMultForZDCstudy[i] )
1778                 fHistZDCenergy[i]->Fill( zdcData->GetZDCEMEnergy(0) );
1779         }
1780     }
1781     if ( fFlagWatchV0 )
1782     {
1783         const AliESDVZERO* vzrData = lESD->GetVZEROData(); //aod the same
1784
1785         
1786         //V0 multiplicity
1787         float lV0mult[64];
1788         //float lV0Amult[64];
1789         //float lV0Cmult[64];
1790         double sumV0mult = 0;
1791         double sumV0Amult = 0;
1792         double sumV0Cmult = 0;
1793         
1794         for (int i = 0; i < 64; i++ )
1795         {
1796             lV0mult[i] = (float)(vzrData->GetMultiplicity(i));
1797             sumV0mult += lV0mult[i];
1798         }
1799         fHistV0multiplicity->Fill( sumV0mult );
1800         
1801         //            for (int i=0; i<32; i++)
1802         //            {
1803         //                lV0Cmult[i] = (float)(vzrData->GetMultiplicityV0C(i));
1804         //                sumV0Cmult += lV0Cmult[i];
1805         //                lV0Amult[i] = (float)(vzrData->GetMultiplicityV0A(i));
1806         //                sumV0Amult += lV0Amult[i];
1807         //            }
1808         sumV0Amult = vzrData->GetMTotV0A();
1809         sumV0Cmult = vzrData->GetMTotV0C();
1810         fHistV0Amultiplicity->Fill( sumV0Amult );
1811         fHistV0Cmultiplicity->Fill( sumV0Cmult );
1812         fHist2DV0ACmultiplicity->Fill( sumV0Amult, sumV0Cmult );
1813         
1814         fHist2DTracksAcceptedVsV0multiplicity->Fill( lNaccept, sumV0Amult );
1815         
1816         //cells
1817         //fHistV0cells->Fill(  );
1818         fHistV0Acells->Fill( vzrData->GetNbPMV0A() );
1819         fHistV0Ccells->Fill( vzrData->GetNbPMV0C() );
1820         fHist2DV0ACcells->Fill( vzrData->GetNbPMV0A(), vzrData->GetNbPMV0C() );
1821         
1822         //rings
1823         for ( int i = 0; i < 4; i++ )
1824         {
1825             int lMultRingV0A = vzrData->GetMRingV0A(i);
1826             int lMultRingV0C = vzrData->GetMRingV0C(i);
1827             fHistV0AmultiplicityRing[i]->Fill( lMultRingV0A );
1828             fHistV0CmultiplicityRing[i]->Fill( lMultRingV0C );
1829             fHist2DV0ACmultiplicityRing[i]->Fill( lMultRingV0A, lMultRingV0C );
1830             
1831             fHist2DTracksAcceptedVsV0AmultiplicityRing[i]->Fill( lNaccept, lMultRingV0A );
1832             fHist2DTracksAcceptedVsV0CmultiplicityRing[i]->Fill( lNaccept, lMultRingV0C );
1833         }
1834         
1835     }
1836     if ( fFlagWatchFMD )
1837     {
1838         // multiplicity
1839         //const AliESDFMD *lFMDdata = lESD->GetFMDData();
1840         
1841         
1842     }
1843     
1844     //        //look at V0s
1845     //        AliESDv0 *lV0 = 0x0;
1846     //        //cout << ">>>> showing v0's:" << endl;
1847     //        for ( int iV0 = 0; iV0 < lESD->GetNumberOfV0s(); iV0++ )
1848     //        {
1849     //            lV0 = lESD->GetV0( iV0 );
1850     
1851     //            int lV0pdgCode = lV0->GetPdgCode();
1852     //            double lMassV0 = lV0->GetEffMass();//ChangeMassHypothesis( lV0pdgCode );  //GetEffMass();
1853     //            //if ( abs(lV0pdgCode) > 10 )//== 3122 )
1854     //            {
1855     //                //cout << ">>>>>>>>>>>>>>>>>>>>>>>>> " << lV0pdgCode << ": mass = " << lMassV0 << endl;
1856     //                fHistV0spectra->Fill(lV0pdgCode);//lMassV0);
1857     //            }
1858     
1859     //            if ( lMassV0 > 1. )// abs(lPdgCode) == 3122 || abs(lPdgCode) == 421 ) //abs(lPdgCode) % 1000 == 122 )
1860     //            {
1861     //                //cout << ">>>>>>>>>>>>>>>>>>>>>>>>> " << lV0pdgCode << ": mass = " << lMassV0 << endl;
1862     //            }
1863     //        }
1864     
1865     
1866     
1867     
1868     
1869     
1870     //Debuging output of track filter
1871     if( fShowPerTrackStats )
1872         Printf("NofTracks= %d , accepted %d , LowPt %d, HighPt %d, LowCharge %d,  ESD cut fail %d , Decay Filer %d", event->GetNumberOfTracks(), lNaccept, lPtUnder, lPtOver ,lNoCharge , lCutfail  , lDecay);
1873     
1874     // Post output data
1875     PostData(1, fOutList);
1876     
1877     for( Int_t i = 0; i < fLRCproc.GetEntries(); i++)
1878     {
1879         PostData( Proc(i)->GetOutputSlotNumber(),Proc(i)->CreateOutput() );
1880     }
1881     
1882 }
1883
1884
1885 void AliAnalysisTaskLRC::SetParticleTypeToProcessors( int windowId, TString strPid ) //char* strPid )//AliLRCBase::LRCparticleType whichParticleToFill )
1886 {
1887     //dummy actions
1888     windowId+=1;
1889     if(0) printf("%s", strPid.Data());
1890
1891     //Set pid types for LRC processors (windowId = 0 - fwd window, =1 - backward window)
1892     //    Int_t lLrcNum=fLRCproc.GetEntries();
1893     
1894     //    for(Int_t i=0; i < lLrcNum; i++)
1895     //    {
1896     //        AliLRCBase *lrcProc = dynamic_cast<AliLRCBase*> (fLRCproc.At(i));
1897     //        if(lrcProc)
1898     //        {
1899     //            //tmp lines!
1900     //            windowId++;
1901     //            if (0)
1902     //                printf("%s",strPid);
1903     //            //commented for this commit
1904     //            //            if ( windowId == 0 )
1905     //            //                lrcProc->SetParticleType( "fwd", strPid );
1906     //            //            else if ( windowId == 1 )
1907     //            //                lrcProc->SetParticleType( "bkwd", strPid );
1908     
1909     //            //fTmpCounter++;
1910     //        }
1911     //        else continue;
1912     //        //fTmpCounter++;
1913     //    }
1914     
1915 }
1916
1917 void AliAnalysisTaskLRC::SetParticleTypeForTask( TString strF, TString strB ) //char* strF, char* strB )
1918 {
1919     // Set PID sensitivity to 'true' and write pids for windows
1920     fPIDsensingFlag = kTRUE;
1921     fStrPIDforFwd = strF;
1922     fStrPIDforBwd = strB;
1923
1924     //    sprintf ( fStrPIDforFwd, "%s", strF );
1925     //    sprintf ( fStrPIDforBwd, "%s", strB );
1926     //cout << "we just set fwd win to " << fStrPIDforFwd << " and bwd win to " << fStrPIDforBwd << endl;
1927 }
1928
1929
1930
1931 //________________________________________________________________________
1932 void AliAnalysisTaskLRC::Terminate(Option_t *)
1933 {
1934     // Draw result to the screen
1935     // Called once at the end of the query
1936     //fOutList = dynamic_cast<TList*> (GetOutputData(0));
1937     
1938     //lESDTrackCuts->DrawHistograms();
1939     
1940     
1941     fAnalysisTimer->Stop();
1942     Double_t rtime = fAnalysisTimer->RealTime();
1943     Double_t ctime = fAnalysisTimer->CpuTime();
1944     
1945     printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
1946     
1947 }
1948
1949
1950 void AliAnalysisTaskLRC::AddLRCProcess(AliLRCBase *newProc)
1951 {
1952     // Add new AliLRCBase (Main LRC processor per ETA window) to the processing list
1953     // Used to add new ETA window to AnalysisTask
1954     if(!newProc)
1955     {
1956         Printf("ERROR:No AliLRCBase object -  NULL pointer!");
1957         return;
1958     }
1959     
1960     fLRCproc.Add(newProc);
1961     newProc->SetOutputSlotNumber(fLRCproc.GetEntries() + 1 );//( fSetIncludeEventTreeInOutput ? 2 : 1 ) );
1962     DefineOutput(newProc->GetOutputSlotNumber(),TList::Class());
1963     return ;
1964 }
1965
1966
1967 Double_t AliAnalysisTaskLRC::GetEventPlane(AliVEvent *event)
1968 {
1969     // Get the event plane
1970     
1971     if ( !event )
1972         return -1;
1973     //TString gAnalysisLevel = fBalance->GetAnalysisLevel();
1974     
1975     Float_t lVZEROEventPlane    = -10.;
1976     Float_t lReactionPlane      = -10.;
1977     Double_t qxTot = 0.0, qyTot = 0.0;
1978     
1979     //MC: from reaction plane
1980     if( fRunKine )//gAnalysisLevel == "MC"){
1981     {
1982         AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader());
1983         if (headerH)
1984         {
1985             lReactionPlane = headerH->ReactionPlaneAngle();
1986             //lReactionPlane *= TMath::RadToDeg();
1987         }
1988     }//MC
1989     
1990     // AOD,ESD,ESDMC: from VZERO Event Plane
1991     else
1992     {
1993         AliEventplane *ep = event->GetEventplane();
1994         if( ep )
1995         {
1996             lVZEROEventPlane = ep->CalculateVZEROEventPlane( event, 10, 2, qxTot, qyTot);
1997             if( lVZEROEventPlane < 0. )
1998                 lVZEROEventPlane += TMath::Pi();
1999             //lReactionPlane = lVZEROEventPlane*TMath::RadToDeg();
2000             lReactionPlane = lVZEROEventPlane;
2001             //            cout << "lReactionPlane: " << lReactionPlane << endl;
2002             //            int aa;
2003             //            cin >> aa;
2004         }
2005     }//AOD,ESD,ESDMC
2006     
2007     return lReactionPlane;
2008 }
2009
2010
2011 void AliAnalysisTaskLRC::AddTrackCutForBits(AliESDtrackCuts*  const cuts, TString cutsName )
2012 {
2013     //TString *strPtrCutsName = new TString(cutsName.Data());
2014 //    cout << "ADDING CUTS " << cutsName << endl;
2015     //    fArrTrackCuts[fNumberOfCutsToRemember] = cuts;
2016     int lNumberOfCutsToRemember = fArrTrackCuts.GetEntries();
2017     
2018     fArrCutsNames[lNumberOfCutsToRemember] = cutsName;
2019     fArrTrackCuts.Add( cuts );// = cuts;
2020     //    fArrCutsNames.Add( strPtrCutsName );
2021     //    fNumberOfCutsToRemember++;
2022 }
2023
2024 void AliAnalysisTaskLRC::UseToyEvents()
2025 {
2026     //generate all events here and fill LRC processors
2027     for ( Int_t toyEventId = 0; toyEventId < fNumberOfToyEvents; toyEventId++ )
2028     {
2029         if ( toyEventId % 10000 == 0 )
2030             printf("Processing %d event...\n", toyEventId );
2031
2032         //        cout << "toyEventId=" << toyEventId << endl;
2033         int nToyTracks =  0;
2034         const double kExpParam = 0.42;
2035         const double kEtaGausSigma = 0.3;
2036         const double kPhiGausSigma = TMath::PiOver4();
2037
2038         //try to do "jets" with some particles
2039         const int nJets = 2;
2040
2041         //jets eta angles
2042         double lJets_Eta[nJets];
2043         lJets_Eta[0] = fRand->Uniform(0, 0.8);
2044         lJets_Eta[1] = -lJets_Eta[0];
2045
2046         //jets phi angles
2047         double lJets_Phi[nJets];
2048         lJets_Phi[0] = fRand->Uniform(0, 2*TMath::Pi());
2049         lJets_Phi[1] = lJets_Phi[0] + TMath::Pi();
2050         FixAngleInTwoPi( lJets_Phi[1] );
2051
2052
2053         //jets loop
2054         for ( Int_t jetId = 0; jetId < nJets; jetId++ )
2055         {
2056             int nPartInJet = fRand->Gaus(4,1);
2057             //jet tracks params
2058             for ( Int_t toyTrackId = 0; toyTrackId < nPartInJet; toyTrackId++ )
2059             {
2060                 fArrayTracksPt[nToyTracks] = fRand->Exp( kExpParam );
2061                 fArrayTracksEta[nToyTracks] = fRand->Gaus( lJets_Eta[jetId], kEtaGausSigma );
2062                 fArrayTracksPhi[nToyTracks] = fRand->Gaus( lJets_Phi[jetId], kPhiGausSigma );
2063                 FixAngleInTwoPi( fArrayTracksPhi[nToyTracks] );
2064
2065                 nToyTracks++;
2066             }
2067         }
2068
2069         FillLRCProcessors( nToyTracks, 0 );
2070
2071         //        for ( Int_t toyTrackId = 0; toyTrackId < nToyTracks; toyTrackId++ )
2072         //        {
2073         //            // ########### ProfilePhiByHand
2074         //            ProfilePhiByHand( nToyTracks);
2075         //        }
2076
2077     }
2078 }
2079
2080
2081 void AliAnalysisTaskLRC::FillLRCProcessors( int numberOfAcceptedTracksForLRC, Double_t eventCentrality )
2082 {
2083     //pass signal to LRC-based analysis to start event
2084     Int_t lLrcNum = fLRCproc.GetEntries(); // Number of processors attached
2085     
2086     //prepare phi rotation step
2087     Double_t phiStep = 2 * TMath::Pi();
2088     if ( fNumberOfPhiSectors > 1 )
2089         phiStep /= fNumberOfPhiSectors;
2090     Double_t lPhiRotatedExtra = 0; //additional phi rotation of all tracks in case
2091     
2092     
2093     //start LRC filling
2094     for ( Int_t sectorId = 0; sectorId < fNumberOfPhiSectors; sectorId++ )
2095     {
2096         if ( lLrcNum > kMaxParticlesNumber )
2097         {
2098             AliDebug(AliLog::kError, "lLrcNum is too large... Mistake???");
2099             break;
2100         }
2101         
2102         //pass signal to LRC-based analysis to start event
2103         for( Int_t lrcProcessorId = 0; lrcProcessorId < lLrcNum; lrcProcessorId++ )
2104         {
2105             AliLRCBase *lrcBase = dynamic_cast<AliLRCBase*> (fLRCproc.At(lrcProcessorId));
2106             if (!lrcBase)
2107             {
2108                 AliDebug(AliLog::kError, "can't cast lrcBase to AliLRCBase!" );
2109                 break;
2110             }
2111             lrcBase->StartEvent();
2112             //pass the centrality
2113             lrcBase->SetEventCentrality( eventCentrality );
2114             //        }
2115             
2116             //pass track data to LRC-based analysis
2117             for ( Int_t trackId = 0; trackId < numberOfAcceptedTracksForLRC; trackId++ )
2118             {
2119                 if ( fEtInsteadOfPt && (fArrayTracksPt[trackId]  == 0 ) ) //in Et-mode we didn't get Et! exit
2120                 {
2121                     AliDebug(AliLog::kError, "pT=0 Mistake???" );
2122                     break;
2123                 }               //if ( fEtInsteadOfPt && lMostProbablePIDPure == -1 ) //don't have pure PID
2124                 //      continue;
2125                 
2126                 //rotate track phi
2127                 Double_t lPhi = fArrayTracksPhi[trackId] + lPhiRotatedExtra;
2128                 FixAngleInTwoPi( lPhi );
2129                 
2130                 fHistPhiLRCrotationsCheck->Fill( lPhi );
2131                 
2132                 //            for(Int_t i = 0; i < lLrcNum; i++ )
2133                 //            {
2134                 lrcBase->AddTrackPtEta(
2135                             fArrayTracksPt[trackId]
2136                             , fArrayTracksEta[trackId]
2137                             , lPhi
2138                             , fArrayTracksCharge[trackId]
2139                             , fArrayTracksPID[trackId]
2140                             );
2141                 //            }
2142                 
2143             }
2144             
2145             //take event only if at least 1 track in this event fulfill the requirements! //21.11.11
2146             Bool_t lDontTakeEventDecision =  kFALSE; //lNaccept > 0 ? kFALSE : kTRUE;
2147             //pass signal to LRC-based analysis to finish the event
2148             //        for( Int_t i = 0; i < lLrcNum; i++ )
2149             //        {
2150             lrcBase->FinishEvent( lDontTakeEventDecision );
2151         }
2152         
2153         lPhiRotatedExtra += phiStep; //increase phi step to rotate tracks
2154     }
2155 }
2156
2157 void AliAnalysisTaskLRC::ProfilePhiByHand( int numberOfAcceptedTracksForLRC )
2158 {
2159     //profiling tracks in phi BY HAND if requested
2160     
2161     double lFictionEventPlanePhi = fRand->Uniform(0,2*TMath::Pi());
2162     fHistPhiArtificialEvPlane->Fill( lFictionEventPlanePhi );
2163     TF1 lFictionEventPlaneFunc( "fictionEventPlane", "0.75 + 0.25*TMath::Cos(x)", 0, 2*TMath::Pi() );
2164     for ( Int_t trackId = 0; trackId < numberOfAcceptedTracksForLRC; trackId++ )
2165     {
2166         //check phi wrt event plane
2167         Double_t lProfiledPhiWrtEvPlane = lFictionEventPlaneFunc.GetRandom();
2168         //FixAngleInTwoPi( lProfiledPhiWrtEvPlane );
2169         fHistPhiArtificialProfilingCheckWrtEvPlane->Fill( lProfiledPhiWrtEvPlane );
2170         
2171         //  double lOppositePhi = ( ( fRand->Uniform(0,1) > 0.5 ) ? -TMath::Pi() : 0 );
2172         //  double lRandomConeAngleEta = fRand->Gaus(0,0.5);
2173         //  double lRandomConeAnglePhi = fRand->Gaus(0,TMath::PiOver2());
2174         
2175         //double lFictionPhi = fRand->Uniform(0,2*TMath::Pi());
2176         //double lPhiSign = ( ( fRand->Uniform(0,1) > 0.5 ) ? TMath::Pi() : 0 );
2177         
2178         //add event plane phi angle
2179         Double_t lProfiledPhi = lProfiledPhiWrtEvPlane + lFictionEventPlanePhi;
2180         FixAngleInTwoPi( lProfiledPhi );
2181         fHistPhiArtificialProfilingCheck->Fill( lProfiledPhi );
2182         
2183         fArrayTracksPhi[trackId] = lProfiledPhi; //lPhi + lRandomConeAnglePhi; //lPhi;
2184         
2185         //        fArrayTracksPhi[lNumberOfAcceptedTracksForLRC]      = lPhi - lRandomConeAnglePhi; //lPhi;
2186         //        if ( fArrayTracksPhi[lNumberOfAcceptedTracksForLRC] > 2 * TMath::Pi() )
2187         //            fArrayTracksPhi[lNumberOfAcceptedTracksForLRC] -= 2 * TMath::Pi();
2188     }
2189 }