]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/hfe/AliAnalysisTaskFlowITSTPCTOFQCSP.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix of sigmaZ for crossing tracklets from Alex
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliAnalysisTaskFlowITSTPCTOFQCSP.cxx
CommitLineData
4c12b1ec
R		 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16
17////////////////////////////////////////////////////////////////////////
18// //
19// Task for Heavy Flavour Electron Flow with TPC plus TOF //
20// Non-Photonic Electron identified with Invariant mass //
21// analysis methos in function SelectPhotonicElectron //
22// //
23// //
24// Author: Andrea Dubla (Utrecht University) //
25// //
26// //
27////////////////////////////////////////////////////////////////////////
28
29#include "TChain.h"
30#include "TTree.h"
31#include "TH2F.h"
32#include "TMath.h"
33#include "TCanvas.h"
34#include "THnSparse.h"
35#include "TLorentzVector.h"
36#include "TString.h"
37#include "TFile.h"
38#include "AliAnalysisTask.h"
39#include "AliAnalysisManager.h"
40#include "AliESDEvent.h"
41#include "AliESDHandler.h"
42#include "AliAODEvent.h"
43#include "AliAODHandler.h"
44#include "AliAnalysisTaskFlowITSTPCTOFQCSP.h"
45#include "TGeoGlobalMagField.h"
46#include "AliLog.h"
47#include "AliAnalysisTaskSE.h"
48#include "TRefArray.h"
49#include "TVector.h"
50#include "AliESDInputHandler.h"
51#include "AliESDpid.h"
52#include "AliAODInputHandler.h"
53#include "AliAODPid.h"
54#include "AliESDtrackCuts.h"
55#include "AliPhysicsSelection.h"
56#include "AliCentralitySelectionTask.h"
57#include "AliESDCaloCluster.h"
58#include "AliAODCaloCluster.h"
59#include "AliESDCaloTrigger.h"
60#include "AliGeomManager.h"
61#include "stdio.h"
62#include "TGeoManager.h"
63#include "iostream"
64#include "fstream"
65#include "AliEMCALTrack.h"
66//#include "AliEMCALTracker.h"
67#include "AliMagF.h"
68#include "AliKFParticle.h"
69#include "AliKFVertex.h"
70#include "AliHFEcontainer.h"
71#include "AliHFEcuts.h"
72#include "AliHFEpid.h"
73#include "AliHFEpidBase.h"
74#include "AliHFEpidQAmanager.h"
75#include "AliHFEtools.h"
76#include "AliCFContainer.h"
77#include "AliCFManager.h"
78#include "AliKFParticle.h"
79#include "AliKFVertex.h"
80#include "AliCentrality.h"
81#include "AliVEvent.h"
82#include "AliStack.h"
83#include "AliMCEvent.h"
84#include "TProfile.h"
85#include "AliFlowCandidateTrack.h"
86#include "AliFlowTrackCuts.h"
87#include "AliFlowEventSimple.h"
88#include "AliFlowCommonConstants.h"
89#include "AliFlowEvent.h"
90#include "TVector3.h"
91#include "TRandom2.h"
92#include "AliESDVZERO.h"
93#include "AliAODVZERO.h"
94#include "AliPID.h"
95#include "AliPIDResponse.h"
96#include "AliFlowTrack.h"
97#include "AliAnalysisTaskVnV0.h"
98#include "AliSelectNonHFE.h"
99
100
101class AliFlowTrackCuts;
102
103using namespace std;
104
105ClassImp(AliAnalysisTaskFlowITSTPCTOFQCSP)
106//________________________________________________________________________
107AliAnalysisTaskFlowITSTPCTOFQCSP::AliAnalysisTaskFlowITSTPCTOFQCSP(const char *name)
108: AliAnalysisTaskSE(name)
109,fDebug(0)
110,fAOD(0)
111,fOutputList(0)
112,fCuts(0)
113,fIdentifiedAsOutInz(kFALSE)
114,fPassTheEventCut(kFALSE)
115,fCFM(0)
116,fPID(0)
117,fPIDqa(0)
118,fCutsRP(0) // track cuts for reference particles
119,fNullCuts(0) // dummy cuts for flow event tracks
120,fFlowEvent(0) //! flow events (one for each inv mass band)
121,fkCentralityMethod(0)
122,fCentrality(0)
123,fCentralityMin(0)
124,fCentralityMax(0)
125,fInvmassCut(0)
126,fpTCutmin(0)
127,fpTCutmax(5)
128,fTrigger(0)
129,fPhi(0)
130,fEta(0)
131,fVZEROA(0)
132,fVZEROC(0)
133,fTPCM(0)
134,fNoEvents(0)
135,fInclusiveElecPt(0)
136,fTPCnsigma(0)
137,fTPCnsigmaAft(0)
138,fITSnsigmaAft(0)
139,fTPCnsigmaVSptAft(0)
140,fTOFns(0)
141,fTOFnsAft(0)
142,fTOFBetaAft(0)
143,fCentralityPass(0)
144,fCentralityNoPass(0)
145,fInvmassLS1(0)
146,fInvmassULS1(0)
147,fPhotoElecPt(0)
148,fSemiInclElecPt(0)
149,fULSElecPt(0)
150,fLSElecPt(0)
151,fMultCorAfterCuts(0)
152,fMultvsCentr(0)
153,fSubEventDPhiv2(0)
154,EPVzA(0)
155,EPVzC(0)
156,EPTPC(0)
157,fV2Phi(0)
158,fvertex(0)
159,fMultCorBeforeCuts(0)
160,fQAPid(0)
161,fminITSnsigmaLowpT(-1)
162,fmaxITSnsigmaLowpT(1)
163,fminITSnsigmaHighpT(-2)
164,fmaxITSnsigmaHighpT(2)
165,fminTPCnsigmaLowpT(-1)
166,fmaxTPCnsigmaLowpT(3)
167,fminTPCnsigmaHighpT(0)
168,fmaxTPCnsigmaHighpT(3)
169//,fQAPIDSparse(kFALSE)
170,fminTOFnSigma(-2)
171,fmaxTOFnSigma(2)
172,fQAPidSparse(0)
173,fTPCS(0)
174,fVz(0)
175,fOpeningAngleCut(0.1)
176,fOP_angle(0)
177,fOpeningAngleLS(0)
178,fOpeningAngleULS(0)
179,fNonHFE(new AliSelectNonHFE)
180,fDCA(0)
181,fITSnsigma(0)
182,fTPCnsigmaAftITSTOF(0)
183,fTPCnsigmaAftTOF(0)
184,fITSnsigmaAftTOF(0)
185,fITSvsTOF(0)
186,fTPCvsITS(0)
187,fTPCvsITSafterTOF(0)
188,fTPCvsTOF(0)
189,fMultCut(0)
190,fMultCorAfterCentrBeforeCuts(0)
191,fMultCorAfterVZTRKComp(0)
192,fCentralityBeforePileup(0)
193,fCentralityAfterVZTRK(0)
194,fCentralityAfterCorrCut(0)
195,fMultCorAfterCorrCut(0)
196,EPVz(0)
197,EPTPCp(0)
198,EPTPCn(0)
199,fSubEventDPhiv2new(0)
200,fV2Phivzerotot(0)
201,fHistCentrDistr(0x0)
202,fCentralityNoPassForFlattening(0)
164d5623 203,fptminAsso(0)
7fcd6b0f 204,fSparsephipsiULS(0)
205,fSparsephipsiLS(0)
206,fSparseMassULS(0)
207,fSparseMassLS(0)
208,fAssoTPCCluster(0)
209,fAssoITSRefit(0)
4c12b1ec
R		 210{
211 //Named constructor
212
213 fPID = new AliHFEpid("hfePid");
214 // Define input and output slots here
215 // Input slot #0 works with a TChain
216 DefineInput(0, TChain::Class());
217 // Output slot #0 id reserved by the base class for AOD
218 // Output slot #1 writes into a TH1 container
219 // DefineOutput(1, TH1I::Class());
220 DefineOutput(1, TList::Class());
221 DefineOutput(2, AliFlowEventSimple::Class());
222}
223
224//________________________________________________________________________
225AliAnalysisTaskFlowITSTPCTOFQCSP::AliAnalysisTaskFlowITSTPCTOFQCSP()
226: AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisElectFlow")
227,fDebug(0)
228,fAOD(0)
229,fOutputList(0)
230,fCuts(0)
231,fIdentifiedAsOutInz(kFALSE)
232,fPassTheEventCut(kFALSE)
233,fCFM(0)
234,fPID(0)
235,fPIDqa(0)
236,fCutsRP(0) // track cuts for reference particles
237,fNullCuts(0) // dummy cuts for flow event tracks
238,fFlowEvent(0) //! flow events (one for each inv mass band)
239,fkCentralityMethod(0)
240,fCentrality(0)
241,fCentralityMin(0)
242,fCentralityMax(0)
243,fInvmassCut(0)
244,fpTCutmin(0)
245,fpTCutmax(5)
246,fTrigger(0)
247,fPhi(0)
248,fEta(0)
249,fVZEROA(0)
250,fVZEROC(0)
251,fTPCM(0)
252,fNoEvents(0)
253,fInclusiveElecPt(0)
254,fTPCnsigma(0)
255,fTPCnsigmaAft(0)
256,fITSnsigmaAft(0)
257,fTPCnsigmaVSptAft(0)
258,fTOFns(0)
259,fTOFnsAft(0)
260,fTOFBetaAft(0)
261,fCentralityPass(0)
262,fCentralityNoPass(0)
263,fInvmassLS1(0)
264,fInvmassULS1(0)
265,fPhotoElecPt(0)
266,fSemiInclElecPt(0)
267,fULSElecPt(0)
268,fLSElecPt(0)
269,fMultCorAfterCuts(0)
270,fMultvsCentr(0)
271,fSubEventDPhiv2(0)
272,EPVzA(0)
273,EPVzC(0)
274,EPTPC(0)
275,fV2Phi(0)
276,fvertex(0)
277,fMultCorBeforeCuts(0)
278,fQAPid(0)
279,fminITSnsigmaLowpT(-1)
280,fmaxITSnsigmaLowpT(1)
281,fminITSnsigmaHighpT(-2)
282,fmaxITSnsigmaHighpT(2)
283,fminTPCnsigmaLowpT(-1)
284,fmaxTPCnsigmaLowpT(3)
285,fminTPCnsigmaHighpT(0)
286,fmaxTPCnsigmaHighpT(3)
287//,fQAPIDSparse(kFALSE)
288,fminTOFnSigma(-2)
289,fmaxTOFnSigma(2)
290,fQAPidSparse(0)
291,fTPCS(0)
292,fVz(0)
293,fOpeningAngleCut(0.1)
294,fOP_angle(0)
295,fOpeningAngleLS(0)
296,fOpeningAngleULS(0)
297,fNonHFE(new AliSelectNonHFE)
298,fDCA(0)
299,fITSnsigma(0)
300,fTPCnsigmaAftITSTOF(0)
301,fTPCnsigmaAftTOF(0)
302,fITSnsigmaAftTOF(0)
303,fITSvsTOF(0)
304,fTPCvsITS(0)
305,fTPCvsITSafterTOF(0)
306,fTPCvsTOF(0)
307,fMultCut(0)
308,fMultCorAfterCentrBeforeCuts(0)
309,fMultCorAfterVZTRKComp(0)
310,fCentralityBeforePileup(0)
311,fCentralityAfterVZTRK(0)
312,fCentralityAfterCorrCut(0)
313,fMultCorAfterCorrCut(0)
314,EPVz(0)
315,EPTPCp(0)
316,EPTPCn(0)
317,fSubEventDPhiv2new(0)
318,fV2Phivzerotot(0)
319,fHistCentrDistr(0x0)
320,fCentralityNoPassForFlattening(0)
164d5623 321,fptminAsso(0)
7fcd6b0f 322,fSparsephipsiULS(0)
323,fSparsephipsiLS(0)
324,fSparseMassULS(0)
325,fSparseMassLS(0)
326,fAssoTPCCluster(0)
327,fAssoITSRefit(0)
4c12b1ec
R		 328{
329 //Default constructor
330 fPID = new AliHFEpid("hfePid");
331 // Constructor
332 // Define input and output slots here
333 // Input slot #0 works with a TChain
334 DefineInput(0, TChain::Class());
335 // Output slot #0 id reserved by the base class for AOD
336 // Output slot #1 writes into a TH1 container
337 // DefineOutput(1, TH1I::Class());
338 DefineOutput(1, TList::Class());
339 DefineOutput(2, AliFlowEventSimple::Class());
340 //DefineOutput(3, TTree::Class());
341}
342//_________________________________________
343
344AliAnalysisTaskFlowITSTPCTOFQCSP::~AliAnalysisTaskFlowITSTPCTOFQCSP()
345{
346 //Destructor
347
348 delete fOutputList;
349 delete fPID;
350 // delete fPIDResponse;
351 delete fCFM;
352 delete fPIDqa;
353 if (fOutputList) delete fOutputList;
354 if (fFlowEvent) delete fFlowEvent;
355 delete fNonHFE;
356}
357//_________________________________________
358
359void AliAnalysisTaskFlowITSTPCTOFQCSP::UserExec(Option_t*)
360{
361 //Main loop
362 //Called for each event
363
364 // create pointer to event
365
366 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
367
368
369 if (!fAOD)
370 {
371 printf("ERROR: fAOD not available\n");
372 return;
373 }
374
375 if(!fCuts)
376 {
377 AliError("HFE cuts not available");
378 return;
379 }
380
381 if(!fPID->IsInitialized())
382 {
383 // Initialize PID with the given run number
384 AliWarning("PID not initialised, get from Run no");
385 fPID->InitializePID(fAOD->GetRunNumber());
386 }
387
388 // cout << "kTrigger == " << fTrigger <<endl;
389
390 if(fTrigger==0){
391 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral)) return;
392 }
393 if(fTrigger==1){
394 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral)) return;
395 }
396 if(fTrigger==2){
397 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMCEGA)) return;
398 }
399 if(fTrigger==3){
400 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB)) return;
401 }
402 if(fTrigger==4){
403 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kCentral | AliVEvent::kSemiCentral))) return;
404 }
405
406
407
408 //---------------CENTRALITY AND EVENT SELECTION-----------------------
409
410
411
412 Int_t fNOtrks = fAOD->GetNumberOfTracks();
413 Float_t vtxz = -999;
414 const AliAODVertex* trkVtx = fAOD->GetPrimaryVertex();
415 if (!trkVtx || trkVtx->GetNContributors()<=0)return;
416 TString vtxTtl = trkVtx->GetTitle();
417 if (!vtxTtl.Contains("VertexerTracks"))return;
418 const AliAODVertex* spdVtx = fAOD->GetPrimaryVertexSPD();
419 if (!spdVtx || spdVtx->GetNContributors()<=0)return;
420 if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5)return;
421 vtxz = trkVtx->GetZ();
422 if(TMath::Abs(vtxz)>fVz)return;
423
424 // Event cut
425 if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fAOD)) return;
426 if(fNOtrks<2) return;
427
428
429 Bool_t pass = kFALSE; //to select centrality
430 CheckCentrality(fAOD,pass);
431 if(!pass)return;
432
433 fvertex->Fill(vtxz);
434
435
436 fNoEvents->Fill(0);
437 PlotVZeroMultiplcities(fAOD);
438
439 SetNullCuts(fAOD);
440 PrepareFlowEvent(fAOD->GetNumberOfTracks(),fFlowEvent); //Calculate event plane Qvector and EP resolution for inclusive
441
442 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
443 if(!pidResponse)
444 {
445 AliDebug(1, "Using default PID Response");
446 pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
447 }
448
449 fPID->SetPIDResponse(pidResponse);
450
451 fCFM->SetRecEventInfo(fAOD);
452
453 // Look for kink mother
454 Int_t numberofvertices = fAOD->GetNumberOfVertices();
455 Double_t listofmotherkink[numberofvertices];
456 Int_t numberofmotherkink = 0;
457 for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
458 AliAODVertex *aodvertex = fAOD->GetVertex(ivertex);
459 if(!aodvertex) continue;
460 if(aodvertex->GetType()==AliAODVertex::kKink) {
461 AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
462 if(!mother) continue;
463 Int_t idmother = mother->GetID();
464 listofmotherkink[numberofmotherkink] = idmother;
465 //printf("ID %d\n",idmother);
466 numberofmotherkink++;
467 }
468 }
469
470 //=============================================V0EP from Alex======================================================================
471 Double_t qxEPa = 0, qyEPa = 0;
472 Double_t qxEPc = 0, qyEPc = 0;
473 Double_t qxEP = 0, qyEP = 0;
474
475 Double_t evPlAngV0A = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 8, 2, qxEPa, qyEPa);
476 Double_t evPlAngV0C = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 9, 2, qxEPc, qyEPc);
477 Double_t evPlAngV0 = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 10, 2, qxEP, qyEP);
478
479
480 Double_t Qx2 = 0, Qy2 = 0;
481 Double_t Qx2p = 0, Qy2p = 0;
482 Double_t Qx2n = 0, Qy2n = 0;
483
484 for (Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++){
485
486 AliAODTrack* aodTrack = fAOD->GetTrack(iT);
487
488 if (!aodTrack)
489 continue;
490
491 if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < 70) || (aodTrack->Pt() >= 20.0))
492 continue;
493
494 if (!aodTrack->TestFilterBit(128))
495 continue;
496
497
498 if(aodTrack->Eta()>0 && aodTrack->Eta()<0.8){
499
500 Qx2p += TMath::Cos(2*aodTrack->Phi());
501 Qy2p += TMath::Sin(2*aodTrack->Phi());
502 }
503 if(aodTrack->Eta()<0 && aodTrack->Eta()> -0.8){
504
505 Qx2n += TMath::Cos(2*aodTrack->Phi());
506 Qy2n += TMath::Sin(2*aodTrack->Phi());
507 }
508
509
510 Qx2 += TMath::Cos(2*aodTrack->Phi());
511 Qy2 += TMath::Sin(2*aodTrack->Phi());
512
513
514
515
516 }
517
518 Double_t evPlAngTPC = TMath::ATan2(Qy2, Qx2)/2.;
519 Double_t evPlAngTPCn = TMath::ATan2(Qy2n, Qx2n)/2.;
520 Double_t evPlAngTPCp = TMath::ATan2(Qy2p, Qx2p)/2.;
521
522 EPVzA->Fill(evPlAngV0A);
523 EPVzC->Fill(evPlAngV0C);
524 EPTPC->Fill(evPlAngTPC);
525
526 EPTPCn->Fill(evPlAngTPCn);
527 EPTPCp->Fill(evPlAngTPCp);
528 EPVz->Fill(evPlAngV0);
529
530 fSubEventDPhiv2->Fill(0.5, TMath::Cos(2.*(evPlAngV0A-evPlAngTPC))); // vzeroa - tpc
531 fSubEventDPhiv2->Fill(1.5, TMath::Cos(2.*(evPlAngV0A-evPlAngV0C))); // vzeroa - vzeroc
532 fSubEventDPhiv2->Fill(2.5, TMath::Cos(2.*(evPlAngV0C-evPlAngTPC))); // tpc - vzeroc
533
534
535 fSubEventDPhiv2new->Fill(0.5, TMath::Cos(2.*(evPlAngV0-evPlAngTPCp))); // vzero - tpcp
536 fSubEventDPhiv2new->Fill(1.5, TMath::Cos(2.*(evPlAngV0-evPlAngTPCn))); // vzero - tpcn
537 fSubEventDPhiv2new->Fill(2.5, TMath::Cos(2.*(evPlAngTPCp-evPlAngTPCn))); // tpcp - tpcn
538
539
540 //====================================================================================================================
541
542 AliAODTrack *track = NULL;
543
544 // Track loop
545 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++)
546 {
547 track = fAOD->GetTrack(iTracks);
548 if (!track)
549 {
550 printf("ERROR: Could not receive track %d\n", iTracks);
551 continue;
552 }
553 if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue; // TESTBIT FOR AOD double Counting
554
555 //--------------------------------------hfe begin-----------------------------------------------------------
556 //==========================================================================================================
557 //======================================track cuts==========================================================
558 if(track->Eta()<-0.8 || track->Eta()>0.8) continue; //eta cuts on candidates
559
560 // RecKine: ITSTPC cuts
561 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
562
563 // Reject kink mother
564 Bool_t kinkmotherpass = kTRUE;
565 for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
566 if(track->GetID() == listofmotherkink[kinkmother]) {
567 kinkmotherpass = kFALSE;
568 continue;
569 }
570 }
571 if(!kinkmotherpass) continue;
572
573 // RecPrim
574 // if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue; //deleted for DCA absence
575 // HFEcuts: ITS layers cuts
576 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
577 // HFE cuts: TPC PID cleanup
578 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
579 //==========================================================================================================
580 Double_t eta = track->Eta();
581 Double_t phi = track->Phi();
582 Double_t pt = track->Pt(); //pt track after cuts
583 if(pt<fpTCutmin || pt>fpTCutmax) continue;
584 //==========================================================================================================
585 //=========================================PID==============================================================
586 if(track->GetTPCsignalN() < fTPCS) continue;
587 Float_t fITSnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasITS(track, AliPID::kElectron) : 1000;
588 Float_t fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
589 Float_t fTOFnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTOF(track, AliPID::kElectron) : 1000;
590 // Float_t eDEDX = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kTRUE);
591 fITSnsigma->Fill(track->P(),fITSnSigma);
592 fTPCnsigma->Fill(track->P(),fTPCnSigma);
593 fTOFns->Fill(track->P(),fTOFnSigma);
594 fITSvsTOF->Fill(fTOFnSigma,fITSnSigma);
595 fTPCvsITS->Fill(fTPCnSigma,fITSnSigma);
596 fTPCvsTOF->Fill(fTPCnSigma,fTOFnSigma);
597
598 if( pt >= 0.3){
599 if(fTOFnSigma < fminTOFnSigma || fTOFnSigma > fmaxTOFnSigma) continue;
600 }//cuts on nsigma tof full pt range
601
602 fITSnsigmaAftTOF->Fill(track->P(),fITSnSigma);
603 fTPCnsigmaAftTOF->Fill(track->P(),fTPCnSigma);
604 fTPCvsITSafterTOF->Fill(fTPCnSigma,fITSnSigma);
605
606 Double_t valPidSparse[3] = {
607 pt,
608 fITSnSigma,
609 fTPCnSigma,
610 };
611 fQAPidSparse->Fill(valPidSparse);
612
613
614 if( pt < 1.5){
615 if(fITSnSigma < fminITSnsigmaLowpT || fITSnSigma > fmaxITSnsigmaLowpT)continue;
616 }//cuts on nsigma its low pt
617 if( pt >= 1.5){
618 if(fITSnSigma < fminITSnsigmaHighpT || fITSnSigma > fmaxITSnsigmaHighpT)continue;
619 }//cuts on nsigma its high pt
620 fTPCnsigmaAftITSTOF->Fill(track->P(),fTPCnSigma);
621 if(pt >= 0.25 && pt < 1.5){
622 if(fTPCnSigma < fminTPCnsigmaLowpT || fTPCnSigma > fmaxTPCnsigmaLowpT) continue;
623 }//cuts on nsigma tpc lowpt
624 if( pt >= 1.5){
625 if(fTPCnSigma < fminTPCnsigmaHighpT || fTPCnSigma > fmaxTPCnsigmaHighpT) continue;
626 }//cuts on nsigma tpc high pt
627 fTPCnsigmaAft->Fill(track->P(),fTPCnSigma);
628 fTPCnsigmaVSptAft->Fill(pt,fTPCnSigma);
629
630 //==========================================================================================================
631 //=========================================QA PID SPARSE====================================================
632 Float_t timeTOF = track->GetTOFsignal();
633 Double_t intTime[5] = {-99., -99., -99., -99., -99.};
634 track->GetIntegratedTimes(intTime);
635 Float_t timeElec = intTime[0];
636 Float_t intLength = 2.99792458e-2* timeElec;
637 Double_t beta = 0.1;
638 if ((intLength > 0) && (timeTOF > 0))
639 beta = intLength/2.99792458e-2/timeTOF;
640
641 // if(fQAPIDSparse){
642 // Double_t valPid[4] = {
643 // track->P(),
644 // track->GetTPCsignal(),
645 // beta,
646 // track->Charge()
647 // };
648 // fQAPid->Fill(valPid);
649 // }
650
651
652 fITSnsigmaAft->Fill(track->P(),fITSnSigma);
653 fTPCnsigmaAft->Fill(track->P(),fTPCnSigma);
654 fTOFnsAft->Fill(track->P(),fTOFnSigma);
655 fTOFBetaAft->Fill(track->P(),beta);
656 fInclusiveElecPt->Fill(pt);
657 fPhi->Fill(phi);
658 fEta->Fill(eta);
659 //=========================================================================================================
660 //----------------------Flow of Inclusive Electrons--------------------------------------------------------
661 AliFlowTrack *sTrack = new AliFlowTrack();
662 sTrack->Set(track);
663 sTrack->SetID(track->GetID());
664 sTrack->SetForRPSelection(kTRUE);
665 sTrack->SetForPOISelection(kTRUE);
666 sTrack->SetMass(263732);
667 for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs)
668 {
669 // cout << " no of rps " << iRPs << endl;
670 AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
671 if (!iRP) continue;
672 if (!iRP->InRPSelection()) continue;
673 if( sTrack->GetID() == iRP->GetID())
674 {
675 if(fDebug) printf(" was in RP set");
676 // cout << sTrack->GetID() <<" == " << iRP->GetID() << " was in RP set====REMOVED" <<endl;
677 iRP->SetForRPSelection(kFALSE);
678 // fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
679 }
680 } //end of for loop on RPs
681 fFlowEvent->InsertTrack(((AliFlowTrack*) sTrack));
682 fFlowEvent->SetNumberOfPOIs(fFlowEvent->GetNumberOfPOIs()+1);
683 //============================Event Plane Method with V0====================================================
684 Double_t v2PhiV0A = TMath::Cos(2*(phi - evPlAngV0A));
685 Double_t v2PhiV0C = TMath::Cos(2*(phi - evPlAngV0C));
686 Double_t v2Phi[3] = {
687 v2PhiV0A,
688 v2PhiV0C,
689 pt};
690 fV2Phi->Fill(v2Phi);
691
692 Double_t v2PhiVz = TMath::Cos(2*(phi - evPlAngV0));
693 Double_t v2PhiV0tot[2] = {
694 v2PhiVz,
695 pt};
696 fV2Phivzerotot->Fill(v2PhiV0tot);
697
698 //==========================================================================================================
699 //=========================================================================================================
700
701
702 if(fDCA){
703 fNonHFE = new AliSelectNonHFE();
704 fNonHFE->SetAODanalysis(kTRUE);
705 fNonHFE->SetInvariantMassCut(fInvmassCut);
706 if(fOP_angle) fNonHFE->SetOpeningAngleCut(fOpeningAngleCut);
707 //fNonHFE->SetChi2OverNDFCut(fChi2Cut);
708 //if(fDCAcutFlag) fNonHFE->SetDCACut(fDCAcut);
709 fNonHFE->SetAlgorithm("DCA"); //KF
710 fNonHFE->SetPIDresponse(pidResponse);
711 fNonHFE->SetTrackCuts(-3,3);
712
713 fNonHFE->SetHistAngleBack(fOpeningAngleLS);
714 fNonHFE->SetHistAngle(fOpeningAngleULS);
715 //fNonHFE->SetHistDCABack(fDCABack);
716 //fNonHFE->SetHistDCA(fDCA);
717 fNonHFE->SetHistMassBack(fInvmassLS1);
718 fNonHFE->SetHistMass(fInvmassULS1);
719
720 fNonHFE->FindNonHFE(iTracks,track,fAOD);
721
722 // Int_t *fUlsPartner = fNonHFE->GetPartnersULS();
723 // Int_t *fLsPartner = fNonHFE->GetPartnersLS();
724 // Bool_t fUlsIsPartner = kFALSE;
725 // Bool_t fLsIsPartner = kFALSE;
726 if(fNonHFE->IsULS()){
727 for(Int_t kULS =0; kULS < fNonHFE->GetNULS(); kULS++){
728 fULSElecPt->Fill(track->Pt());
729 }
730 }
731
732 if(fNonHFE->IsLS()){
733 for(Int_t kLS =0; kLS < fNonHFE->GetNLS(); kLS++){
734 fLSElecPt->Fill(track->Pt());
735 }
736 }
737 }
738 if(!fDCA){
739 //=========================================================================================================
740 //----------------------Selection and Flow of Photonic Electrons-----------------------------
741 Bool_t fFlagPhotonicElec = kFALSE;
7fcd6b0f 742 SelectPhotonicElectron(iTracks,track,fTPCnSigma,evPlAngV0,fFlagPhotonicElec);
4c12b1ec
R		 743 if(fFlagPhotonicElec){fPhotoElecPt->Fill(pt);}
744 // Semi inclusive electron
745 if(!fFlagPhotonicElec){fSemiInclElecPt->Fill(pt);}
746 }
747 //-------------------------------------------------------------------------------------------
748
749 }//end loop on track
750 PostData(1, fOutputList);
751 PostData(2, fFlowEvent);
752
753 //----------hfe end---------
754}
755//_________________________________________
7fcd6b0f 756void AliAnalysisTaskFlowITSTPCTOFQCSP::SelectPhotonicElectron(Int_t itrack,const AliAODTrack *track,Float_t fTPCnSigma,Double_t evPlAngV0, Bool_t &fFlagPhotonicElec)
4c12b1ec
R		 757{
758
759 //Identify non-heavy flavour electrons using Invariant mass method
760 Bool_t flagPhotonicElec = kFALSE;
761
762 for(Int_t jTracks = 0; jTracks<fAOD->GetNumberOfTracks(); jTracks++){
763 AliAODTrack *trackAsso = fAOD->GetTrack(jTracks);
764 if (!trackAsso) {
765 printf("ERROR: Could not receive track %d\n", jTracks);
766 continue;
767 }
768 // if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue; // TESTBIT FOR AOD double Counting
769 if(!trackAsso->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
7fcd6b0f 770
771 if(fAssoITSRefit){
772 if(!(trackAsso->GetStatus()&AliESDtrack::kITSrefit)) continue;
773 }
774
775 if(!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)) continue;
776
777// if((!(trackAsso->GetStatus()&AliESDtrack::kITSrefit)|| (!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
4c12b1ec
R		 778
779
780 if(jTracks == itrack) continue;
781 Double_t ptAsso=-999., nsigma=-999.0;
782 Double_t mass=-999., width = -999;
783 Double_t openingAngle = -999.;
784 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
785
786
787 ptAsso = trackAsso->Pt();
788 Short_t chargeAsso = trackAsso->Charge();
789 Short_t charge = track->Charge();
790 // nsigma = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
791 nsigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
792
793
7fcd6b0f 794 //if(trackAsso->GetTPCNcls() < 80) continue;
795 if(trackAsso->GetTPCNcls() < fAssoTPCCluster) continue;
4c12b1ec
R		 796 if(nsigma < -3 || nsigma > 3) continue;
797 if(trackAsso->Eta()<-0.9 || trackAsso->Eta()>0.9) continue;
164d5623 798 if(ptAsso < fptminAsso) continue;
4c12b1ec
R		 799
800 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
801 if(charge>0) fPDGe1 = -11;
802 if(chargeAsso>0) fPDGe2 = -11;
803
804 if(charge == chargeAsso) fFlagLS = kTRUE;
805 if(charge != chargeAsso) fFlagULS = kTRUE;
806
807 AliKFParticle::SetField(fAOD->GetMagneticField());
808 AliKFParticle ge1 = AliKFParticle(*track, fPDGe1);
809 AliKFParticle ge2 = AliKFParticle(*trackAsso, fPDGe2);
810 AliKFParticle recg(ge1, ge2);
811
812 if(recg.GetNDF()<1) continue;
813 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
814 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
815 recg.GetMass(mass,width);
816
817 openingAngle = ge1.GetAngle(ge2);
818 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
819 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
820 if(fOP_angle){if(openingAngle > fOpeningAngleCut) continue;}
821
822
823 if(fFlagLS) fInvmassLS1->Fill(mass);
824 if(fFlagULS) fInvmassULS1->Fill(mass);
825
7fcd6b0f 826 if(fFlagULS){
827 Double_t MassSparseULS[3] = {
828 track->Pt(),
829 mass
830 };
831 fSparseMassULS->Fill(MassSparseULS);
832 }
833 if(fFlagLS){
834 Double_t MassSparseLS[3] = {
835 track->Pt(),
836 mass
837 };
838 fSparseMassLS->Fill(MassSparseLS);
839 }
840
841
4c12b1ec
R		 842 if(mass<fInvmassCut){
843 if(fFlagULS){fULSElecPt->Fill(track->Pt());}
844 if(fFlagLS){fLSElecPt->Fill(track->Pt());}
845 }
846
7fcd6b0f 847
848 Double_t phi = track->Phi();
849 Float_t DeltaPhi_eEP = TVector2::Phi_0_2pi(phi - evPlAngV0);
850 if(DeltaPhi_eEP > TMath::Pi()) {DeltaPhi_eEP = DeltaPhi_eEP - TMath::Pi();}
851
852
853 if(mass<fInvmassCut){
854 if(fFlagULS){
855 Double_t ulsSparse[3] = {
856 track->Pt(),
857 fTPCnSigma,
858 DeltaPhi_eEP
859 };
860 fSparsephipsiULS->Fill(ulsSparse);
861 }
862 if(fFlagLS){
863 Double_t lsSparse[3] = {
864 track->Pt(),
865 fTPCnSigma,
866 DeltaPhi_eEP
867 };
868 fSparsephipsiLS->Fill(lsSparse);
869 }
870 }
871
872
873
4c12b1ec
R		 874 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
875 flagPhotonicElec = kTRUE;
876 }
877 }//track loop
878
879 fFlagPhotonicElec = flagPhotonicElec;
880}
881//___________________________________________
882void AliAnalysisTaskFlowITSTPCTOFQCSP::UserCreateOutputObjects()
883{
884
885 //Create histograms
886 //----------hfe initialising begin---------
887 fNullCuts = new AliFlowTrackCuts("null_cuts");
888
889 AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
890 cc->SetNbinsMult(10000);
891 cc->SetMultMin(0);
892 cc->SetMultMax(10000);
893
894 cc->SetNbinsPt(100);
895 cc->SetPtMin(0);
896 cc->SetPtMax(5);
897
898 cc->SetNbinsPhi(180);
899 cc->SetPhiMin(0.0);
900 cc->SetPhiMax(TMath::TwoPi());
901
902 cc->SetNbinsEta(30);
903 cc->SetEtaMin(-8.0);
904 cc->SetEtaMax(+8.0);
905
906 cc->SetNbinsQ(500);
907 cc->SetQMin(0.0);
908 cc->SetQMax(3.0);
909
910
911 // AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
912 // AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
913 // if (!inputHandler){
914 // AliFatal("Input handler needed");
915 // }
916 // else{
917 // fPIDResponse=inputHandler->GetPIDResponse();
918 // }
919 //pid response object
920 // if (!fPIDResponse)AliError("PIDResponse object was not created");
921
922
923 //--------Initialize PID
924 fPID->SetHasMCData(kFALSE);
925 if(!fPID->GetNumberOfPIDdetectors())
926 {
927 fPID->AddDetector("ITS", 0);
928 fPID->AddDetector("TOF", 1);
929 fPID->AddDetector("TPC", 2);
930
931 }
932
933 fPID->SortDetectors();
934 fPIDqa = new AliHFEpidQAmanager();
935 fPIDqa->Initialize(fPID);
936
937
938
939 //--------Initialize correction Framework and Cuts
940 fCFM = new AliCFManager;
941 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
942 fCFM->SetNStepParticle(kNcutSteps);
943 for(Int_t istep = 0; istep < kNcutSteps; istep++)
944 fCFM->SetParticleCutsList(istep, NULL);
945
946 if(!fCuts){
947 AliWarning("Cuts not available. Default cuts will be used");
948 fCuts = new AliHFEcuts;
949 fCuts->CreateStandardCuts();
950 }
951
952 fCuts->SetAOD();
953 fCuts->Initialize(fCFM);
954 //----------hfe initialising end--------
955 //---------Output Tlist
956 fOutputList = new TList();
957 fOutputList->SetOwner();
958 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
959
960 fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
961 fOutputList->Add(fNoEvents);
962
963 fITSnsigma = new TH2F("fITSnsigma", "ITS - n sigma before HFE pid",600,0,6,400,-20,20);
964 fOutputList->Add(fITSnsigma);
965
966 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma before HFE pid",600,0,6,400,-20,20);
967 fOutputList->Add(fTPCnsigma);
968
969 fITSnsigmaAft = new TH2F("fITSnsigmaAft", "ITS - n sigma after HFE pid",1000,0,10,300,-10,20);
970 fOutputList->Add(fITSnsigmaAft);
971 fITSvsTOF = new TH2F("fITSvsTOF", "ITS tof",400,-20,20,400,-20,20);
972 fOutputList->Add(fITSvsTOF);
973 fTPCvsITS = new TH2F("TPCvsITS", "TPC ITS",400,-20,20,400,-20,20);
974 fOutputList->Add(fTPCvsITS);
975 fTPCvsTOF = new TH2F("TPCvsTOF", "TPC TOF",400,-20,20,400,-20,20);
976 fOutputList->Add(fTPCvsTOF);
977 fTPCvsITSafterTOF = new TH2F("TPCvsITSafterTOF", "TPC ITS",400,-20,20,400,-20,20);
978 fOutputList->Add(fTPCvsITSafterTOF);
979
980
981 fITSnsigmaAftTOF = new TH2F("fITSnsigmaAftTOF", "ITS - n sigma after HFE pid",600,0,6,400,-20,20);
982 fOutputList->Add(fITSnsigmaAftTOF);
983
984 fTPCnsigmaAft = new TH2F("fTPCnsigmaAft", "TPC - n sigma after HFE pid",600,0,6,400,-20,20);
985 fOutputList->Add(fTPCnsigmaAft);
986
987 fTPCnsigmaVSptAft = new TH2F("fTPCnsigmaVSptAft", "TPC - n sigma after HFE pid",600,0,6,400,-20,20);
988 fOutputList->Add(fTPCnsigmaVSptAft);
989
990 fTPCnsigmaAftITSTOF = new TH2F("fTPCnsigmaAftITSTOF", "TPC - n sigma after HFE pid",600,0,6,400,-20,20);
991 fOutputList->Add(fTPCnsigmaAftITSTOF);
992
993 fTPCnsigmaAftTOF = new TH2F("fTPCnsigmaAftTOF", "TPC - n sigma after HFE pid",600,0,6,400,-20,20);
994 fOutputList->Add(fTPCnsigmaAftTOF);
995
996 fTOFns = new TH2F("fTOFns","track TOFnSigma",600,0,6,400,-20,20);
997 fOutputList->Add(fTOFns);
998
999 fTOFnsAft = new TH2F("fTOFnsAft","track TOFnSigma",600,0,6,400,-20,20);
1000 fOutputList->Add(fTOFnsAft);
1001
1002 fTOFBetaAft = new TH2F("fTOFBetaAft","track TOFBeta",600,0,6,120,0,1.2);
1003 fOutputList->Add(fTOFBetaAft);
1004
1005 fInclusiveElecPt = new TH1F("fInclElecPt", "Inclusive electron pt",100,0,5);
1006 fOutputList->Add(fInclusiveElecPt);
1007
1008 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,5);
1009 fOutputList->Add(fPhotoElecPt);
1010
1011 fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",100,0,5);
1012 fOutputList->Add(fSemiInclElecPt);
1013
1014 fULSElecPt = new TH1F("fULSElecPt", "ULS electron pt",100,0,5);
1015 fOutputList->Add(fULSElecPt);
1016
1017 fLSElecPt = new TH1F("fLSElecPt", "LS electron pt",100,0,5);
1018 fOutputList->Add(fLSElecPt);
1019
1020 fInvmassLS1 = new TH1F("fInvmassLS1", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
1021 fOutputList->Add(fInvmassLS1);
1022
1023 fInvmassULS1 = new TH1F("fInvmassULS1", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
1024 fOutputList->Add(fInvmassULS1);
1025
1026 fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
1027 fOutputList->Add(fCentralityPass);
1028
1029 fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
1030 fOutputList->Add(fCentralityNoPass);
1031
1032 fCentralityNoPassForFlattening = new TH1F("fCentralityNoPassForFlattening", "Centrality No Pass for flattening", 101, -1, 100);
1033 fOutputList->Add(fCentralityNoPassForFlattening);
1034
1035 fCentralityBeforePileup = new TH1F("fCentralityBeforePileup", "fCentralityBeforePileup Pass", 101, -1, 100);
1036 fOutputList->Add(fCentralityBeforePileup);
1037
1038 fCentralityAfterVZTRK = new TH1F("fCentralityAfterVZTRK", "fCentralityAfterVZTRK Pass", 101, -1, 100);
1039 fOutputList->Add(fCentralityAfterVZTRK);
1040
1041 fCentralityAfterCorrCut = new TH1F("fCentralityAfterCorrCut", "fCentralityAfterCorrCut Pass", 101, -1, 100);
1042 fOutputList->Add(fCentralityAfterCorrCut);
1043
1044 fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
1045 fOutputList->Add(fPhi);
1046
1047 fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
1048 fOutputList->Add(fEta);
1049
1050 fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
1051 fOutputList->Add(fVZEROA);
1052
1053 fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
1054 fOutputList->Add(fVZEROC);
1055
1056 fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
1057 fOutputList->Add(fTPCM);
1058
1059 fvertex = new TH1D("fvertex", "vertex distribution", 300, -15,15);
1060 fOutputList->Add(fvertex);
1061
1062 fMultCorBeforeCuts = new TH2F("fMultCorBeforeCuts", "TPC vs Global multiplicity (Before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1063 fOutputList->Add(fMultCorBeforeCuts);
1064
1065 fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1066 fOutputList->Add(fMultCorAfterCuts);
1067
1068 fMultCorAfterCentrBeforeCuts = new TH2F("fMultCorAfterCentrBeforeCuts", "TPC vs Global multiplicity (After CC before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1069 fOutputList->Add(fMultCorAfterCentrBeforeCuts);
1070
1071 fMultCorAfterVZTRKComp = new TH2F("fMultCorAfterVZTRKComp", "TPC vs Global multiplicity (After V0-TRK); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1072 fOutputList->Add(fMultCorAfterVZTRKComp);
1073
1074 fMultCorAfterCorrCut = new TH2F("fMultCorAfterCorrCut", "TPC vs Global multiplicity (After CorrCut); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1075 fOutputList->Add(fMultCorAfterCorrCut);
1076
1077 fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 100, 0., 100, 100, 0, 3000);
1078 fOutputList->Add(fMultvsCentr);
1079
1080 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
1081 fOutputList->Add(fOpeningAngleLS);
1082
1083 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
1084 fOutputList->Add(fOpeningAngleULS);
1085
1086
1087
1088 //----------------------------------------------------------------------------
1089 EPVzA = new TH1D("EPVzA", "EPVzA", 80, -2, 2);
1090 fOutputList->Add(EPVzA);
1091 EPVzC = new TH1D("EPVzC", "EPVzC", 80, -2, 2);
1092 fOutputList->Add(EPVzC);
1093 EPTPC = new TH1D("EPTPC", "EPTPC", 80, -2, 2);
1094 fOutputList->Add(EPTPC);
1095 EPVz = new TH1D("EPVz", "EPVz", 80, -2, 2);
1096 fOutputList->Add(EPVz);
1097 EPTPCp = new TH1D("EPTPCp", "EPTPCp", 80, -2, 2);
1098 fOutputList->Add(EPTPCp);
1099 EPTPCn = new TH1D("EPTPCn", "EPTPCn", 80, -2, 2);
1100 fOutputList->Add(EPTPCn);
1101
1102
1103 //----------------------------------------------------------------------------
1104 fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 3, 0, 3);
1105 fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
1106 fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
1107 fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
1108 fOutputList->Add(fSubEventDPhiv2);
1109
1110 fSubEventDPhiv2new = new TProfile("fSubEventDPhiv2new", "fSubEventDPhiv2new", 3, 0, 3);
1111 fSubEventDPhiv2new->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
1112 fSubEventDPhiv2new->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
1113 fSubEventDPhiv2new->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
1114 fOutputList->Add(fSubEventDPhiv2new);
1115
1116 //================================Event Plane with VZERO=====================
1117 const Int_t nPtBins = 12;
1118 Double_t binsPt[nPtBins+1] = {0, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 2, 2.5, 3, 3.5, 4, 5};
1119 // v2A, v2C, pt
1120 Int_t bins[3] = { 50, 50, nPtBins};
1121 Double_t xmin[3] = { -1., -1., 0};
1122 Double_t xmax[3] = { 1., 1., 5};
1123 fV2Phi = new THnSparseF("fV2Phi", "v2A:v2C:pt", 3, bins, xmin, xmax);
1124 // Set bin limits for axes which are not standard binned
1125 fV2Phi->SetBinEdges(2, binsPt);
1126 // set axes titles
1127 fV2Phi->GetAxis(0)->SetTitle("v_{2} (V0A)");
1128 fV2Phi->GetAxis(1)->SetTitle("v_{2} (V0C)");
1129 fV2Phi->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
1130 fOutputList->Add(fV2Phi);
1131 //================================Event Plane with VZERO=====================
1132 Int_t binsV[2] = { 50, 100};
1133 Double_t xminV[2] = { -1., 0};
1134 Double_t xmaxV[2] = { 1., 5};
1135 fV2Phivzerotot = new THnSparseF("fV2Phivzerotot", "v2:pt", 2, binsV, xminV, xmaxV);
1136 // Set bin limits for axes which are not standard binned
1137 //fV2Phivzerotot->SetBinEdges(1, binsPt);
1138 // set axes titles
1139 fV2Phivzerotot->GetAxis(0)->SetTitle("v_{2} (V0)");
1140 fV2Phivzerotot->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
1141 fOutputList->Add(fV2Phivzerotot);
1142
1143
1144
1145 //----------------------------------------------------------------------------
1146 //----------------------------------------------------------------------------
1147 // if(fQAPIDSparse){
1148 // Int_t binsQA[4] = { 150, 100, 120, 3};
1149 // Double_t xminQA[4] = { 0., 50, 0, -1.5};
1150 // Double_t xmaxQA[4] = { 15., 150, 1.2, 1.5};
1151 // fQAPid = new THnSparseF("fQAPid", "p:dEdx:beta:ch", 4, binsQA, xminQA, xmaxQA);
1152 // fQAPid->GetAxis(0)->SetTitle("p (Gev/c");
1153 // fQAPid->GetAxis(1)->SetTitle("dE/dx");
1154// fQAPid->GetAxis(2)->SetTitle("#beta (TOF)");
1155// fQAPid->GetAxis(3)->SetTitle("charge");
1156// fOutputList->Add(fQAPid);
1157// }
1158 //===========================================================================
1159 Int_t binsQA2[3] = { 100, 40, 150/*, 60*/};
1160 Double_t xminQA2[3] = { 0., -2, -15/*, -3*/};
1161 Double_t xmaxQA2[3] = { 5., 2, 15/*, 3*/};
1162 fQAPidSparse = new THnSparseF("fQAPidSparse", "pt:itsnsigma:tpcnsigma", 3, binsQA2, xminQA2, xmaxQA2);
1163 fQAPidSparse->GetAxis(0)->SetTitle("pt (Gev/c)");
1164 fQAPidSparse->GetAxis(1)->SetTitle("itsnsigma");
1165 fQAPidSparse->GetAxis(2)->SetTitle("tpcnsigma");
1166 fOutputList->Add(fQAPidSparse);
1167 //===========================================================================
7fcd6b0f 1168
1169
1170
1171 Int_t binsphipsi[3] = { 100, 150, 6};
1172 Double_t xminphipsi[3] = { 0., -15, 0};
1173 Double_t xmaxphipsi[3] = { 5., 15, TMath::Pi()};
1174 fSparsephipsiULS = new THnSparseF("fSparsephipsiULS", "pt:tpcnsigma:DeltaPhiULS", 3, binsphipsi, xminphipsi, xmaxphipsi);
1175 fSparsephipsiULS->GetAxis(0)->SetTitle("pt (Gev/c)");
1176 fSparsephipsiULS->GetAxis(1)->SetTitle("tpcnsigma");
1177 fSparsephipsiULS->GetAxis(2)->SetTitle("DeltaPhiULS");
1178 fOutputList->Add(fSparsephipsiULS);
1179
1180 fSparsephipsiLS = new THnSparseF("fSparsephipsiLS", "pt:tpcnsigma:DeltaPhiLS", 3, binsphipsi, xminphipsi, xmaxphipsi);
1181 fSparsephipsiLS->GetAxis(0)->SetTitle("pt (Gev/c)");
1182 fSparsephipsiLS->GetAxis(1)->SetTitle("tpcnsigma");
1183 fSparsephipsiLS->GetAxis(2)->SetTitle("DeltaPhiLS");
1184 fOutputList->Add(fSparsephipsiLS);
1185
1186
1187 Int_t binsmass[2] = { 100, 200};
1188 Double_t xminmass[2] = { 0., 0};
1189 Double_t xmaxmass[2] = { 5., 1.};
1190 fSparseMassULS = new THnSparseF("fSparseMassULS", "pt:mass (GeV/c^{2})", 2, binsmass, xminmass, xmaxmass);
1191 fSparseMassULS->GetAxis(0)->SetTitle("pt (Gev/c)");
1192 fSparseMassULS->GetAxis(1)->SetTitle("mass");
1193 fOutputList->Add(fSparseMassULS);
1194
1195 fSparseMassLS = new THnSparseF("fSparseMassLS", "pt:mass (GeV/c^{2})", 2, binsmass, xminmass, xmaxmass);
1196 fSparseMassLS->GetAxis(0)->SetTitle("pt (Gev/c)");
1197 fSparseMassLS->GetAxis(1)->SetTitle("mass");
1198 fOutputList->Add(fSparseMassLS);
1199
1200
4c12b1ec
R		 1201 PostData(1,fOutputList);
1202 // create and post flowevent
1203 fFlowEvent = new AliFlowEvent(10000);
1204 PostData(2, fFlowEvent);
1205
1206}
1207//________________________________________________________________________
1208void AliAnalysisTaskFlowITSTPCTOFQCSP::Terminate(Option_t *)
1209{
1210 // Info("Terminate");
1211 AliAnalysisTaskSE::Terminate();
1212}
1213//_____________________________________________________________________________
1214template <typename T> void AliAnalysisTaskFlowITSTPCTOFQCSP::PlotVZeroMultiplcities(const T* event) const
1215{
1216 // QA multiplicity plots
1217 fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
1218 fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
1219}
1220//_____________________________________________________________________________
1221template <typename T> void AliAnalysisTaskFlowITSTPCTOFQCSP::SetNullCuts(T* event)
1222{
1223 //Set null cuts
1224 if (fDebug) cout << " fCutsRP " << fCutsRP << endl;
1225 fCutsRP->SetEvent(event, MCEvent());
1226 fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
1227 fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
1228 fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
1229 fNullCuts->SetEvent(event, MCEvent());
1230}
1231//_____________________________________________________________________________
1232void AliAnalysisTaskFlowITSTPCTOFQCSP::PrepareFlowEvent(Int_t iMulti, AliFlowEvent *FlowEv) const
1233{
1234 //Prepare flow events
1235 FlowEv->ClearFast();
1236 FlowEv->Fill(fCutsRP, fNullCuts);
1237 FlowEv->SetReferenceMultiplicity(iMulti);
1238 FlowEv->DefineDeadZone(0, 0, 0, 0);
1239 // FlowEv->TagSubeventsInEta(-0.7, 0, 0, 0.7);
1240}
1241//_____________________________________________________________________________
1242Bool_t AliAnalysisTaskFlowITSTPCTOFQCSP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1243{
1244 // Check single track cuts for a given cut step
1245 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1246 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1247 return kTRUE;
1248}
1249//_________________________________________
1250void AliAnalysisTaskFlowITSTPCTOFQCSP::CheckCentrality(AliAODEvent* event, Bool_t &centralitypass)
1251{
1252 //============================Multiplicity TPV vs Global===============================================================================
1253 const Int_t nGoodTracks = event->GetNumberOfTracks();
1254 Float_t multTPC(0.); // tpc mult estimate
1255 Float_t multGlob(0.); // global multiplicity
1256 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
1257 AliAODTrack* trackAOD = event->GetTrack(iTracks);
1258 if (!trackAOD) continue;
1259 if (!(trackAOD->TestFilterBit(1))) continue;
1260 if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.2)) continue;
1261 multTPC++;
1262 }
1263 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
1264 AliAODTrack* trackAOD = event->GetTrack(iTracks);
1265 if (!trackAOD) continue;
1266 if (!(trackAOD->TestFilterBit(16))) continue;
1267 if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.1)) continue;
1268 Double_t b[2] = {-99., -99.};
1269 Double_t bCov[3] = {-99., -99., -99.};
1270 if (!(trackAOD->PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
1271 if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
1272 multGlob++;
1273 } //track loop
1274 fMultCorBeforeCuts->Fill(multGlob, multTPC);//before all cuts...even before centrality selectrion
1275 //============================================================================================================================
1276 // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
1277 if (!fkCentralityMethod) AliFatal("No centrality method set! FATAL ERROR!");
1278 fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethod);
1279 // cout << "--------------Centrality evaluated-------------------------"<<endl;
1280 if ((fCentrality <= fCentralityMin) || (fCentrality > fCentralityMax))
1281 {
1282 fCentralityNoPass->Fill(fCentrality);
1283 // cout << "--------------Fill no pass-----"<< fCentrality <<"--------------------"<<endl;
1284 centralitypass = kFALSE;
1285 }else
1286 {
1287 // cout << "--------------Fill pass----"<< fCentrality <<"---------------------"<<endl;
1288 centralitypass = kTRUE;
1289 }
1290 if (centralitypass){
1291 fMultCorAfterCentrBeforeCuts->Fill(multGlob, multTPC);
1292 fCentralityBeforePileup->Fill(fCentrality);
1293 }//...after centrality selectrion
1294 //============================================================================================================================
1295 //to remove the bias introduced by multeplicity outliers---------------------
1296 Float_t centTrk = event->GetCentrality()->GetCentralityPercentile("TRK");
1297 Float_t centv0 = event->GetCentrality()->GetCentralityPercentile("V0M");
1298 if (TMath::Abs(centv0 - centTrk) > 5.0){
1299 centralitypass = kFALSE;
1300 fCentralityNoPass->Fill(fCentrality);
1301 }
1302 if (centralitypass){
1303 fMultCorAfterVZTRKComp->Fill(multGlob, multTPC);
1304 fCentralityAfterVZTRK->Fill(fCentrality);
1305 }//...after centrality selectrion
1306 //============================================================================================================================
1307 if(fMultCut){
1308 if(fTrigger==1 || fTrigger==4){
1309 if(! (multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))){
1310 // cout <<" Trigger ==" <<fTrigger<< endl;
1311 centralitypass = kFALSE;
1312 fCentralityNoPass->Fill(fCentrality);
1313 }//2011 Semicentral
1314 }
1315 if(fTrigger==0){
1316 if(! (multTPC > (77.9 + 1.395*multGlob) && multTPC < (187.3 + 1.665*multGlob))){
1317 // cout <<" Trigger ==" <<fTrigger<< endl;
1318 centralitypass = kFALSE;
1319 fCentralityNoPass->Fill(fCentrality);
1320 }//2011
1321 }//2011 Central
1322 }
1323 if (centralitypass){
1324 fMultCorAfterCorrCut->Fill(multGlob, multTPC);
1325 fCentralityAfterCorrCut->Fill(fCentrality);
1326 }//...after CORR CUT
1327 //=================================All cuts are passed==================++++==================================================
1328 //=================================Now Centrality flattening for central trigger==================++++==================================================
1329 if(fTrigger==0 || fTrigger==4){
1330 if(!IsEventSelectedForCentrFlattening(fCentrality)){
1331 centralitypass = kFALSE;
1332 fCentralityNoPassForFlattening->Fill(fCentrality);
1333 }
1334 }
1335 //==============================fill histo after all cuts==============================++++==================================================
1336 if(centralitypass){
1337 fCentralityPass->Fill(fCentrality);
1338 fMultCorAfterCuts->Fill(multGlob, multTPC);
1339 fMultvsCentr->Fill(fCentrality, multTPC);
1340 }
1341}
1342//_____________________________________________________________________________
1343void AliAnalysisTaskFlowITSTPCTOFQCSP::SetCentralityParameters(Double_t CentralityMin, Double_t CentralityMax, const char* CentralityMethod)
1344{
1345 // Set a centrality range ]min, max] and define the method to use for centrality selection
1346 fCentralityMin = CentralityMin;
1347 fCentralityMax = CentralityMax;
1348 fkCentralityMethod = CentralityMethod;
1349}
1350//_____________________________________________________________________________
1351void AliAnalysisTaskFlowITSTPCTOFQCSP::SetIDCuts(Double_t minTOFnSigma, Double_t maxTOFnSigma, Double_t minITSnsigmalowpt, Double_t maxITSnsigmalowpt, Double_t minITSnsigmahighpt, Double_t maxITSnsigmahighpt, Double_t minTPCnsigmalowpt, Double_t maxTPCnsigmalowpt, Double_t minTPCnsigmahighpt, Double_t maxTPCnsigmahighpt)
1352{
1353 //Set PID cuts
1354 fminTOFnSigma = minTOFnSigma;
1355 fmaxTOFnSigma = maxTOFnSigma;
1356 fminITSnsigmaLowpT = minITSnsigmalowpt;
1357 fmaxITSnsigmaLowpT = maxITSnsigmalowpt;
1358 fminITSnsigmaHighpT = minITSnsigmahighpt;
1359 fmaxITSnsigmaHighpT = maxITSnsigmahighpt;
1360 fminTPCnsigmaLowpT = minTPCnsigmalowpt;
1361 fmaxTPCnsigmaLowpT = maxTPCnsigmalowpt;
1362 fminTPCnsigmaHighpT = minTPCnsigmahighpt;
1363 fmaxTPCnsigmaHighpT = maxTPCnsigmahighpt;
1364
1365}
1366//_____________________________________________________________________________
1367//_____________________________________________________________________________
1368void AliAnalysisTaskFlowITSTPCTOFQCSP::SetpTCuttrack(Double_t ptmin, Double_t ptmax)
1369{
1370 //Set pt cuts
1371 fpTCutmin = ptmin;
1372 fpTCutmax = ptmax;
164d5623 1373}
4c12b1ec
R		 1374//_____________________________________________________________________________
1375//_____________________________________________________________________________
1376void AliAnalysisTaskFlowITSTPCTOFQCSP::SetHistoForCentralityFlattening(TH1F *h,Double_t minCentr,Double_t maxCentr,Double_t centrRef,Int_t switchTRand){
1377 // set the histo for centrality flattening
1378 // the centrality is flatten in the range minCentr,maxCentr
1379 // if centrRef is zero, the minimum in h within (minCentr,maxCentr) defines the reference
1380 // positive, the value of h(centrRef) defines the reference (-> the centrality distribution might be not flat in the whole desired range)
1381 // negative, h(bin with max in range)*centrRef is used to define the reference (-> defines the maximum loss of events, also in this case the distribution might be not flat)
1382 // switchTRand is used to set the unerflow bin of the histo: if it is < -1 in the analysis the random event selection will be done on using TRandom
1383
1384 if(maxCentr<minCentr){
1385 AliWarning("AliAnalysisCheckCorrdist::Wrong centralities values while setting the histogram for centrality flattening");
1386 }
1387
1388 if(fHistCentrDistr)delete fHistCentrDistr;
1389 fHistCentrDistr=(TH1F*)h->Clone("hCentralityFlat");
1390 fHistCentrDistr->SetTitle("Reference histo for centrality flattening");
1391 Int_t minbin=fHistCentrDistr->FindBin(minCentr*1.00001); // fast if fix bin width
1392 Int_t maxbin=fHistCentrDistr->FindBin(maxCentr*0.9999);
1393 fHistCentrDistr->GetXaxis()->SetRange(minbin,maxbin);
1394 Double_t ref=0.,bincont=0.,binrefwidth=1.;
1395 Int_t binref=0;
1396 if(TMath::Abs(centrRef)<0.0001){
1397 binref=fHistCentrDistr->GetMinimumBin();
1398 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1399 ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
1400 }
1401 else if(centrRef>0.){
1402 binref=h->FindBin(centrRef);
1403 if(binref<1||binref>h->GetNbinsX()){
1404 AliWarning("AliRDHFCuts::Wrong centrality reference value while setting the histogram for centrality flattening");
1405 }
1406 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1407 ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
1408 }
1409 else{
1410 if(centrRef<-1) AliWarning("AliRDHFCuts: with this centrality reference no flattening will be applied");
1411 binref=fHistCentrDistr->GetMaximumBin();
1412 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1413 ref=fHistCentrDistr->GetMaximum()*TMath::Abs(centrRef)/binrefwidth;
1414 }
1415
1416 for(Int_t j=1;j<=h->GetNbinsX();j++){// Now set the "probabilities"
1417 if(h->GetBinLowEdge(j)*1.0001>=minCentr&&h->GetBinLowEdge(j+1)*0.9999<=maxCentr){
1418 bincont=h->GetBinContent(j);
1419 fHistCentrDistr->SetBinContent(j,ref/bincont*h->GetBinWidth(j));
1420 fHistCentrDistr->SetBinError(j,h->GetBinError(j)*ref/bincont);
1421 }
1422 else{
1423 h->SetBinContent(j,1.1);// prob > 1 to assure that events will not be rejected
1424 }
1425 }
1426
1427 fHistCentrDistr->SetBinContent(0,switchTRand);
1428 return;
1429
1430}
1431
1432//-------------------------------------------------
1433Bool_t AliAnalysisTaskFlowITSTPCTOFQCSP::IsEventSelectedForCentrFlattening(Float_t centvalue){
1434 //
1435 // Random event selection, based on fHistCentrDistr, to flatten the centrality distribution
1436 // Can be faster if it was required that fHistCentrDistr covers
1437 // exactly the desired centrality range (e.g. part of the lines below should be done during the
1438 // setting of the histo) and TH1::SetMinimum called
1439 //
1440
1441 if(!fHistCentrDistr) return kTRUE;
1442 // Int_t maxbin=fHistCentrDistr->FindBin(fMaxCentrality*0.9999);
1443 // if(maxbin>fHistCentrDistr->GetNbinsX()){
1444 // AliWarning("AliRDHFCuts: The maximum centrality exceeds the x-axis limit of the histogram for centrality flattening");
1445 // }
1446
1447 Int_t bin=fHistCentrDistr->FindBin(centvalue); // Fast if the histo has a fix bin
1448 Double_t bincont=fHistCentrDistr->GetBinContent(bin);
1449 Double_t centDigits=centvalue-(Int_t)(centvalue*100.)/100.;// this is to extract a random number between 0 and 0.01
1450
1451 if(fHistCentrDistr->GetBinContent(0)<-0.9999){
1452 if(gRandom->Uniform(1.)<bincont)return kTRUE;
1453 return kFALSE;
1454 }
1455
1456 if(centDigits*100.<bincont)return kTRUE;
1457 return kFALSE;
1458
1459}
1460//---------------------------------------------------------------------------
1461
1462
1463//_____________________________________________________________________________
1464