updated for nonHFE reconstruction study
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliAnalysisTaskFlowTPCEMCalQCSP.cxx
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 EMCal          //
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 "AliAnalysisTaskFlowTPCEMCalQCSP.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 "AliEMCALRecoUtils.h"
61 #include "AliEMCALGeometry.h"
62 #include "AliGeomManager.h"
63 #include "stdio.h"
64 #include "TGeoManager.h"
65 #include "iostream"
66 #include "fstream"
67 #include "AliEMCALTrack.h"
68 //#include "AliEMCALTracker.h"
69 #include "AliMagF.h"
70 #include "AliKFParticle.h"
71 #include "AliKFVertex.h"
72 #include "AliPID.h"
73 #include "AliPIDResponse.h"
74 #include "AliHFEcontainer.h"
75 #include "AliHFEcuts.h"
76 #include "AliHFEpid.h"
77 #include "AliHFEpidBase.h"
78 #include "AliHFEpidQAmanager.h"
79 #include "AliHFEtools.h"
80 #include "AliCFContainer.h"
81 #include "AliCFManager.h"
82 #include "AliKFParticle.h"
83 #include "AliKFVertex.h"
84 #include "AliCentrality.h"
85 #include "AliVEvent.h"
86 #include "AliStack.h"
87 #include "AliMCEvent.h"
88 #include "TProfile.h"
89 #include "AliFlowCandidateTrack.h"
90 #include "AliFlowTrackCuts.h"
91 #include "AliFlowEventSimple.h"
92 #include "AliFlowCommonConstants.h"
93 #include "AliFlowEvent.h"
94 #include "TVector3.h"
95 #include "TRandom2.h"
96 #include "AliESDVZERO.h"
97 #include "AliAODVZERO.h"
98 #include "AliPID.h"
99 #include "AliPIDResponse.h"
100 #include "AliFlowTrack.h"
101 #include "AliAnalysisTaskVnV0.h"
102
103
104 class AliFlowTrackCuts;
105
106 using namespace std;
107
108 ClassImp(AliAnalysisTaskFlowTPCEMCalQCSP)
109 //________________________________________________________________________
110   AliAnalysisTaskFlowTPCEMCalQCSP::AliAnalysisTaskFlowTPCEMCalQCSP(const char *name) 
111   : AliAnalysisTaskSE(name)
112 ,fDebug(0)
113 ,fAOD(0)
114 ,fGeom(0)
115 ,fOutputList(0)
116 ,fCuts(0)
117 ,fIdentifiedAsOutInz(kFALSE)
118 ,fPassTheEventCut(kFALSE)
119 ,fCFM(0)
120 ,fPID(0)
121 ,fPIDqa(0)
122 ,fCutsRP(0)     // track cuts for reference particles
123 ,fNullCuts(0) // dummy cuts for flow event tracks
124 ,fFlowEvent(0) //! flow events (one for each inv mass band)
125 ,fkCentralityMethod(0)
126 ,fCentrality(0)
127 ,fCentralityMin(0)
128 ,fCentralityMax(0)
129 ,fInvmassCut(0)
130 ,fpTCut(0)
131 ,fTrigger(0)
132 ,fPhi(0)
133 ,fEta(0)
134 ,fVZEROA(0)
135 ,fVZEROC(0)
136 ,fTPCM(0)
137 ,fNoEvents(0)
138 ,fTrkEovPBef(0)
139 //,fdEdxBef(0)
140 ,fInclusiveElecPt(0)
141 ,fTPCnsigma(0)
142 ,fTPCnsigmaAft(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 ,fminTPC(-1)
152 ,fmaxTPC(3)
153 ,fminEovP(0.8)
154 ,fmaxEovP(1.2)
155 ,fminM20(0.03)
156 ,fmaxM20(0.3)
157 ,fminM02(0.03)
158 ,fmaxM02(0.5)
159 ,fDispersion(1)
160 ,fMultCorAfterCuts(0)
161 ,fMultvsCentr(0)
162 ,fSubEventDPhiv2(0)
163 ,EPVzA(0)
164 ,EPVzC(0)
165 ,EPTPC(0)
166 ,fV2Phi(0)
167 ,fSparseElectronHadron(0)
168 ,fvertex(0)
169 ,fMultCorBeforeCuts(0)
170 ,fSideBandsFlow(kFALSE)
171 ,fPhiminusPsi(kFALSE)
172 ,fFlowEventCont(0) //! flow events (one for each inv mass band)
173 ,fpurity(kFALSE)
174 ,fSparseElectronpurity(0)
175 {
176   //Named constructor
177
178   fPID = new AliHFEpid("hfePid");
179   // Define input and output slots here
180   // Input slot #0 works with a TChain
181   DefineInput(0, TChain::Class());
182   // Output slot #0 id reserved by the base class for AOD
183   // Output slot #1 writes into a TH1 container
184   // DefineOutput(1, TH1I::Class());
185   DefineOutput(1, TList::Class());
186   DefineOutput(2, AliFlowEventSimple::Class());
187 if(fSideBandsFlow){
188     DefineOutput(3, AliFlowEventSimple::Class());
189 }
190  //  DefineOutput(3, TTree::Class());
191 }
192
193 //________________________________________________________________________
194 AliAnalysisTaskFlowTPCEMCalQCSP::AliAnalysisTaskFlowTPCEMCalQCSP() 
195   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskFlowTPCEMCalQCSP")
196 ,fDebug(0)
197 ,fAOD(0)
198 ,fGeom(0)
199 ,fOutputList(0)
200 ,fCuts(0)
201 ,fIdentifiedAsOutInz(kFALSE)
202 ,fPassTheEventCut(kFALSE)
203 ,fCFM(0)
204 ,fPID(0)
205 ,fPIDqa(0)
206 ,fCutsRP(0)     // track cuts for reference particles
207 ,fNullCuts(0) // dummy cuts for flow event tracks
208 ,fFlowEvent(0) //! flow events (one for each inv mass band)
209 ,fkCentralityMethod(0)
210 ,fCentrality(0)
211 ,fCentralityMin(0)
212 ,fCentralityMax(0)
213 ,fInvmassCut(0)
214 ,fpTCut(0)
215 ,fTrigger(0)
216 ,fPhi(0)
217 ,fEta(0)
218 ,fVZEROA(0)
219 ,fVZEROC(0)
220 ,fTPCM(0)
221 ,fNoEvents(0)
222 ,fTrkEovPBef(0)
223 //,fdEdxBef(0)
224 ,fInclusiveElecPt(0)
225 ,fTPCnsigma(0)
226 ,fTPCnsigmaAft(0)
227 ,fCentralityPass(0)
228 ,fCentralityNoPass(0)
229 ,fInvmassLS1(0)
230 ,fInvmassULS1(0)
231 ,fPhotoElecPt(0)
232 ,fSemiInclElecPt(0)
233 ,fULSElecPt(0)
234 ,fLSElecPt(0)
235 ,fminTPC(-1)
236 ,fmaxTPC(3)
237 ,fminEovP(0.8)
238 ,fmaxEovP(1.2)
239 ,fminM20(0.03)
240 ,fmaxM20(0.3)
241 ,fminM02(0.03)
242 ,fmaxM02(0.5)
243 ,fDispersion(1)
244 ,fMultCorAfterCuts(0)
245 ,fMultvsCentr(0)
246 ,fSubEventDPhiv2(0)
247 ,EPVzA(0)
248 ,EPVzC(0)
249 ,EPTPC(0)
250 ,fV2Phi(0)
251 ,fSparseElectronHadron(0)
252 ,fvertex(0)
253 ,fMultCorBeforeCuts(0)
254 ,fSideBandsFlow(kFALSE)
255 ,fPhiminusPsi(kFALSE)
256 ,fFlowEventCont(0) //! flow events (one for each inv mass band)
257 ,fpurity(kFALSE)
258 ,fSparseElectronpurity(0)
259 {
260   //Default constructor
261   fPID = new AliHFEpid("hfePid");
262   // Constructor
263   // Define input and output slots here
264   // Input slot #0 works with a TChain
265   DefineInput(0, TChain::Class());
266   // Output slot #0 id reserved by the base class for AOD
267   // Output slot #1 writes into a TH1 container
268   // DefineOutput(1, TH1I::Class());
269   DefineOutput(1, TList::Class());
270   DefineOutput(2, AliFlowEventSimple::Class());
271     //  DefineOutput(3, TTree::Class());
272 if(fSideBandsFlow){
273     DefineOutput(3, AliFlowEventSimple::Class());
274 }
275     //DefineOutput(3, TTree::Class());
276 }
277 //_________________________________________
278
279 AliAnalysisTaskFlowTPCEMCalQCSP::~AliAnalysisTaskFlowTPCEMCalQCSP()
280 {
281   //Destructor 
282
283   delete fOutputList;
284   delete fGeom;
285   delete fPID;
286   delete fCFM;
287   delete fPIDqa; 
288   if (fOutputList) delete fOutputList;
289   if (fFlowEvent) delete fFlowEvent;
290   if (fFlowEventCont) delete fFlowEventCont;
291
292 }
293 //_________________________________________
294
295 void AliAnalysisTaskFlowTPCEMCalQCSP::UserExec(Option_t*)
296 {
297   //Main loop
298   //Called for each event
299
300   // create pointer to event
301
302   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
303   if (!fAOD)
304   {
305     printf("ERROR: fAOD not available\n");
306     return;
307   }
308
309   if(!fCuts)
310   {
311     AliError("HFE cuts not available");
312     return;
313   }
314
315   if(!fPID->IsInitialized())
316   { 
317     // Initialize PID with the given run number
318     AliWarning("PID not initialised, get from Run no");
319     fPID->InitializePID(fAOD->GetRunNumber());
320   }
321     
322   // cout << "kTrigger   ==   " << fTrigger <<endl;
323     
324     if(fTrigger==0){
325 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral) ) return;
326     }
327     if(fTrigger==1){
328 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kSemiCentral))) return;
329     }
330     if(fTrigger==2){
331 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMCEGA) ) return;
332     }
333     if(fTrigger==3){
334 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB) ) return;
335     }
336   
337     
338 //---------------CENTRALITY AND EVENT SELECTION-----------------------
339    Int_t fNOtrks =  fAOD->GetNumberOfTracks();
340     Float_t vtxz = -999;
341     const AliAODVertex* trkVtx = fAOD->GetPrimaryVertex();
342     if (!trkVtx || trkVtx->GetNContributors()<=0)return;
343     TString vtxTtl = trkVtx->GetTitle();
344     if (!vtxTtl.Contains("VertexerTracks"))return;
345     const AliAODVertex* spdVtx = fAOD->GetPrimaryVertexSPD();
346     if (!spdVtx || spdVtx->GetNContributors()<=0)return;
347     if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5)return;
348     vtxz = trkVtx->GetZ();
349     if(TMath::Abs(vtxz)>10)return;
350     fvertex->Fill(vtxz);
351
352 // Event cut
353     if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fAOD)) return;
354     if(fNOtrks<2) return;
355     
356     Bool_t pass = kFALSE; //to select centrality
357     CheckCentrality(fAOD,pass);
358     if(!pass)return;
359     
360     
361   fNoEvents->Fill(0);
362   PlotVZeroMultiplcities(fAOD);
363
364   SetNullCuts(fAOD);
365   PrepareFlowEvent(fAOD->GetNumberOfTracks(),fFlowEvent);    //Calculate event plane Qvector and EP resolution for inclusive
366
367     if(fSideBandsFlow){
368   PrepareFlowEvent(fAOD->GetNumberOfTracks(),fFlowEventCont);    //Calculate event plane Qvector and EP resolution for inclusive
369     }
370     
371   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
372   if(!pidResponse)
373   {
374     AliDebug(1, "Using default PID Response");
375     pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class()); 
376   }
377
378   fPID->SetPIDResponse(pidResponse);
379   fCFM->SetRecEventInfo(fAOD);
380
381   // Look for kink mother
382   Int_t numberofvertices = fAOD->GetNumberOfVertices();
383   Double_t listofmotherkink[numberofvertices];
384   Int_t numberofmotherkink = 0;
385   for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
386     AliAODVertex *aodvertex = fAOD->GetVertex(ivertex);
387     if(!aodvertex) continue;
388     if(aodvertex->GetType()==AliAODVertex::kKink) {
389       AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
390       if(!mother) continue;
391       Int_t idmother = mother->GetID();
392       listofmotherkink[numberofmotherkink] = idmother;
393       //printf("ID %d\n",idmother);
394       numberofmotherkink++;
395     }
396   }
397     
398 //=============================================V0EP from Alex======================================================================
399     Double_t qxEPa = 0, qyEPa = 0;
400     Double_t qxEPc = 0, qyEPc = 0;
401     
402     Double_t evPlAngV0A = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 8, 2, qxEPa, qyEPa);
403     Double_t evPlAngV0C = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 9, 2, qxEPc, qyEPc);
404     
405     
406     Double_t Qx2 = 0, Qy2 = 0;
407     
408     for (Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++){
409         
410         AliAODTrack* aodTrack = fAOD->GetTrack(iT);
411         
412         if (!aodTrack)
413             continue;
414         
415         if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < 70) || (aodTrack->Pt() >= 20.0))
416             continue;
417         
418         if (!aodTrack->TestFilterBit(128))
419             continue;
420         
421         Qx2 += TMath::Cos(2*aodTrack->Phi());
422         Qy2 += TMath::Sin(2*aodTrack->Phi());
423     }
424     
425     Double_t evPlAngTPC = TMath::ATan2(Qy2, Qx2)/2.;
426     
427     EPVzA->Fill(evPlAngV0A);
428     EPVzC->Fill(evPlAngV0C);
429     EPTPC->Fill(evPlAngTPC);
430     
431     fSubEventDPhiv2->Fill(0.5, TMath::Cos(2.*(evPlAngV0A-evPlAngTPC))); // vzeroa - tpc
432     fSubEventDPhiv2->Fill(1.5, TMath::Cos(2.*(evPlAngV0A-evPlAngV0C))); // vzeroa - vzeroc
433     fSubEventDPhiv2->Fill(2.5, TMath::Cos(2.*(evPlAngV0C-evPlAngTPC))); // tpc - vzeroc
434 //====================================================================================================================
435     
436     
437  AliAODTrack *track = NULL;
438
439 // Track loop 
440  for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) 
441  {
442    track = fAOD->GetTrack(iTracks);
443    if (!track) 
444    {
445      printf("ERROR: Could not receive track %d\n", iTracks);
446      continue;
447    }
448      
449    if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue;  // TESTBIT FOR AOD double Counting
450 //----------hfe begin---------
451    if(track->Eta()<-0.7 || track->Eta()>0.7)    continue;    //eta cuts on candidates
452
453    // RecKine: ITSTPC cuts  
454    if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
455
456    // Reject kink mother
457    Bool_t kinkmotherpass = kTRUE;
458    for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
459      if(track->GetID() == listofmotherkink[kinkmother]) {
460        kinkmotherpass = kFALSE;
461        continue;
462      }
463    }
464    if(!kinkmotherpass) continue;
465
466    // RecPrim
467  //  if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;  //deleted for DCA absence
468    // HFEcuts: ITS layers cuts
469    if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
470    // HFE cuts: TPC PID cleanup
471    if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
472      
473    Double_t fClsE = -999, p = -999, fEovP=-999, pt = -999, fTPCnSigma=0;
474    // Track extrapolation
475    Int_t fClsId = track->GetEMCALcluster();
476    if(fClsId < 0) continue;
477    AliAODCaloCluster *cluster = fAOD->GetCaloCluster(fClsId);
478    if(TMath::Abs(cluster->GetTrackDx()) > 0.05 || TMath::Abs(cluster->GetTrackDz()) > 0.05) continue;
479
480    pt = track->Pt();         //pt track after cuts
481    if(pt<fpTCut) continue;
482    fClsE = cluster->E();
483    p = track->P();
484   // dEdx = track->GetTPCsignal();
485    fEovP = fClsE/p;
486    fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
487    Double_t m20 =cluster->GetM20();
488    Double_t m02 =cluster->GetM02();
489    Double_t disp=cluster->GetDispersion();
490    if(fTPCnSigma >= -1 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,fEovP);
491    fTPCnsigma->Fill(p,fTPCnSigma);
492 //   fdEdxBef->Fill(p,dEdx);
493    Double_t eta = track->Eta(); 
494    Double_t phi = track->Phi();
495 //-----------------------Phiminupsi method to remove the contamination-----------------------------------------------
496 //-----------------------fTPCnSigma < -3.5 hadrons will be selected from this region--------------------------
497      Float_t dPhi_aeh = TVector2::Phi_0_2pi(phi - evPlAngV0A);
498      if(dPhi_aeh > TMath::Pi()) dPhi_aeh = dPhi_aeh - TMath::Pi();
499      Float_t dPhi_ceh = TVector2::Phi_0_2pi(phi - evPlAngV0C);
500      if(dPhi_ceh > TMath::Pi()) dPhi_ceh = dPhi_ceh - TMath::Pi();
501
502      if(fPhiminusPsi){
503          Double_t valueElh[8] = {
504              pt,
505              fEovP,
506              fTPCnSigma,
507              m20,
508              m02,
509              disp,
510              dPhi_aeh,
511              dPhi_ceh};
512          fSparseElectronHadron->Fill(valueElh);
513      }
514 //----------------------------------------------------------------------------------------------------------
515 //---------------------------From here usual electron selection---------------------------------------------
516 //----------------------------------------------------------------------------------------------------------
517 if(m20 < fminM20 || m20 > fmaxM20) continue;
518 if(m02 < fminM02 || m02 > fmaxM02) continue;
519 if(disp > fDispersion ) continue;
520 //---------------------------------for purity---------------------------------------------------------------
521      if(fpurity){
522          Double_t valuepurity[3] = {
523              pt,
524              fEovP,
525              fTPCnSigma};
526          fSparseElectronpurity->Fill(valuepurity);
527      }
528 //----------------------------------------------------------------------------------------------------------
529 //----------------------------------------------------------------------------------------------------------
530 if(fTPCnSigma < fminTPC || fTPCnSigma > fmaxTPC) continue; //cuts on nsigma tpc and EoP
531 //===============================Flow Event for Contamination=============================================
532      if(fSideBandsFlow){
533      if(fEovP>0 && fEovP<0.6){
534          AliFlowTrack *sTrackCont = new AliFlowTrack();
535          sTrackCont->Set(track);
536          sTrackCont->SetID(track->GetID());
537          sTrackCont->SetForRPSelection(kFALSE);
538          sTrackCont->SetForPOISelection(kTRUE);
539          sTrackCont->SetMass(2637);
540          for(int iRPs=0; iRPs!=fFlowEventCont->NumberOfTracks(); ++iRPs)
541          {
542              //   cout << " no of rps " << iRPs << endl;
543              AliFlowTrack *iRPCont = dynamic_cast<AliFlowTrack*>(fFlowEventCont->GetTrack( iRPs ));
544              if (!iRPCont) continue;
545              if (!iRPCont->InRPSelection()) continue;
546              if( sTrackCont->GetID() == iRPCont->GetID())
547              {
548                  if(fDebug) printf(" was in RP set");
549                  //       cout << sTrack->GetID() <<"   ==  " << iRP->GetID() << " was in RP set" <<endl;
550                  iRPCont->SetForRPSelection(kFALSE);
551                  fFlowEventCont->SetNumberOfRPs(fFlowEventCont->GetNumberOfRPs() - 1);
552              }
553          } //end of for loop on RPs
554          fFlowEventCont->InsertTrack(((AliFlowTrack*) sTrackCont));
555      }
556      }
557 //==========================================================================================================
558 //===============================From here eovP cut is used fro QC, SP and EPV0=============================
559 if(fEovP < fminEovP || fEovP >fmaxEovP) continue;
560 //==========================================================================================================
561 //============================Event Plane Method with V0====================================================
562      Double_t v2PhiV0A = TMath::Cos(2*(phi - evPlAngV0A));
563      Double_t v2PhiV0C = TMath::Cos(2*(phi - evPlAngV0C));
564      Double_t v2Phi[3] = {
565          v2PhiV0A,
566          v2PhiV0C,
567          pt};
568      fV2Phi->Fill(v2Phi);
569 //=========================================================================================================
570    fTPCnsigmaAft->Fill(p,fTPCnSigma);
571    fInclusiveElecPt->Fill(pt); 
572    fPhi->Fill(phi); 
573    fEta->Fill(eta); 
574 //----------------------Flow of Inclusive Electrons--------------------------------------------------------
575    AliFlowTrack *sTrack = new AliFlowTrack();
576      sTrack->Set(track);
577      sTrack->SetID(track->GetID());
578      sTrack->SetForRPSelection(kFALSE);
579      sTrack->SetForPOISelection(kTRUE);
580      sTrack->SetMass(263732);
581      for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs)
582      {
583       //   cout << " no of rps " << iRPs << endl;
584          AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
585          if (!iRP) continue;
586          if (!iRP->InRPSelection()) continue;
587         if( sTrack->GetID() == iRP->GetID())
588         {
589           if(fDebug) printf(" was in RP set");
590      //       cout << sTrack->GetID() <<"   ==  " << iRP->GetID() << " was in RP set" <<endl;
591             iRP->SetForRPSelection(kFALSE);
592             fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
593         }
594       } //end of for loop on RPs
595    fFlowEvent->InsertTrack(((AliFlowTrack*) sTrack));
596     
597 //----------------------Selection and Flow of Photonic Electrons-----------------------------
598    Bool_t fFlagPhotonicElec = kFALSE;
599    SelectPhotonicElectron(iTracks,track,fFlagPhotonicElec);
600    if(fFlagPhotonicElec){fPhotoElecPt->Fill(pt);}
601    // Semi inclusive electron 
602    if(!fFlagPhotonicElec){fSemiInclElecPt->Fill(pt);}
603      
604  }//end loop on track
605
606  PostData(1, fOutputList);
607  PostData(2, fFlowEvent);
608     if(fSideBandsFlow){
609  PostData(3, fFlowEventCont);
610     }
611
612  //----------hfe end---------
613 }
614 //_________________________________________
615 void AliAnalysisTaskFlowTPCEMCalQCSP::SelectPhotonicElectron(Int_t itrack,const AliAODTrack *track, Bool_t &fFlagPhotonicElec)
616 {
617   //Identify non-heavy flavour electrons using Invariant mass method
618
619   Bool_t flagPhotonicElec = kFALSE;
620
621   for(Int_t jTracks = 0; jTracks<fAOD->GetNumberOfTracks(); jTracks++){
622     AliAODTrack *trackAsso = fAOD->GetTrack(jTracks);
623     if (!trackAsso) {
624       printf("ERROR: Could not receive track %d\n", jTracks);
625       continue;
626     }
627     //  if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue;  // TESTBIT FOR AOD double Counting
628       if(!trackAsso->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
629       if((!(trackAsso->GetStatus()&AliESDtrack::kITSrefit)|| (!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
630
631       
632     if(jTracks == itrack) continue;
633     Double_t ptAsso=-999., nsigma=-999.0;
634     Double_t mass=-999., width = -999;
635     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
636
637   
638     ptAsso = trackAsso->Pt();
639     Short_t chargeAsso = trackAsso->Charge();
640     Short_t charge = track->Charge();
641     nsigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
642     
643     if(trackAsso->GetTPCNcls() < 80) continue;
644     if(nsigma < -3 || nsigma > 3) continue;
645     if(trackAsso->Eta()<-0.9 || trackAsso->Eta()>0.9) continue;
646     if(ptAsso <0.3) continue;
647
648     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
649     if(charge>0) fPDGe1 = -11;
650     if(chargeAsso>0) fPDGe2 = -11;
651
652     if(charge == chargeAsso) fFlagLS = kTRUE;
653     if(charge != chargeAsso) fFlagULS = kTRUE;
654
655     AliKFParticle ge1 = AliKFParticle(*track, fPDGe1);
656     AliKFParticle ge2 = AliKFParticle(*trackAsso, fPDGe2);
657     AliKFParticle recg(ge1, ge2);
658
659     if(recg.GetNDF()<1) continue;
660     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
661     if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
662     recg.GetMass(mass,width);
663
664     if(fFlagLS) fInvmassLS1->Fill(mass);
665     if(fFlagULS) fInvmassULS1->Fill(mass);
666            
667        if(mass<fInvmassCut){
668        if(fFlagULS){fULSElecPt->Fill(track->Pt());}
669        if(fFlagLS){fLSElecPt->Fill(track->Pt());}
670        }
671     
672     if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
673       flagPhotonicElec = kTRUE;
674     }
675   }//track loop
676   fFlagPhotonicElec = flagPhotonicElec;
677 }
678 //___________________________________________
679 void AliAnalysisTaskFlowTPCEMCalQCSP::UserCreateOutputObjects()
680 {
681   //Create histograms
682
683   //----------hfe initialising begin---------
684    fNullCuts = new AliFlowTrackCuts("null_cuts");
685  
686    AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
687   cc->SetNbinsMult(10000);
688   cc->SetMultMin(0);
689   cc->SetMultMax(10000);
690
691   cc->SetNbinsPt(100);
692   cc->SetPtMin(0);
693   cc->SetPtMax(50);
694
695   cc->SetNbinsPhi(180);
696   cc->SetPhiMin(0.0);
697   cc->SetPhiMax(TMath::TwoPi());
698
699   cc->SetNbinsEta(200);
700   cc->SetEtaMin(-7.0);
701   cc->SetEtaMax(+7.0);
702
703   cc->SetNbinsQ(500);
704   cc->SetQMin(0.0);
705   cc->SetQMax(3.0);
706
707   //--------Initialize PID
708   fPID->SetHasMCData(kFALSE);
709   if(!fPID->GetNumberOfPIDdetectors()) 
710   {
711     fPID->AddDetector("TPC", 0);
712     fPID->AddDetector("EMCAL", 1);
713   }
714
715   fPID->SortDetectors(); 
716   fPIDqa = new AliHFEpidQAmanager();
717   fPIDqa->Initialize(fPID);
718
719   //--------Initialize correction Framework and Cuts
720   fCFM = new AliCFManager;
721   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
722   fCFM->SetNStepParticle(kNcutSteps);
723   for(Int_t istep = 0; istep < kNcutSteps; istep++)
724     fCFM->SetParticleCutsList(istep, NULL);
725
726   if(!fCuts){
727     AliWarning("Cuts not available. Default cuts will be used");
728     fCuts = new AliHFEcuts;
729     fCuts->CreateStandardCuts();
730   }
731
732   fCuts->SetAOD(); 
733   fCuts->Initialize(fCFM);
734   //----------hfe initialising end--------
735   //---------Output Tlist
736   fOutputList = new TList();
737   fOutputList->SetOwner();
738   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
739
740   fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
741   fOutputList->Add(fNoEvents);
742
743   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma before HFE pid",1000,0,50,200,-10,10);
744   fOutputList->Add(fTPCnsigma);
745
746   fTPCnsigmaAft = new TH2F("fTPCnsigmaAft", "TPC - n sigma after HFE pid",1000,0,50,200,-10,10);
747   fOutputList->Add(fTPCnsigmaAft);
748
749   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",1000,0,50,100,0,2);
750   fOutputList->Add(fTrkEovPBef);
751
752 //    fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",1000,0,50,150,0,150);
753  //     fOutputList->Add(fdEdxBef);
754    
755   fInclusiveElecPt = new TH1F("fInclElecPt", "Inclusive electron pt",1000,0,100);
756   fOutputList->Add(fInclusiveElecPt);
757
758   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",1000,0,100);
759   fOutputList->Add(fPhotoElecPt);
760   //    
761   fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",1000,0,100);
762   fOutputList->Add(fSemiInclElecPt);
763
764      fULSElecPt = new TH1F("fULSElecPt", "ULS electron pt",1000,0,100);
765      fOutputList->Add(fULSElecPt);
766
767     fLSElecPt = new TH1F("fLSElecPt", "LS electron pt",1000,0,100);
768     fOutputList->Add(fLSElecPt);
769
770   fInvmassLS1 = new TH1F("fInvmassLS1", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
771   fOutputList->Add(fInvmassLS1);
772
773   fInvmassULS1 = new TH1F("fInvmassULS1", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
774   fOutputList->Add(fInvmassULS1);
775
776   fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
777   fOutputList->Add(fCentralityPass);
778
779   fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
780   fOutputList->Add(fCentralityNoPass);
781
782   fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
783   fOutputList->Add(fPhi);
784
785   fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
786   fOutputList->Add(fEta);
787
788   fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
789   fOutputList->Add(fVZEROA);
790
791   fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
792   fOutputList->Add(fVZEROC);
793
794   fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
795   fOutputList->Add(fTPCM);
796
797   fvertex = new TH1D("fvertex", "vertex distribution", 300, -15,15);
798   fOutputList->Add(fvertex);
799   
800   fMultCorBeforeCuts = new TH2F("fMultCorBeforeCuts", "TPC vs Global multiplicity (Before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
801   fOutputList->Add(fMultCorBeforeCuts);
802
803   fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
804   fOutputList->Add(fMultCorAfterCuts);
805
806   fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 100, 0., 100, 100, 0, 3000);
807   fOutputList->Add(fMultvsCentr);
808     
809 //----------------------------------------------------------------------------
810     EPVzA = new TH1D("EPVzA", "EPVzA", 80, -2, 2);
811     fOutputList->Add(EPVzA);
812     EPVzC = new TH1D("EPVzC", "EPVzC", 80, -2, 2);
813     fOutputList->Add(EPVzC);
814     EPTPC = new TH1D("EPTPC", "EPTPC", 80, -2, 2);
815     fOutputList->Add(EPTPC);
816 //----------------------------------------------------------------------------
817     fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 3, 0, 3);
818     fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
819     fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
820     fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
821     fOutputList->Add(fSubEventDPhiv2);
822 //================================Event Plane with VZERO=====================
823     const Int_t nPtBins = 10;
824     Double_t binsPt[nPtBins+1] = {0, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
825     // v2A, v2C, pt
826     Int_t    bins[3] = {  50,  50, nPtBins};
827     Double_t xmin[3] = { -1., -1.,   0};
828     Double_t xmax[3] = {  1.,  1.,   8};
829     fV2Phi = new THnSparseF("fV2Phi", "v2A:v2C:pt", 3, bins, xmin, xmax);
830     // Set bin limits for axes which are not standard binned
831     fV2Phi->SetBinEdges(2, binsPt);
832     // set axes titles
833     fV2Phi->GetAxis(0)->SetTitle("v_{2} (V0A)");
834     fV2Phi->GetAxis(1)->SetTitle("v_{2} (V0C)");
835     fV2Phi->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
836     fOutputList->Add(fV2Phi);
837 //----------------------------------------------------------------------------
838     if(fPhiminusPsi){
839     Int_t binsvElectH[8]={ 600,  200, 200 ,100,  100,  100,   10,          10}; //pt, E/p,TPCnSigma,M20,M02,Disp Phi-psiV0A ,Phi-PsiV0C,eta (commented)
840     Double_t xminvElectH[8]={0,    0, -10 ,  0,    0,    0,    0,           0};
841     Double_t xmaxvElectH[8]={20,   2,  10 ,  2,    2,    2,  TMath::Pi(), TMath::Pi()};
842     fSparseElectronHadron = new THnSparseD("ElectronHadron","ElectronHadron",8,binsvElectH,xminvElectH,xmaxvElectH);
843     fSparseElectronHadron->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
844     fSparseElectronHadron->GetAxis(1)->SetTitle("EovP");
845     fSparseElectronHadron->GetAxis(2)->SetTitle("TPCnSigma");
846     fSparseElectronHadron->GetAxis(3)->SetTitle("M20");
847     fSparseElectronHadron->GetAxis(4)->SetTitle("M02");
848     fSparseElectronHadron->GetAxis(5)->SetTitle("Disp");
849     fSparseElectronHadron->GetAxis(6)->SetTitle("phiminuspsi V0A");
850     fSparseElectronHadron->GetAxis(7)->SetTitle("phiminuspsi V0C");
851     fOutputList->Add(fSparseElectronHadron);
852     }
853 //----------------------------------------------------------------------------
854 //----------------------------------------------------------------------------
855     if(fpurity){
856     Int_t binsvpurity[3]={   600,200, 200}; //pt, E/p,TPCnSigma
857     Double_t xminvpurity[3]={0,    0, -10};
858     Double_t xmaxvpurity[3]={20,   2,  10};
859     fSparseElectronpurity = new THnSparseD("Electronpurity","Electronpurity",3,binsvpurity,xminvpurity,xmaxvpurity);
860     fSparseElectronpurity->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
861     fSparseElectronpurity->GetAxis(1)->SetTitle("EovP");
862     fSparseElectronpurity->GetAxis(2)->SetTitle("TPCnSigma");
863     fOutputList->Add(fSparseElectronpurity);
864     }
865 //----------------------------------------------------------------------------
866
867   PostData(1,fOutputList);
868  // create and post flowevent
869   fFlowEvent = new AliFlowEvent(10000);
870   PostData(2, fFlowEvent);
871     
872     if(fSideBandsFlow){
873     fFlowEventCont = new AliFlowEvent(10000);
874     PostData(3, fFlowEventCont);
875     }
876  }
877
878 //________________________________________________________________________
879 void AliAnalysisTaskFlowTPCEMCalQCSP::Terminate(Option_t *)
880 {
881   // Info("Terminate");
882   AliAnalysisTaskSE::Terminate();
883 }
884 //_____________________________________________________________________________
885 template <typename T> void AliAnalysisTaskFlowTPCEMCalQCSP::PlotVZeroMultiplcities(const T* event) const
886 {
887   // QA multiplicity plots
888   fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
889   fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
890 }
891 //_____________________________________________________________________________
892 template <typename T> void AliAnalysisTaskFlowTPCEMCalQCSP::SetNullCuts(T* event)
893 {
894  //Set null cuts
895     if (fDebug) cout << " fCutsRP " << fCutsRP << endl;
896     fCutsRP->SetEvent(event, MCEvent());
897     fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
898     fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
899     fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
900     fNullCuts->SetEvent(event, MCEvent());
901 }
902 //_____________________________________________________________________________
903 void AliAnalysisTaskFlowTPCEMCalQCSP::PrepareFlowEvent(Int_t iMulti, AliFlowEvent *FlowEv) const 
904 {
905  //Prepare flow events
906    FlowEv->ClearFast();
907    FlowEv->Fill(fCutsRP, fNullCuts);
908    FlowEv->SetReferenceMultiplicity(iMulti);
909    FlowEv->DefineDeadZone(0, 0, 0, 0);
910  //  FlowEv->TagSubeventsInEta(-0.7, 0, 0, 0.7);
911 }
912 //_____________________________________________________________________________
913 Bool_t AliAnalysisTaskFlowTPCEMCalQCSP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
914 {
915   // Check single track cuts for a given cut step
916   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
917   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
918   return kTRUE;
919 }
920 //_________________________________________
921 void AliAnalysisTaskFlowTPCEMCalQCSP::CheckCentrality(AliAODEvent* event, Bool_t &centralitypass)
922 {
923   // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
924   if (!fkCentralityMethod) AliFatal("No centrality method set! FATAL ERROR!");
925   fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethod);
926 //  cout << "--------------Centrality evaluated-------------------------"<<endl;
927   if ((fCentrality <= fCentralityMin) || (fCentrality > fCentralityMax))
928   {
929     fCentralityNoPass->Fill(fCentrality);
930 //    cout << "--------------Fill no pass-----"<< fCentrality <<"--------------------"<<endl;
931     centralitypass = kFALSE;
932   }else
933   { 
934 //    cout << "--------------Fill pass----"<< fCentrality <<"---------------------"<<endl;
935     centralitypass = kTRUE; 
936   }
937 //to remove the bias introduced by multeplicity outliers---------------------
938     Float_t centTrk = event->GetCentrality()->GetCentralityPercentile("TRK");
939     Float_t centv0 = event->GetCentrality()->GetCentralityPercentile("V0M");
940
941     if (TMath::Abs(centv0 - centTrk) > 5.0){
942         centralitypass = kFALSE;
943         fCentralityNoPass->Fill(fCentrality);
944      }
945     const Int_t nGoodTracks = event->GetNumberOfTracks();
946     
947     Float_t multTPC(0.); // tpc mult estimate
948     Float_t multGlob(0.); // global multiplicity
949     for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
950         AliAODTrack* trackAOD = event->GetTrack(iTracks);
951         if (!trackAOD) continue;
952         if (!(trackAOD->TestFilterBit(1))) continue;
953         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;
954         multTPC++;
955     }
956     for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
957         AliAODTrack* trackAOD = event->GetTrack(iTracks);
958         if (!trackAOD) continue;
959         if (!(trackAOD->TestFilterBit(16))) continue;
960         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;
961         Double_t b[2] = {-99., -99.};
962         Double_t bCov[3] = {-99., -99., -99.};
963         if (!(trackAOD->PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
964         if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
965         multGlob++;
966     } //track loop
967     //     printf(" mult TPC %.2f, mult Glob %.2f \n", multTPC, multGlob);
968  //   if(! (multTPC > (-40.3+1.22*multGlob) && multTPC < (32.1+1.59*multGlob))){  2010
969     if(! (multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))){ 
970         centralitypass = kFALSE;
971         fCentralityNoPass->Fill(fCentrality);
972     }//2011
973     fMultCorBeforeCuts->Fill(multGlob, multTPC);
974
975     if(centralitypass){
976     fCentralityPass->Fill(fCentrality);
977     fMultCorAfterCuts->Fill(multGlob, multTPC);
978     fMultvsCentr->Fill(fCentrality, multTPC);
979     }
980 }
981 //_____________________________________________________________________________
982 void AliAnalysisTaskFlowTPCEMCalQCSP::SetCentralityParameters(Double_t CentralityMin, Double_t CentralityMax, const char* CentralityMethod)
983 {
984   // Set a centrality range ]min, max] and define the method to use for centrality selection
985   fCentralityMin = CentralityMin;
986   fCentralityMax = CentralityMax;
987   fkCentralityMethod = CentralityMethod;
988 }
989 //_____________________________________________________________________________
990 void AliAnalysisTaskFlowTPCEMCalQCSP::SetIDCuts(Double_t minTPC, Double_t maxTPC, Double_t minEovP, Double_t maxEovP, Double_t minM20, Double_t maxM20, Double_t minM02, Double_t maxM02, Double_t Dispersion)
991 {
992     //Set ID cuts
993     fminTPC = minTPC;
994     fmaxTPC = maxTPC;
995     fminEovP = minEovP;
996     fmaxEovP = maxEovP;
997     fminM20 = minM20;
998     fmaxM20 = maxM20;
999     fminM02 = minM02;
1000     fmaxM02 = maxM02;
1001     fDispersion = Dispersion;
1002 }
1003 //_____________________________________________________________________________
1004
1005