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