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