]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG2/FLOW/macros/AddTaskFlow.C
settable rang phi histogram + few other fixes
[u/mrichter/AliRoot.git] / PWG2 / FLOW / macros / AddTaskFlow.C
index 5537c5a4e0860f57f490ae9017208166d2d9d139..444df54a0e97d1c3b51c7dac07850c59d1ecb914 100644 (file)
 //
 /////////////////////////////////////////////////////////////////////////////////////////////
 
-// Define the range for eta subevents (for SP method)
-Double_t minA = -0.9;
-Double_t maxA = -0.01;
-Double_t minB = 0.01;
-Double_t maxB = 0.9;
-
-// use physics selection class
-Bool_t  UsePhysicsSelection = kFALSE;
-
-// SETTING THE CUTS
-
-//----------Event selection----------
-Bool_t UseMultCutforESD = kTRUE;
-//Bool_t UseMultCutforESD = kFALSE;
-const Int_t multminESD = 1;  //used for CORRFW cuts 
-const Int_t multmaxESD = 1000000; //used for CORRFW cuts 
-
-Bool_t requireVtxCuts = kTRUE;
-Double_t vertexXmin = -1.0;
-Double_t vertexXmax = 1.0;
-Double_t vertexYmin = -1.0;
-Double_t vertexYmax = 1.0;
-Double_t vertexZmin = -15.0; //-1.e99;
-Double_t vertexZmax = 15.0; //1.e99;
-
-//Bool_t UseMultCut = kFALSE;
-Bool_t UseMultCut = kTRUE;
-const Int_t multmin = 10;     //used for AliFlowEventSimple (to set the centrality)
-const Int_t multmax = 10000;     //used for AliFlowEventSimple (to set the centrality)
-//const Int_t multmin = 10;     //used for AliFlowEventSimple (to set the centrality)
-//const Int_t multmax = 1000000;     //used for AliFlowEventSimple (to set the centrality)
-
-
-//----------For RP selection----------
-//KINEMATICS (on generated and reconstructed tracks)
-Bool_t UseKineforRP =  kTRUE;
-const Double_t ptminRP = 0.0;
-const Double_t ptmaxRP = 10.0;
-const Double_t etaminRP  = -0.9;
-const Double_t etamaxRP  = 0.9;
-const Int_t    chargeRP = 1;  //not used
-
-//PID (on generated and reconstructed tracks)
-Bool_t UsePIDforRP = kFALSE;
-const Int_t PdgRP = 211;
-
-//TRACK QUALITY (on reconstructed tracks only)
-//see /CORRFW/AliCFTrackQualityCuts class
-Bool_t UseTrackQualityforRP =  kTRUE;
-const Int_t    minClustersTpcRP = 80;           //default = -1; 
-const Double_t maxChi2PerClusterTpcRP = 3.5;    //default = 1.e+09;
-const UShort_t minDedxClusterTpcRP = 0;
-const Int_t    minClustersItsRP = 2;            //panos
-const Double_t maxChi2PerClusterItsRP = 1.e+09; 
-const Int_t    minClustersTrdRP = -1;
-const Int_t    minTrackletTrdRP = -1;
-const Int_t    minTrackletTrdPidRP = -1;
-const Double_t maxChi2PerClusterTrdRP = 1.e+09;
-const ULong_t  statusRP = AliESDtrack::kTPCrefit;   //AliESDtrack::kTPCrefit &  AliESDtrack::kITSrefit 
-
-//PRIMARY (on reconstructed tracks only)
-//see /CORRFW/AliCFTrackIsPrimaryCuts class
-Bool_t UsePrimariesforRP = kTRUE;
-const Bool_t   spdVertexRP = kFALSE;
-const Bool_t   tpcVertexRP = kFALSE;
-const Float_t  minDcaToVertexXyRP = 0.;
-const Float_t  minDcaToVertexZRP = 0.;
-const Float_t  maxDcaToVertexXyRP = 2.4;         //default = 1.e+10;  //2.4;
-const Float_t  maxDcaToVertexZRP = 3.2;          //default = 1.e+10;  //3.2;
-const Bool_t   dcaToVertex2dRP = kFALSE;         //default = kFALSE;
-const Bool_t   absDcaToVertexRP = kTRUE;         //default = kTRUE;
-const Double_t minNSigmaToVertexRP = 0.;
-const Double_t maxNSigmaToVertexRP = 1.e+10; //3.; //1.e+10
-const Double_t maxSigmaDcaXySP = 1.e+10;
-const Double_t maxSigmaDcaZSP = 1.e+10;
-const Bool_t   requireSigmaToVertexSP = kFALSE;
-const Bool_t   acceptKinkDaughtersSP = kFALSE;  //default = kTRUE;
-
-//ACCEPTANCE (on generated tracks only : AliMCParticle)
-//see /CORRFW/AliCFAcceptanceCuts class
-Bool_t UseAcceptanceforRP =  kFALSE; 
-const Int_t  minTrackrefsItsRP = 0;//3;
-const Int_t  minTrackrefsTpcRP = 0;//2;
-const Int_t  minTrackrefsTrdRP = 0; 
-const Int_t  minTrackrefsTofRP = 0; 
-const Int_t  minTrackrefsMuonRP = 0; 
-//default for all is 0
-
-//----------For POI selection----------
-//KINEMATICS (on generated and reconstructed tracks)
-Bool_t UseKineforPOI = kTRUE;
-const Double_t ptminPOI = 0.0;
-const Double_t ptmaxPOI = 10.0;
-const Double_t etaminPOI  = -0.9;
-const Double_t etamaxPOI  = 0.9;
-const Int_t    chargePOI = 1;  //not used
-
-//PID (on generated and reconstructed tracks)
-Bool_t UsePIDforPOI = kFALSE;
-const Int_t PdgPOI = 321;
-
-//TRACK QUALITY (on reconstructed tracks only)
-//see /CORRFW/AliCFTrackQualityCuts class
-Bool_t UseTrackQualityforPOI = kTRUE;
-const Int_t    minClustersTpcPOI = 80;
-const Double_t maxChi2PerClusterTpcPOI = 3.5;    
-const UShort_t minDedxClusterTpcPOI = 0;
-const Int_t    minClustersItsPOI = 2;
-const Double_t maxChi2PerClusterItsPOI = 1.e+09;
-const Int_t    minClustersTrdPOI = -1;
-const Int_t    minTrackletTrdPOI = -1;
-const Int_t    minTrackletTrdPidPOI = -1;
-const Double_t maxChi2PerClusterTrdPOI = 1.e+09;
-const ULong_t  statusPOI = AliESDtrack::kTPCrefit;   
-
-//PRIMARY (on reconstructed tracks only)
-//see /CORRFW/AliCFTrackIsPrimaryCuts class
-Bool_t UsePrimariesforPOI = kTRUE;
-const Bool_t   spdVertexPOI = kFALSE;
-const Bool_t   tpcVertexPOI = kFALSE;
-const Float_t  minDcaToVertexXyPOI = 0.;
-const Float_t  minDcaToVertexZPOI = 0.;
-const Float_t  maxDcaToVertexXyPOI = 2.4;
-const Float_t  maxDcaToVertexZPOI = 3.2;
-const Bool_t   dcaToVertex2dPOI =  kFALSE;
-const Bool_t   absDcaToVertexPOI = kTRUE;
-const Double_t minNSigmaToVertexPOI = 0.;
-const Double_t maxNSigmaToVertexPOI = 1.e+10;  
-const Double_t maxSigmaDcaXyPOI = 1.e+10;
-const Double_t maxSigmaDcaZPOI = 1.e+10;
-const Bool_t   requireSigmaToVertexPOI = kFALSE;
-const Bool_t   acceptKinkDaughtersPOI = kFALSE;
-
-//ACCEPTANCE (on generated tracks only : AliMCParticle)
-//see /CORRFW/AliCFAcceptanceCuts class
-Bool_t UseAcceptanceforPOI = kFALSE;
-const Int_t minTrackrefsItsPOI = 3;
-const Int_t minTrackrefsTpcPOI = 2;
-const Int_t minTrackrefsTrdPOI = 0; 
-const Int_t minTrackrefsTofPOI = 0; 
-const Int_t minTrackrefsMuonPOI = 0; 
-
-
-//----------For Adding Flow to the Event----------
-const Bool_t AddToEvent = kFALSE;
-Double_t ellipticFlow = 0.05;
-
-
-AliAnalysisTaskFlowEvent* AddTaskFlow(TString type, Bool_t* METHODS, Bool_t QA, Bool_t* WEIGHTS)
+void AddTaskFlow(Float_t centrMin=-1,
+                Float_t centrMax=-1,
+                TString fileNameBase="AnalysisResults" )
 {
-  //boleans for the methods
-  Bool_t SP       = METHODS[0];
-  Bool_t LYZ1SUM  = METHODS[1];
-  Bool_t LYZ1PROD = METHODS[2];
-  Bool_t LYZ2SUM  = METHODS[3];
-  Bool_t LYZ2PROD = METHODS[4];
-  Bool_t LYZEP    = METHODS[5];
-  Bool_t GFC      = METHODS[6];
-  Bool_t QC       = METHODS[7];
-  Bool_t FQD      = METHODS[8];
-  Bool_t MCEP     = METHODS[9];   
-  //for using weights
+  // Define the range for eta subevents (for SP method)
+  Double_t minA = -0.9;
+  Double_t maxA = -0.5;
+  Double_t minB = 0.5;
+  Double_t maxB = 0.9;
+  /*
+  //PMD As RP
+  Double_t minB = 2.3;
+  Double_t maxB = 3.1;
+  Double_t minA = 3.1;
+  Double_t maxA = 3.9;
+  */
+
+  // AFTERBURNER
+  Bool_t useAfterBurner=kFALSE;
+  Double_t v1=0.0;
+  Double_t v2=0.0;
+  Double_t v3=0.0;
+  Double_t v4=0.0;
+  Int_t numberOfTrackClones=0; //non-flow
+
+  // Define a range of the detector to exclude
+  Bool_t ExcludeRegion = kFALSE;
+  Double_t excludeEtaMin = -0.;
+  Double_t excludeEtaMax = 0.;
+  Double_t excludePhiMin = 0.;
+  Double_t excludePhiMax = 0.;
+
+  // use physics selection class
+  Bool_t  UsePhysicsSelection = kTRUE;
+
+  // charge of poi
+  const Int_t    chargePOI = 1; 
+
+  // QA
+  Bool_t runQAtask=kFALSE;
+  Bool_t FillQAntuple=kFALSE;
+  Bool_t DoQAcorrelations=kFALSE;
+
+  // RUN SETTINGS
+  // Flow analysis method can be:(set to kTRUE or kFALSE)
+  Bool_t MCEP     = kTRUE;  // correlation with Monte Carlo reaction plane
+  Bool_t SP       = kTRUE;  // scalar product method (similar to eventplane method)
+  Bool_t GFC      = kTRUE;  // cumulants based on generating function
+  Bool_t QC       = kTRUE;  // cumulants using Q vectors
+  Bool_t FQD      = kTRUE;  // fit of the distribution of the Q vector (only integrated v)
+  Bool_t LYZ1SUM  = kTRUE;  // Lee Yang Zeroes using sum generating function (integrated v)
+  Bool_t LYZ1PROD = kFALSE;  // Lee Yang Zeroes using product generating function (integrated v)
+  Bool_t LYZ2SUM  = kFALSE; // Lee Yang Zeroes using sum generating function (second pass differential v)
+  Bool_t LYZ2PROD = kFALSE; // Lee Yang Zeroes using product generating function (second pass differential v)
+  Bool_t LYZEP    = kFALSE; // Lee Yang Zeroes Event plane using sum generating function (gives eventplane + weight)
+  Bool_t MH       = kTRUE;  // azimuthal correlators in mixed harmonics  
+  Bool_t NL       = kFALSE;  // nested loops (for instance distribution of phi1-phi2 for all distinct pairs)
+
+  Bool_t METHODS[] = {SP,LYZ1SUM,LYZ1PROD,LYZ2SUM,LYZ2PROD,LYZEP,GFC,QC,FQD,MCEP,MH,NL};
+
+  // Boolean to use/not use weights for the Q vector
+  Bool_t WEIGHTS[] = {kFALSE,kFALSE,kFALSE}; //Phi, v'(pt), v'(eta)
+
+  // SETTING THE CUTS
+
+  //---------Data selection----------
+  //kMC, kGlobal, kTPCstandalone, kSPDtracklet, kPMD
+  AliFlowTrackCuts::trackParameterType rptype = AliFlowTrackCuts::kGlobal;
+  AliFlowTrackCuts::trackParameterType poitype = AliFlowTrackCuts::kGlobal;
+
+  //---------Parameter mixing--------
+  //kPure - no mixing, kTrackWithMCkine, kTrackWithMCPID, kTrackWithMCpt
+  AliFlowTrackCuts::trackParameterMix rpmix = AliFlowTrackCuts::kPure;
+  AliFlowTrackCuts::trackParameterMix poimix = AliFlowTrackCuts::kPure;
+
+
+  const char* rptypestr = AliFlowTrackCuts::GetParamTypeName(rptype);
+  const char* poitypestr = AliFlowTrackCuts::GetParamTypeName(poitype);
+
+
+  TString centralityName("");
+  if((centrMin > 0)||(centrMax > 0)) {
+    centralityName+=Form("%.0f",centrMin);
+    centralityName+="-";
+    centralityName+=Form("%.0f",centrMax);
+  }
+
+  TString fileName(fileNameBase);
+  fileName.Append(".root");
+  //===========================================================================
+  printf("CREATE CUTS\n");
+  cout << "Used for RP: "<< rptypestr << endl;  
+  cout << "Used for POI: "<< poitypestr << endl;  
+  // EVENTS CUTS:
+  AliFlowEventCuts* cutsEvent = new AliFlowEventCuts("event cuts");
+  if((centrMin > 0)||(centrMax > 0)) {
+    cutsEvent->SetCentralityPercentileRange(centrMin,centrMax);
+    cutsEvent->SetCentralityPercentileMethod(AliFlowEventCuts::kV0);
+    cutsEvent->SetRefMultMethod(AliFlowEventCuts::kV0);
+    //cutsEvent->SetCentralityPercentileMethod(AliFlowEventCuts::kSPD1tracklets);
+    cutsEvent->SetNContributorsRange(2);
+    cutsEvent->SetPrimaryVertexZrange(-10.,10.);
+    cutsEvent->SetCutSPDvertexerAnomaly(); //"Francesco's cut"
+    cutsEvent->SetCutZDCtiming();
+  }
+  
+  // RP TRACK CUTS:
+  AliFlowTrackCuts* cutsRP = new AliFlowTrackCuts("rp cuts");
+  cutsRP->SetParamType(rptype);
+  cutsRP->SetParamMix(rpmix);
+  cutsRP->SetPtRange(0.2,5.);
+  cutsRP->SetEtaRange(-0.8,0.8);
+  //cutsRP->SetRequireCharge(kTRUE);
+  //cutsRP->SetCharge(chargeRP);
+  //cutsRP->SetPID(PdgRP);
+  if(rptype != AliFlowTrackCuts::kMC) {
+    cutsRP->SetMinNClustersTPC(70);
+    cutsRP->SetMinChi2PerClusterTPC(0.1);
+    cutsRP->SetMaxChi2PerClusterTPC(4.0);
+    cutsRP->SetMinNClustersITS(2);
+    cutsRP->SetRequireITSRefit(kTRUE);
+    cutsRP->SetRequireTPCRefit(kTRUE);
+    //cutsRP->SetMaxChi2PerClusterITS(1.e+09);
+    cutsRP->SetMaxDCAToVertexXY(0.3);
+    cutsRP->SetMaxDCAToVertexZ(0.3);
+    //cutsRP->SetDCAToVertex2D(kTRUE);
+    //cutsRP->SetMaxNsigmaToVertex(1.e+10);
+    //cutsRP->SetRequireSigmaToVertex(kFALSE);
+    cutsRP->SetAcceptKinkDaughters(kFALSE);
+  }
+  if(rptype == AliFlowTrackCuts::kPMD) {
+    cutsRP->SetEtaRange(2.3,3.9);
+    cutsRP->SetPmdDetPlane(0);
+    cutsRP->SetPmdAdc(270);
+    cutsRP->SetPmdNcell(1);
+  }
+
+
+  // POI TRACK CUTS:
+  AliFlowTrackCuts* cutsPOI = new AliFlowTrackCuts("poi cuts");
+  cutsPOI->SetParamType(poitype);
+  cutsPOI->SetParamMix(poimix);
+  cutsPOI->SetPtRange(0.0,10.);
+  cutsPOI->SetEtaRange(-1.2,1.2);
+  if(chargePOI != 0)
+    cutsPOI->SetCharge(chargePOI);
+  //cutsPOI->SetRequireCharge(kTRUE);
+  //cutsPOI->SetPID(PdgRP);
+  if(rptype != AliFlowTrackCuts::kMC) {
+    cutsPOI->SetMinNClustersTPC(80);
+    cutsPOI->SetMinChi2PerClusterTPC(0.1);
+    cutsPOI->SetMaxChi2PerClusterTPC(4.0);
+    cutsPOI->SetRequireITSRefit(kTRUE);
+    cutsPOI->SetRequireTPCRefit(kTRUE);
+    cutsPOI->SetMinNClustersITS(2);
+    //cutsPOI->SetMaxChi2PerClusterITS(1.e+09);
+    cutsPOI->SetMaxDCAToVertexXY(0.3);
+    cutsPOI->SetMaxDCAToVertexZ(0.3);
+    //cutsPOI->SetDCAToVertex2D(kTRUE);
+    //cutsPOI->SetMaxNsigmaToVertex(1.e+10);
+    //cutsPOI->SetRequireSigmaToVertex(kFALSE);
+    cutsPOI->SetAcceptKinkDaughters(kFALSE);
+    //cutsPOI->SetPID(AliPID::kProton, AliFlowTrackCuts::kTOFpid);
+    //cutsPOI->SetPID(AliPID::kPion, AliFlowTrackCuts::kTPCpid);
+    //cutsPOI->SetPID(AliPID::kProton, AliFlowTrackCuts::kTPCTOFpid);
+    //cutsPOI->SetTPCTOFpidCrossOverPt(0.4);
+    //francesco's TPC Bethe Bloch for data:
+    //cutsPOI->GetESDpid().GetTPCResponse().SetBetheBlochParameters(4.36414e-02,1.75977e+01,1.14385e-08,2.27907e+00,3.36699e+00);
+    //cutsPOI->GetESDpid().GetTPCResponse().SetMip(49);
+  }
+  if(poitype == AliFlowTrackCuts::kPMD) {
+    cutsPOI->SetEtaRange(2.3,3.9);
+    cutsPOI->SetPmdDetPlane(0);
+    cutsPOI->SetPmdAdc(270);
+    cutsPOI->SetPmdNcell(1);
+  }
+
+
   Bool_t useWeights  = WEIGHTS[0] || WEIGHTS[1] || WEIGHTS[2];
   if (useWeights) cout<<"Weights are used"<<endl;
   else cout<<"Weights are not used"<<endl;
-
-
+  
   // Get the pointer to the existing analysis manager via the static access method.
   //==============================================================================
   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
   if (!mgr) {
     Error("AddTaskFlowEvent", "No analysis manager to connect to.");
     return NULL;
-  }   
+  }
   
   // Check the analysis type using the event handlers connected to the analysis
-  // manager. The availability of MC handler cann also be checked here.
+  // manager. The availability of MC handler can also be checked here.
   //==============================================================================
   if (!mgr->GetInputEventHandler()) {
     ::Error("AddTaskFlowEvent", "This task requires an input event handler");
     return NULL;
   }  
-    
+
   // Open external input files
   //===========================================================================
   //weights: 
@@ -209,7 +223,7 @@ AliAnalysisTaskFlowEvent* AddTaskFlow(TString type, Bool_t* METHODS, Bool_t QA,
       break;
     } 
   }
-    
+  
   //LYZ2
   if (LYZ2SUM || LYZ2PROD) {
     //read the outputfile of the first run
@@ -222,14 +236,12 @@ AliAnalysisTaskFlowEvent* AddTaskFlow(TString type, Bool_t* METHODS, Bool_t QA,
       cout<<"WARNING: You do not have an output file:"<<endl;
       cout<<"         "<<pwd.Data()<<endl;
       exit(0);
-    } else {
-      outputFile = TFile::Open(pwd.Data(),"READ");
-    }
+    } else { outputFile = TFile::Open(pwd.Data(),"READ");}
     
     if (LYZ2SUM){  
       // read the output directory from LYZ1SUM 
       TString inputFileNameLYZ2SUM = "outputLYZ1SUManalysis" ;
-      inputFileNameLYZ2SUM += type;
+      inputFileNameLYZ2SUM += rptypestr;
       cout<<"The input directory is "<<inputFileNameLYZ2SUM.Data()<<endl;
       TFile* fInputFileLYZ2SUM = (TFile*)outputFile->FindObjectAny(inputFileNameLYZ2SUM.Data());
       if(!fInputFileLYZ2SUM || fInputFileLYZ2SUM->IsZombie()) { 
@@ -246,7 +258,7 @@ AliAnalysisTaskFlowEvent* AddTaskFlow(TString type, Bool_t* METHODS, Bool_t QA,
     if (LYZ2PROD){  
       // read the output directory from LYZ1PROD 
       TString inputFileNameLYZ2PROD = "outputLYZ1PRODanalysis" ;
-      inputFileNameLYZ2PROD += type;
+      inputFileNameLYZ2PROD += rptypestr;
       cout<<"The input directory is "<<inputFileNameLYZ2PROD.Data()<<endl;
       TFile* fInputFileLYZ2PROD = (TFile*)outputFile->FindObjectAny(inputFileNameLYZ2PROD.Data());
       if(!fInputFileLYZ2PROD || fInputFileLYZ2PROD->IsZombie()) { 
@@ -261,7 +273,6 @@ AliAnalysisTaskFlowEvent* AddTaskFlow(TString type, Bool_t* METHODS, Bool_t QA,
     }
   }
 
-
   if (LYZEP) {
     //read the outputfile of the second run
     TString outputFileName = "AnalysisResults2.root";
@@ -276,10 +287,10 @@ AliAnalysisTaskFlowEvent* AddTaskFlow(TString type, Bool_t* METHODS, Bool_t QA,
     } else {
       outputFile = TFile::Open(pwd.Data(),"READ");
     }
-
+    
     // read the output file from LYZ2SUM
     TString inputFileNameLYZEP = "outputLYZ2SUManalysis" ;
-    inputFileNameLYZEP += type;
+    inputFileNameLYZEP += rptypestr;
     cout<<"The input file is "<<inputFileNameLYZEP.Data()<<endl;
     TFile* fInputFileLYZEP = (TFile*)outputFile->FindObjectAny(inputFileNameLYZEP.Data());
     if(!fInputFileLYZEP || fInputFileLYZEP->IsZombie()) { 
@@ -293,360 +304,61 @@ AliAnalysisTaskFlowEvent* AddTaskFlow(TString type, Bool_t* METHODS, Bool_t QA,
     cout<<"LYZEP input file/list read..."<<endl;
   }
   
-
-
-  // Create the task, add it to the manager.
-  //===========================================================================
-  AliAnalysisTaskFlowEvent *taskFE = NULL;
-  if (QA) { 
-    if(AddToEvent) { 
-      taskFE = new AliAnalysisTaskFlowEvent("TaskFlowEvent",kTRUE,1);
-      taskFE->SetEllipticFlowValue(ellipticFlow); }    //TEST
-    else {taskFE = new AliAnalysisTaskFlowEvent("TaskFlowEvent",kTRUE); }
-    taskFE->SetAnalysisType(type);
-    if (UseMultCut) {
-      taskFE->SetMinMult(multmin);
-      taskFE->SetMaxMult(multmax);
-    }
-    taskFE->SetSubeventEtaRange(minA, maxA, minB, maxB);
-    if (UsePhysicsSelection) {
-      taskFE->SelectCollisionCandidates();
-      cout<<"Using Physics Selection"<<endl;
-    }
-    mgr->AddTask(taskFE);
-  }
-  else { 
-    if(AddToEvent) { 
-      taskFE = new AliAnalysisTaskFlowEvent("TaskFlowEvent",kFALSE,1);
-      taskFE->SetEllipticFlowValue(ellipticFlow); }    //TEST
-    else {taskFE = new AliAnalysisTaskFlowEvent("TaskFlowEvent",kFALSE); }
-    taskFE->SetAnalysisType(type);
-    if (UseMultCut) {
-      taskFE->SetMinMult(multmin);
-      taskFE->SetMaxMult(multmax);
-    }
-    taskFE->SetSubeventEtaRange(minA, maxA, minB, maxB);
-    if (UsePhysicsSelection) {
-      taskFE->SelectCollisionCandidates();
-      cout<<"Using Physics Selection"<<endl;
-    }
-    mgr->AddTask(taskFE);
-  }
-  // Create cuts using the correction framework (CORRFW)
-  //===========================================================================
-  if (QA){
-    //Set TList for the QA histograms
-    TList* qaRP  = new TList(); 
-    TList* qaPOI = new TList();
-  }
-
-  //----------Event cuts----------
-  AliCFEventGenCuts* mcEventCuts = new AliCFEventGenCuts("mcEventCuts","MC-level event cuts");
-  mcEventCuts->SetNTracksCut(multminESD,multmaxESD); 
-  mcEventCuts->SetRequireVtxCuts(requireVtxCuts);
-  mcEventCuts->SetVertexXCut(vertexXmin, vertexXmax);
-  mcEventCuts->SetVertexYCut(vertexYmin, vertexYmax);
-  mcEventCuts->SetVertexZCut(vertexZmin, vertexZmax);
-  if (QA) { 
-    mcEventCuts->SetQAOn(qaRP);
-  }
-  AliCFEventRecCuts* recEventCuts = new AliCFEventRecCuts("recEventCuts","rec-level event cuts");
-  recEventCuts->SetNTracksCut(multminESD,multmaxESD); 
-  recEventCuts->SetRequireVtxCuts(requireVtxCuts);
-  recEventCuts->SetVertexXCut(vertexXmin, vertexXmax);
-  recEventCuts->SetVertexYCut(vertexYmin, vertexYmax);
-  recEventCuts->SetVertexZCut(vertexZmin, vertexZmax);
-  if (QA) { 
-    recEventCuts->SetQAOn(qaRP);
-  }
-  
-  //----------Cuts for RP----------
-  //KINEMATICS (MC and reconstructed)
-  AliCFTrackKineCuts* mcKineCutsRP = new AliCFTrackKineCuts("mcKineCutsRP","MC-level kinematic cuts");
-  mcKineCutsRP->SetPtRange(ptminRP,ptmaxRP);
-  mcKineCutsRP->SetEtaRange(etaminRP,etamaxRP);
-  //mcKineCutsRP->SetChargeMC(chargeRP);
-  if (QA) { 
-    mcKineCutsRP->SetQAOn(qaRP);
-  }
-
-  AliCFTrackKineCuts *recKineCutsRP = new AliCFTrackKineCuts("recKineCutsRP","rec-level kine cuts");
-  recKineCutsRP->SetPtRange(ptminRP,ptmaxRP);
-  recKineCutsRP->SetEtaRange(etaminRP,etamaxRP);
-  //recKineCutsRP->SetChargeRec(chargeRP);
-  if (QA) { 
-    recKineCutsRP->SetQAOn(qaRP);
-  }
-
-  //PID (MC and reconstructed)
-  AliCFParticleGenCuts* mcGenCutsRP = new AliCFParticleGenCuts("mcGenCutsRP","MC particle generation cuts for RP");
-  mcGenCutsRP->SetRequireIsPrimary();
-  if (UsePIDforRP) {mcGenCutsRP->SetRequirePdgCode(PdgRP);}
-  if (QA) { 
-    mcGenCutsRP->SetQAOn(qaRP);
-  }
-
-  int n_species = AliPID::kSPECIES ;
-  Double_t* prior = new Double_t[n_species];
   
-  prior[0] = 0.0244519 ;
-  prior[1] = 0.0143988 ;
-  prior[2] = 0.805747  ;
-  prior[3] = 0.0928785 ;
-  prior[4] = 0.0625243 ;
-  
-  AliCFTrackCutPid* cutPidRP = NULL;
-  if(UsePIDforRP) {
-    cutPidRP = new AliCFTrackCutPid("cutPidRP","ESD_PID for RP") ;
-    cutPidRP->SetPriors(prior);
-    cutPidRP->SetProbabilityCut(0.0);
-    cutPidRP->SetDetectors("TPC TOF");
-    switch(TMath::Abs(PDG1)) {
-    case 11   : cutPidRP->SetParticleType(AliPID::kElectron, kTRUE); break;
-    case 13   : cutPidRP->SetParticleType(AliPID::kMuon    , kTRUE); break;
-    case 211  : cutPidRP->SetParticleType(AliPID::kPion    , kTRUE); break;
-    case 321  : cutPidRP->SetParticleType(AliPID::kKaon    , kTRUE); break;
-    case 2212 : cutPidRP->SetParticleType(AliPID::kProton  , kTRUE); break;
-    default   : printf("UNDEFINED PID\n"); break;
-    }
-    if (QA) { 
-      cutPidRP->SetQAOn(qaRP); 
+  // Create the FMD task and add it to the manager
+  //===========================================================================
+  if (rptypestr == "FMD") {
+    AliFMDAnalysisTaskSE *taskfmd = NULL;
+    if (rptypestr == "FMD") {
+      taskfmd = new AliFMDAnalysisTaskSE("TaskFMD");
+      mgr->AddTask(taskfmd);
+      
+      AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+      pars->Init();
+      pars->SetProcessPrimary(kTRUE); //for MC only
+      pars->SetProcessHits(kFALSE);
+      
+      //pars->SetRealData(kTRUE); //for real data
+      //pars->SetProcessPrimary(kFALSE); //for real data
     }
   }
   
-  //TRACK QUALITY
-  AliCFTrackQualityCuts *recQualityCutsRP = new AliCFTrackQualityCuts("recQualityCutsRP","rec-level quality cuts");
-  recQualityCutsRP->SetMinNClusterTPC(minClustersTpcRP);
-  //recQualityCutsRP->SetMinFoundClusterTPC(minFoundClustersTpcRP); //only for internal TPC QA
-  recQualityCutsRP->SetMaxChi2PerClusterTPC(maxChi2PerClusterTpcRP);
-  recQualityCutsRP->SetMinNdEdxClusterTPC(minDedxClusterTpcRP);     //to reject secondaries
-
-  recQualityCutsRP->SetMinNClusterITS(minClustersItsRP);
-  recQualityCutsRP->SetMaxChi2PerClusterITS(maxChi2PerClusterItsRP);
-
-  recQualityCutsRP->SetMinNClusterTRD(minClustersTrdRP);
-  recQualityCutsRP->SetMinNTrackletTRD(minTrackletTrdRP);
-  recQualityCutsRP->SetMinNTrackletTRDpid(minTrackletTrdPidRP);
-  recQualityCutsRP->SetMaxChi2PerTrackletTRD(maxChi2PerClusterTrdRP);
-  recQualityCutsRP->SetStatus(statusRP);  
-  if (QA) { 
-    recQualityCutsRP->SetQAOn(qaRP);
-  }
-
-  /* 
-  //How to set this?
-  void SetMaxCovDiagonalElements(Float_t c1=1.e+09, Float_t c2=1.e+09, Float_t c3=1.e+09, Float_t c4=1.e+09, Float_t c5=1.e+09)
-    {fCovariance11Max=c1;fCovariance22Max=c2;fCovariance33Max=c3;fCovariance44Max=c4;fCovariance55Max=c5;}
-  */
-
-  //PRIMARIES
-  AliCFTrackIsPrimaryCuts *recIsPrimaryCutsRP = new AliCFTrackIsPrimaryCuts("recIsPrimaryCutsRP","rec-level isPrimary cuts");
-  recIsPrimaryCutsRP->UseSPDvertex(spdVertexRP);
-  recIsPrimaryCutsRP->UseTPCvertex(tpcVertexRP);
-  recIsPrimaryCutsRP->SetMinDCAToVertexXY(minDcaToVertexXyRP); 
-  recIsPrimaryCutsRP->SetMinDCAToVertexZ(minDcaToVertexZRP);
-  recIsPrimaryCutsRP->SetMaxDCAToVertexXY(maxDcaToVertexXyRP);
-  recIsPrimaryCutsRP->SetMaxDCAToVertexZ(maxDcaToVertexZRP); 
-  recIsPrimaryCutsRP->SetDCAToVertex2D(dcaToVertex2dRP);
-  recIsPrimaryCutsRP->SetAbsDCAToVertex(absDcaToVertexRP);
-  recIsPrimaryCutsRP->SetMinNSigmaToVertex(minNSigmaToVertexRP); 
-  recIsPrimaryCutsRP->SetMaxNSigmaToVertex(maxNSigmaToVertexRP); 
-  recIsPrimaryCutsRP->SetMaxSigmaDCAxy(maxSigmaDcaXySP);
-  recIsPrimaryCutsRP->SetMaxSigmaDCAz(maxSigmaDcaZSP);
-  recIsPrimaryCutsRP->SetRequireSigmaToVertex(requireSigmaToVertexSP);
-  recIsPrimaryCutsRP->SetAcceptKinkDaughters(acceptKinkDaughtersSP);
-  if (QA) { 
-    recIsPrimaryCutsRP->SetQAOn(qaRP);
-  }
-  
-  //ACCEPTANCE
-  AliCFAcceptanceCuts *mcAccCutsRP = new AliCFAcceptanceCuts("mcAccCutsRP","MC acceptance cuts");
-  mcAccCutsRP->SetMinNHitITS(minTrackrefsItsRP);
-  mcAccCutsRP->SetMinNHitTPC(minTrackrefsTpcRP);
-  mcAccCutsRP->SetMinNHitTRD(minTrackrefsTrdRP); 
-  mcAccCutsRP->SetMinNHitTOF(minTrackrefsTofRP);
-  mcAccCutsRP->SetMinNHitMUON(minTrackrefsMuonRP);
-  if (QA) { 
-    mcAccCutsRP->SetQAOn(qaRP);
-  }
+  // Create the task, add it to the manager.
+  //===========================================================================
+  AliAnalysisTaskFlowEvent *taskFE = NULL;
 
-  
-  //----------Cuts for POI----------
-  //KINEMATICS (MC and reconstructed)
-  AliCFTrackKineCuts* mcKineCutsPOI = new AliCFTrackKineCuts("mcKineCutsPOI","MC-level kinematic cuts");
-  mcKineCutsPOI->SetPtRange(ptminPOI,ptmaxPOI);
-  mcKineCutsPOI->SetEtaRange(etaminPOI,etamaxPOI);
-  //mcKineCutsPOI->SetChargeMC(chargePOI);
-  if (QA) { 
-    mcKineCutsPOI->SetQAOn(qaPOI);
-  }
-  
-  AliCFTrackKineCuts *recKineCutsPOI = new AliCFTrackKineCuts("recKineCutsPOI","rec-level kine cuts");
-  recKineCutsPOI->SetPtRange(ptminPOI,ptmaxPOI);
-  recKineCutsPOI->SetEtaRange(etaminPOI,etamaxPOI);
-  //recKineCutsPOI->SetChargeRec(chargePOI);
-  if (QA) { 
-    recKineCutsPOI->SetQAOn(qaPOI);
-  }
-  
-  //PID (MC and reconstructed)
-  AliCFParticleGenCuts* mcGenCutsPOI = new AliCFParticleGenCuts("mcGenCutsPOI","MC particle generation cuts for POI");
-  mcGenCutsPOI->SetRequireIsPrimary();
-  if (UsePIDforPOI) {mcGenCutsPOI->SetRequirePdgCode(PdgPOI);}
-  if (QA) { 
-    mcGenCutsPOI->SetQAOn(qaPOI);
-  }
-
-  AliCFTrackCutPid* cutPidPOI = NULL;
-  if (UsePIDforPOI) {
-    cutPidPOI = new AliCFTrackCutPid("cutPidPOI","ESD_PID for POI") ;
-    cutPidPOI->SetPriors(prior);
-    cutPidPOI->SetProbabilityCut(0.0);
-    cutPidPOI->SetDetectors("TPC TOF");
-    switch(TMath::Abs(PDG2)) {
-    case 11   : cutPidPOI->SetParticleType(AliPID::kElectron, kTRUE); break;
-    case 13   : cutPidPOI->SetParticleType(AliPID::kMuon    , kTRUE); break;
-    case 211  : cutPidPOI->SetParticleType(AliPID::kPion    , kTRUE); break;
-    case 321  : cutPidPOI->SetParticleType(AliPID::kKaon    , kTRUE); break;
-    case 2212 : cutPidPOI->SetParticleType(AliPID::kProton  , kTRUE); break;
-    default   : printf("UNDEFINED PID\n"); break;
-    }
-    if (QA) { 
-      cutPidPOI->SetQAOn(qaPOI);
+  if(useAfterBurner)
+    { 
+      taskFE = new AliAnalysisTaskFlowEvent("TaskFlowEvent","",kFALSE,1);
+      taskFE->SetFlow(v1,v2,v3,v4); 
+      taskFE->SetNonFlowNumberOfTrackClones(numberOfTrackClones);
+      taskFE->SetAfterburnerOn();
     }
+  else {taskFE = new AliAnalysisTaskFlowEvent("TaskFlowEvent","",kFALSE); }
+  if (ExcludeRegion) {
+    taskFE->DefineDeadZone(excludeEtaMin, excludeEtaMax, excludePhiMin, excludePhiMax); 
   }
-
-  //TRACK QUALITY
-  AliCFTrackQualityCuts *recQualityCutsPOI = new AliCFTrackQualityCuts("recQualityCutsPOI","rec-level quality cuts");
-  recQualityCutsPOI->SetMinNClusterTPC(minClustersTpcPOI);
-  //recQualityCutsPOI->SetMinFoundClusterTPC(minFoundClustersTpcPOI); //only for internal TPC QA
-  recQualityCutsPOI->SetMaxChi2PerClusterTPC(maxChi2PerClusterTpcPOI);
-  recQualityCutsPOI->SetMinNdEdxClusterTPC(minDedxClusterTpcPOI);     //to reject secondaries
-
-  recQualityCutsPOI->SetMinNClusterITS(minClustersItsPOI);
-  recQualityCutsPOI->SetMaxChi2PerClusterITS(maxChi2PerClusterItsPOI);
-
-  recQualityCutsPOI->SetMinNClusterTRD(minClustersTrdPOI);
-  recQualityCutsPOI->SetMinNTrackletTRD(minTrackletTrdPOI);
-  recQualityCutsPOI->SetMinNTrackletTRDpid(minTrackletTrdPidPOI);
-  recQualityCutsPOI->SetMaxChi2PerTrackletTRD(maxChi2PerClusterTrdPOI);
-  recQualityCutsPOI->SetStatus(statusPOI); 
-  if (QA) { 
-    recQualityCutsPOI->SetQAOn(qaPOI);
-  }
-
-  //PRIMARIES
-  AliCFTrackIsPrimaryCuts *recIsPrimaryCutsPOI = new AliCFTrackIsPrimaryCuts("recIsPrimaryCutsPOI","rec-level isPrimary cuts");
-  recIsPrimaryCutsPOI->UseSPDvertex(spdVertexPOI);
-  recIsPrimaryCutsPOI->UseTPCvertex(tpcVertexPOI);
-  recIsPrimaryCutsPOI->SetMinDCAToVertexXY(minDcaToVertexXyPOI); 
-  recIsPrimaryCutsPOI->SetMinDCAToVertexZ(minDcaToVertexZPOI);
-  recIsPrimaryCutsPOI->SetMaxDCAToVertexXY(maxDcaToVertexXyPOI);
-  recIsPrimaryCutsPOI->SetMaxDCAToVertexZ(maxDcaToVertexZPOI); 
-  recIsPrimaryCutsPOI->SetDCAToVertex2D(dcaToVertex2dPOI);
-  recIsPrimaryCutsPOI->SetAbsDCAToVertex(absDcaToVertexPOI);
-  recIsPrimaryCutsPOI->SetMinNSigmaToVertex(minNSigmaToVertexPOI); 
-  recIsPrimaryCutsPOI->SetMaxNSigmaToVertex(maxNSigmaToVertexPOI); 
-  recIsPrimaryCutsPOI->SetMaxSigmaDCAxy(maxSigmaDcaXyPOI);
-  recIsPrimaryCutsPOI->SetMaxSigmaDCAz(maxSigmaDcaZPOI);
-  recIsPrimaryCutsPOI->SetRequireSigmaToVertex(requireSigmaToVertexPOI);
-  recIsPrimaryCutsPOI->SetAcceptKinkDaughters(acceptKinkDaughtersPOI);
-  if (QA) { 
-    recIsPrimaryCutsPOI->SetQAOn(qaPOI);
-  }
-
-  //ACCEPTANCE
-  AliCFAcceptanceCuts *mcAccCutsPOI = new AliCFAcceptanceCuts("mcAccCutsPOI","MC acceptance cuts");
-  mcAccCutsPOI->SetMinNHitITS(minTrackrefsItsPOI);
-  mcAccCutsPOI->SetMinNHitTPC(minTrackrefsTpcPOI);
-  mcAccCutsPOI->SetMinNHitTRD(minTrackrefsTrdPOI); 
-  mcAccCutsPOI->SetMinNHitTOF(minTrackrefsTofPOI);
-  mcAccCutsPOI->SetMinNHitMUON(minTrackrefsMuonPOI);
-  if (QA) { 
-    mcAccCutsPOI->SetQAOn(qaPOI);
-  }
-
-     
-  
-  //----------Create Cut Lists----------
-  printf("CREATE EVENT CUTS\n");
-  TObjArray* mcEventList = new TObjArray(0);  
-  if (UseMultCutforESD) mcEventList->AddLast(mcEventCuts);//cut on mult and vertex
-    
-  TObjArray* recEventList = new TObjArray(0);
-  if (UseMultCutforESD) recEventList->AddLast(recEventCuts);//cut on mult and vertex
-
-  printf("CREATE MC KINE CUTS\n");
-  TObjArray* mcListRP = new TObjArray(0);
-  if (UseKineforRP) mcListRP->AddLast(mcKineCutsRP); //cut on pt/eta/phi
-  mcListRP->AddLast(mcGenCutsRP); //cut on primary and if (UsePIDforRP) MC PID
-  
-  TObjArray* mcListPOI = new TObjArray(0);
-  if (UseKineforPOI) mcListPOI->AddLast(mcKineCutsPOI); //cut on pt/eta/phi
-  mcListPOI->AddLast(mcGenCutsPOI); //cut on primary and if (UsePIDforPOI) MC PID
-  
-  printf("CREATE MC ACCEPTANCE CUTS\n");
-  TObjArray* accListRP = new TObjArray(0) ;
-  if (UseAcceptanceforRP) accListRP->AddLast(mcAccCutsRP); //cut on number of track references
-  
-  TObjArray* accListPOI = new TObjArray(0) ;
-  if (UseAcceptanceforPOI) accListPOI->AddLast(mcAccCutsPOI); //cut on number of track references
-  
-  printf("CREATE ESD RECONSTRUCTION CUTS\n");
-  TObjArray* recListRP = new TObjArray(0) ;
-  if (UseKineforRP)         recListRP->AddLast(recKineCutsRP); //cut on pt/eta/phi
-  if (UseTrackQualityforRP) recListRP->AddLast(recQualityCutsRP);
-  if (UsePrimariesforRP)    recListRP->AddLast(recIsPrimaryCutsRP); //cut if it is a primary
-  
-  TObjArray* recListPOI = new TObjArray(0) ;
-  if (UseKineforPOI)         recListPOI->AddLast(recKineCutsPOI); //cut on pt/eta/phi
-  if (UseTrackQualityforPOI) recListPOI->AddLast(recQualityCutsPOI);
-  if (UsePrimariesforPOI)    recListPOI->AddLast(recIsPrimaryCutsPOI); //cut if it is a primary
-  
-  printf("CREATE ESD PID CUTS\n");
-  TObjArray* fPIDCutListRP = new TObjArray(0) ;
-  if(UsePIDforRP) {fPIDCutListRP->AddLast(cutPidRP);} //cut on ESD PID
-  
-  TObjArray* fPIDCutListPOI = new TObjArray(0) ;
-  if (UsePIDforPOI)  {fPIDCutListPOI->AddLast(cutPidPOI);} //cut on ESD PID
-  
-
-  //----------Add Cut Lists to the CF Manager----------
-  printf("CREATE INTERFACE AND CUTS\n");
-  AliCFManager* cfmgrRP = new AliCFManager();
-  cfmgrRP->SetNStepEvent(3);
-  cfmgrRP->SetEventCutsList(AliCFManager::kEvtGenCuts,mcEventList); 
-  cfmgrRP->SetEventCutsList(AliCFManager::kEvtRecCuts,recEventList); 
-  cfmgrRP->SetNStepParticle(4); 
-  cfmgrRP->SetParticleCutsList(AliCFManager::kPartGenCuts,mcListRP);
-  cfmgrRP->SetParticleCutsList(AliCFManager::kPartAccCuts,accListRP);
-  cfmgrRP->SetParticleCutsList(AliCFManager::kPartRecCuts,recListRP);
-  cfmgrRP->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutListRP);
-  
-  AliCFManager* cfmgrPOI = new AliCFManager();
-  cfmgrPOI->SetNStepEvent(3);
-  cfmgrPOI->SetEventCutsList(AliCFManager::kEvtGenCuts,mcEventList); 
-  cfmgrPOI->SetEventCutsList(AliCFManager::kEvtRecCuts,recEventList); 
-  cfmgrPOI->SetNStepParticle(4); 
-  cfmgrPOI->SetParticleCutsList(AliCFManager::kPartGenCuts,mcListPOI);
-  cfmgrPOI->SetParticleCutsList(AliCFManager::kPartAccCuts,accListPOI);
-  cfmgrPOI->SetParticleCutsList(AliCFManager::kPartRecCuts,recListPOI);
-  cfmgrPOI->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutListPOI);
-  
-  if (QA) {
-    taskFE->SetQAList1(qaRP);
-    taskFE->SetQAList2(qaPOI);
+  taskFE->SetSubeventEtaRange(minA, maxA, minB, maxB);
+  if (UsePhysicsSelection) {
+    //taskFE->SelectCollisionCandidates(AliVEvent::kUserDefined);
+    taskFE->SelectCollisionCandidates(AliVEvent::kMB);
+    cout<<"Using Physics Selection"<<endl;
   }
-  taskFE->SetCFManager1(cfmgrRP);
-  taskFE->SetCFManager2(cfmgrPOI);
-
+  mgr->AddTask(taskFE);
+  
+  // Pass cuts for RPs and POIs to the task:
+  taskFE->SetCutsEvent(cutsEvent);
 
+  // Pass cuts for RPs and POIs to the task:
+  taskFE->SetCutsRP(cutsRP);
+  taskFE->SetCutsPOI(cutsPOI);
 
   // Create the analysis tasks, add them to the manager.
   //===========================================================================
   if (SP){
     AliAnalysisTaskScalarProduct *taskSP = new AliAnalysisTaskScalarProduct("TaskScalarProduct",WEIGHTS[0]);
+    taskSP->SetRelDiffMsub(1.0);
+    taskSP->SetApplyCorrectionForNUA(kTRUE);
     mgr->AddTask(taskSP);
   }
   if (LYZ1SUM){
@@ -689,175 +401,294 @@ AliAnalysisTaskFlowEvent* AddTaskFlow(TString type, Bool_t* METHODS, Bool_t QA,
     taskQC->SetUsePhiWeights(WEIGHTS[0]); 
     taskQC->SetUsePtWeights(WEIGHTS[1]);
     taskQC->SetUseEtaWeights(WEIGHTS[2]); 
+    taskQC->SetnBinsMult(10000);
+    taskQC->SetMinMult(0.);
+    taskQC->SetMaxMult(10000.);
+    taskQC->SetApplyCorrectionForNUA(kTRUE);
     mgr->AddTask(taskQC);
   }
   if (FQD){
     AliAnalysisTaskFittingQDistribution *taskFQD = new AliAnalysisTaskFittingQDistribution("TaskFittingQDistribution",kFALSE);
     taskFQD->SetUsePhiWeights(WEIGHTS[0]); 
+    taskFQD->SetqMin(0.);
+    taskFQD->SetqMax(1000.);
+    taskFQD->SetqNbins(10000);
     mgr->AddTask(taskFQD);
   }
   if (MCEP){
     AliAnalysisTaskMCEventPlane *taskMCEP = new AliAnalysisTaskMCEventPlane("TaskMCEventPlane");
     mgr->AddTask(taskMCEP);
   }
+  if (MH){
+    AliAnalysisTaskMixedHarmonics *taskMH = new AliAnalysisTaskMixedHarmonics("TaskMixedHarmonics",useWeights);
+   taskMH->SetHarmonic(1); // n in cos[n(phi1+phi2-2phi3)] and cos[n(psi1+psi2-2phi3)]
+    taskMH->SetNoOfMultipicityBins(10000);
+    taskMH->SetMultipicityBinWidth(1.);
+    taskMH->SetMinMultiplicity(1.);
+    taskMH->SetCalculateVsM(kTRUE);
+    taskMH->SetCorrectForDetectorEffects(kTRUE);
+    taskMH->SetEvaluateDifferential3pCorrelator(kTRUE); // evaluate <<cos[n(psi1+psi2-2phi3)]>> (Remark: two nested loops)    
+    if(chargePOI == 0)
+      taskMH->SetOppositeChargesPOI(kTRUE);    
+    else
+      taskMH->SetOppositeChargesPOI(kFALSE); // POIs psi1 and psi2 in cos[n(psi1+psi2-2phi3)] will have opposite charges  
+    mgr->AddTask(taskMH);
+  }  
+  if (NL){
+    AliAnalysisTaskNestedLoops *taskNL = new AliAnalysisTaskNestedLoops("TaskNestedLoops",useWeights);
+    taskNL->SetHarmonic(1); // n in cos[n(phi1+phi2-2phi3)] and cos[n(psi1+psi2-2phi3)]
+    taskNL->SetEvaluateNestedLoopsForRAD(kTRUE); // RAD = Relative Angle Distribution
+    taskNL->SetEvaluateNestedLoopsForMH(kTRUE); // evalaute <<cos[n(phi1+phi2-2phi3)]>> (Remark: three nested loops)   
+    taskNL->SetEvaluateDifferential3pCorrelator(kFALSE); // evaluate <<cos[n(psi1+psi2-2phi3)]>>  (Remark: three nested loops)   
+    taskNL->SetOppositeChargesPOI(kFALSE); // POIs psi1 and psi2 in cos[n(psi1+psi2-2phi3)] will have opposite charges  
+    mgr->AddTask(taskNL);
+  }
 
   // Create the output container for the data produced by the task
   // Connect to the input and output containers
   //===========================================================================
   AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
-  AliAnalysisDataContainer *coutputFE = mgr->CreateContainer("cobjFlowEventSimple",  AliFlowEventSimple::Class(),AliAnalysisManager::kExchangeContainer);
+  
+  if (rptypestr == "FMD") {
+    AliAnalysisDataContainer *coutputFMD = 
+      mgr->CreateContainer(Form("BackgroundCorrected_%s",centralityName.Data()), TList::Class(), AliAnalysisManager::kExchangeContainer);
+    //input and output taskFMD     
+    mgr->ConnectInput(taskfmd, 0, cinput1);
+    mgr->ConnectOutput(taskfmd, 1, coutputFMD);
+    //input into taskFE
+    mgr->ConnectInput(taskFE,1,coutputFMD);
+  }
+  
+  AliAnalysisDataContainer *coutputFE = 
+  mgr->CreateContainer(Form("cobjFlowEventSimple_%s",centralityName.Data()),AliFlowEventSimple::Class(),AliAnalysisManager::kExchangeContainer);
   mgr->ConnectInput(taskFE,0,cinput1); 
   mgr->ConnectOutput(taskFE,1,coutputFE);
-
-  if (QA) { 
-    TString qaNameRPFE = AliAnalysisManager::GetCommonFileName();
-    qaNameRPFE += ":QAforRP_FE_";
-    qaNameRPFE += type;
-
-    AliAnalysisDataContainer *coutputQA1FE =
-      mgr->CreateContainer("QARPFE", TList::Class(),AliAnalysisManager::kOutputContainer,qaNameRPFE); 
-    
-    TString qaNamePOIFE = AliAnalysisManager::GetCommonFileName();
-    qaNamePOIFE += ":QAforPOI_FE_";
-    qaNamePOIFE += type;
-        
-    AliAnalysisDataContainer *coutputQA2FE =
-      mgr->CreateContainer("QAPOIFE", TList::Class(),AliAnalysisManager::kOutputContainer,qaNamePOIFE); 
-
-    mgr->ConnectOutput(taskFE,2,coutputQA1FE); 
-    mgr->ConnectOutput(taskFE,3,coutputQA2FE); 
-  }
+  
 
   // Create the output containers for the data produced by the analysis tasks
   // Connect to the input and output containers
   //===========================================================================
   if (useWeights) {    
-    AliAnalysisDataContainer *cinputWeights = mgr->CreateContainer("cobjWeights",TList::Class(),AliAnalysisManager::kInputContainer); 
+    AliAnalysisDataContainer *cinputWeights = mgr->CreateContainer(Form("cobjWeights_%s",centralityName.Data()),
+                                                                  TList::Class(),AliAnalysisManager::kInputContainer); 
   }
 
   if(SP) {
-    TString outputSP = AliAnalysisManager::GetCommonFileName();
+    TString outputSP = fileName;
     outputSP += ":outputSPanalysis";
-    outputSP+= type;
-    
-    AliAnalysisDataContainer *coutputSP = mgr->CreateContainer("cobjSP", TList::Class(),AliAnalysisManager::kOutputContainer,outputSP); 
+    outputSP+= rptypestr;
+    TString cobjSP = "cobjSP";
+    if((centrMin > 0)||(centrMax > 0)) {
+      cobjSP += "_"; cobjSP += centralityName.Data(); }
+    AliAnalysisDataContainer *coutputSP = mgr->CreateContainer(cobjSP.Data(), 
+                                                              TList::Class(),AliAnalysisManager::kOutputContainer,outputSP); 
     mgr->ConnectInput(taskSP,0,coutputFE); 
     mgr->ConnectOutput(taskSP,1,coutputSP); 
     if (WEIGHTS[0]) {
       mgr->ConnectInput(taskSP,1,cinputWeights);
       cinputWeights->SetData(weightsList);
-    } 
+    }
   }
   if(LYZ1SUM) {
-    TString outputLYZ1SUM = AliAnalysisManager::GetCommonFileName();
+    TString outputLYZ1SUM = fileName;
     outputLYZ1SUM += ":outputLYZ1SUManalysis";
-    outputLYZ1SUM+= type;
-    
-    AliAnalysisDataContainer *coutputLYZ1SUM = mgr->CreateContainer("cobjLYZ1SUM", TList::Class(),AliAnalysisManager::kOutputContainer,outputLYZ1SUM); 
-    mgr->ConnectInput(taskLYZ1SUM,0,coutputFE); 
-    mgr->ConnectOutput(taskLYZ1SUM,0,coutputLYZ1SUM); 
+    outputLYZ1SUM+= rptypestr;
+    TString cobjLYZ1SUM = "cobjLYZ1SUM";
+    if((centrMin > 0)||(centrMax > 0)) {
+      cobjLYZ1SUM += "_"; cobjLYZ1SUM += centralityName.Data(); }
+    AliAnalysisDataContainer *coutputLYZ1SUM = mgr->CreateContainer(cobjLYZ1SUM.Data(),
+                                                                   TList::Class(),AliAnalysisManager::kOutputContainer,outputLYZ1SUM); 
+    mgr->ConnectInput(taskLYZ1SUM,0,coutputFE);
+    mgr->ConnectOutput(taskLYZ1SUM,1,coutputLYZ1SUM);
   }
   if(LYZ1PROD) {
-    TString outputLYZ1PROD = AliAnalysisManager::GetCommonFileName();
+    TString outputLYZ1PROD = fileName;
     outputLYZ1PROD += ":outputLYZ1PRODanalysis";
-    outputLYZ1PROD+= type;
-    
-    AliAnalysisDataContainer *coutputLYZ1PROD = mgr->CreateContainer("cobjLYZ1PROD", TList::Class(),AliAnalysisManager::kOutputContainer,outputLYZ1PROD); 
+    outputLYZ1PROD+= rptypestr;
+    TString cobjLYZ1PROD = "cobjLYZ1PROD";
+    if((centrMin > 0)||(centrMax > 0)) {
+      cobjLYZ1PROD += "_"; cobjLYZ1PROD += centralityName.Data(); }
+    AliAnalysisDataContainer *coutputLYZ1PROD = mgr->CreateContainer(cobjLYZ1PROD.Data(),
+                                                                    TList::Class(),AliAnalysisManager::kOutputContainer,outputLYZ1PROD); 
     mgr->ConnectInput(taskLYZ1PROD,0,coutputFE); 
-    mgr->ConnectOutput(taskLYZ1PROD,0,coutputLYZ1PROD);
+    mgr->ConnectOutput(taskLYZ1PROD,1,coutputLYZ1PROD);
   }
   if(LYZ2SUM) {
-    AliAnalysisDataContainer *cinputLYZ2SUM = mgr->CreateContainer("cobjLYZ2SUMin",TList::Class(),AliAnalysisManager::kInputContainer);
-    TString outputLYZ2SUM = AliAnalysisManager::GetCommonFileName();
+    AliAnalysisDataContainer *cinputLYZ2SUM = mgr->CreateContainer(Form("cobjLYZ2SUMin_%s",centralityName.Data()),
+                                                                  TList::Class(),AliAnalysisManager::kInputContainer);
+    TString outputLYZ2SUM = fileName;
     outputLYZ2SUM += ":outputLYZ2SUManalysis";
-    outputLYZ2SUM+= type;
-    
-    AliAnalysisDataContainer *coutputLYZ2SUM = mgr->CreateContainer("cobjLYZ2SUM", TList::Class(),AliAnalysisManager::kOutputContainer,outputLYZ2SUM); 
+    outputLYZ2SUM+= rptypestr;
+    TString cobjLYZ2SUM = "cobjLYZ2SUM";
+    if((centrMin > 0)||(centrMax > 0)) {
+      cobjLYZ2SUM += "_"; cobjLYZ2SUM += centralityName.Data(); }
+    AliAnalysisDataContainer *coutputLYZ2SUM = mgr->CreateContainer(cobjLYZ2SUM.Data(),
+                                                                   TList::Class(),AliAnalysisManager::kOutputContainer,outputLYZ2SUM); 
     mgr->ConnectInput(taskLYZ2SUM,0,coutputFE); 
     mgr->ConnectInput(taskLYZ2SUM,1,cinputLYZ2SUM);
-    mgr->ConnectOutput(taskLYZ2SUM,0,coutputLYZ2SUM); 
+    mgr->ConnectOutput(taskLYZ2SUM,1,coutputLYZ2SUM); 
     cinputLYZ2SUM->SetData(fInputListLYZ2SUM);
   }
   if(LYZ2PROD) {
-    AliAnalysisDataContainer *cinputLYZ2PROD = mgr->CreateContainer("cobjLYZ2PRODin",TList::Class(),AliAnalysisManager::kInputContainer);
-    TString outputLYZ2PROD = AliAnalysisManager::GetCommonFileName();
+    AliAnalysisDataContainer *cinputLYZ2PROD = mgr->CreateContainer(Form("cobjLYZ2PRODin_%s",centralityName.Data()),
+                                                                   TList::Class(),AliAnalysisManager::kInputContainer);
+    TString outputLYZ2PROD = fileName;
     outputLYZ2PROD += ":outputLYZ2PRODanalysis";
-    outputLYZ2PROD+= type;
-    
-    AliAnalysisDataContainer *coutputLYZ2PROD = mgr->CreateContainer("cobjLYZ2PROD", TList::Class(),AliAnalysisManager::kOutputContainer,outputLYZ2PROD); 
+    outputLYZ2PROD+= rptypestr;
+    TString cobjLYZ2PROD = "cobjLYZ2PROD";
+    if((centrMin > 0)||(centrMax > 0)) {
+      cobjLYZ2PROD += "_"; cobjLYZ2PROD += centralityName.Data(); }   
+    AliAnalysisDataContainer *coutputLYZ2PROD = mgr->CreateContainer(cobjLYZ2PROD.Data(),
+                                                                    TList::Class(),AliAnalysisManager::kOutputContainer,outputLYZ2PROD); 
     mgr->ConnectInput(taskLYZ2PROD,0,coutputFE); 
     mgr->ConnectInput(taskLYZ2PROD,1,cinputLYZ2PROD);
-    mgr->ConnectOutput(taskLYZ2PROD,0,coutputLYZ2PROD); 
+    mgr->ConnectOutput(taskLYZ2PROD,1,coutputLYZ2PROD); 
     cinputLYZ2PROD->SetData(fInputListLYZ2PROD);
   }
   if(LYZEP) {
-    AliAnalysisDataContainer *cinputLYZEP = mgr->CreateContainer("cobjLYZEPin",TList::Class(),AliAnalysisManager::kInputContainer);
-    TString outputLYZEP = AliAnalysisManager::GetCommonFileName();
+    AliAnalysisDataContainer *cinputLYZEP = mgr->CreateContainer(Form("cobjLYZEPin_%s",centralityName.Data()),
+                                                                TList::Class(),AliAnalysisManager::kInputContainer);
+    TString outputLYZEP = fileName;
     outputLYZEP += ":outputLYZEPanalysis";
-    outputLYZEP+= type;
-    
-    AliAnalysisDataContainer *coutputLYZEP = mgr->CreateContainer("cobjLYZEP", TList::Class(),AliAnalysisManager::kOutputContainer,outputLYZEP); 
+    outputLYZEP+= rptypestr;
+    TString cobjLYZEP = "cobjLYZEP";
+    if((centrMin > 0)||(centrMax > 0)) {
+      cobjLYZEP += "_"; cobjLYZEP += centralityName.Data(); }       
+    AliAnalysisDataContainer *coutputLYZEP = mgr->CreateContainer(cobjLYZEP.Data(),
+                                                                 TList::Class(),AliAnalysisManager::kOutputContainer,outputLYZEP); 
     mgr->ConnectInput(taskLYZEP,0,coutputFE); 
     mgr->ConnectInput(taskLYZEP,1,cinputLYZEP);
-    mgr->ConnectOutput(taskLYZEP,0,coutputLYZEP); 
+    mgr->ConnectOutput(taskLYZEP,1,coutputLYZEP); 
     cinputLYZEP->SetData(fInputListLYZEP);
   }
   if(GFC) {
-    TString outputGFC = AliAnalysisManager::GetCommonFileName();
+    TString outputGFC = fileName;
     outputGFC += ":outputGFCanalysis";
-    outputGFC+= type;
-    
-    AliAnalysisDataContainer *coutputGFC = mgr->CreateContainer("cobjGFC", TList::Class(),AliAnalysisManager::kOutputContainer,outputGFC); 
+    outputGFC+= rptypestr;
+    TString cobjGFC = "cobjLYZ1PROD";
+    if((centrMin > 0)||(centrMax > 0)) {
+      cobjGFC += "_"; cobjGFC += centralityName.Data(); }           
+    AliAnalysisDataContainer *coutputGFC = mgr->CreateContainer(cobjGFC.Data(),
+                                                               TList::Class(),AliAnalysisManager::kOutputContainer,outputGFC); 
     mgr->ConnectInput(taskGFC,0,coutputFE); 
-    mgr->ConnectOutput(taskGFC,0,coutputGFC);
+    mgr->ConnectOutput(taskGFC,1,coutputGFC);
     if (useWeights) {
       mgr->ConnectInput(taskGFC,1,cinputWeights);
       cinputWeights->SetData(weightsList);
     } 
   }
   if(QC) {
-    TString outputQC = AliAnalysisManager::GetCommonFileName();
+    TString outputQC = fileName;
     outputQC += ":outputQCanalysis";
-    outputQC+= type;
-
-    AliAnalysisDataContainer *coutputQC = mgr->CreateContainer("cobjQC", TList::Class(),AliAnalysisManager::kOutputContainer,outputQC); 
+    outputQC+= rptypestr;
+    TString cobjQC = "cobjQC";
+    if((centrMin > 0)||(centrMax > 0)) {
+      cobjQC += "_"; cobjQC += centralityName.Data(); }       
+    AliAnalysisDataContainer *coutputQC = mgr->CreateContainer(cobjQC.Data(),
+                                                              TList::Class(),AliAnalysisManager::kOutputContainer,outputQC); 
     mgr->ConnectInput(taskQC,0,coutputFE); 
-    mgr->ConnectOutput(taskQC,0,coutputQC);
+    mgr->ConnectOutput(taskQC,1,coutputQC);
     if (useWeights) {
       mgr->ConnectInput(taskQC,1,cinputWeights);
       cinputWeights->SetData(weightsList);
-    } 
+    }
   }
   if(FQD) {
-    TString outputFQD = AliAnalysisManager::GetCommonFileName();
+    TString outputFQD = fileName;
     outputFQD += ":outputFQDanalysis";
-    outputFQD+= type;
-    
-    AliAnalysisDataContainer *coutputFQD = mgr->CreateContainer("cobjFQD", TList::Class(),AliAnalysisManager::kOutputContainer,outputFQD); 
+    outputFQD+= rptypestr;
+    TString cobjFQD = "cobjFQD";
+    if((centrMin > 0)||(centrMax > 0)) {
+      cobjFQD += "_"; cobjFQD += centralityName.Data(); }       
+    AliAnalysisDataContainer *coutputFQD = mgr->CreateContainer(cobjFQD.Data(),
+                                                               TList::Class(),AliAnalysisManager::kOutputContainer,outputFQD); 
     mgr->ConnectInput(taskFQD,0,coutputFE); 
-    mgr->ConnectOutput(taskFQD,0,coutputFQD);
+    mgr->ConnectOutput(taskFQD,1,coutputFQD);
     if(useWeights) {
       mgr->ConnectInput(taskFQD,1,cinputWeights);
       cinputWeights->SetData(weightsList);
     } 
   }
   if(MCEP) {
-    TString outputMCEP = AliAnalysisManager::GetCommonFileName();
+    TString outputMCEP = fileName;
     outputMCEP += ":outputMCEPanalysis";
-    outputMCEP+= type;
+    outputMCEP+= rptypestr;
+    TString cobjMCEP = "cobjMCEP";
+    if((centrMin > 0)||(centrMax > 0)) {
+      cobjMCEP += "_"; cobjMCEP += centralityName.Data(); }           
+    AliAnalysisDataContainer *coutputMCEP = mgr->CreateContainer(cobjMCEP.Data(),
+                                                                TList::Class(),AliAnalysisManager::kOutputContainer,outputMCEP); 
+    mgr->ConnectInput(taskMCEP,0,coutputFE);
+    mgr->ConnectOutput(taskMCEP,1,coutputMCEP); 
+  }
+  if(MH) {
+    TString outputMH = fileName;
+    outputMH += ":outputMHanalysis";
+    outputMH += rptypestr;
+    TString cobjMH = "cobjMH";
+    if((centrMin > 0)||(centrMax > 0)) {
+      cobjMH += "_"; cobjMH += centralityName.Data(); }                   
+    AliAnalysisDataContainer *coutputMH = mgr->CreateContainer(cobjMH.Data(),
+                                                              TList::Class(),AliAnalysisManager::kOutputContainer,outputMH); 
+    mgr->ConnectInput(taskMH,0,coutputFE); 
+    mgr->ConnectOutput(taskMH,1,coutputMH); 
+    //if (useWeights) {
+    //  mgr->ConnectInput(taskMH,1,cinputWeights);
+    //  cinputWeights->SetData(weightsList);
+    //} 
+  }
+  if(NL) {
+    TString outputNL = fileName;
+    outputNL += ":outputNLanalysis";
+    outputNL += rptypestr;
+    TString cobjNL = "cobjNL";
+    if((centrMin > 0)||(centrMax > 0)) {
+      cobjNL += "_"; cobjNL += centralityName.Data(); }                   
+    AliAnalysisDataContainer *coutputNL = mgr->CreateContainer(cobjNL.Data(),
+                                                              TList::Class(),AliAnalysisManager::kOutputContainer,outputNL); 
+    mgr->ConnectInput(taskNL,0,coutputFE);
+    mgr->ConnectOutput(taskNL,1,coutputNL);
+    //if (useWeights) {
+    //  mgr->ConnectInput(taskNL,1,cinputWeights);
+    //  cinputWeights->SetData(weightsList);
+    //} 
+  }
+
+  ///////////////////////////////////////////////////////////////////////////////////////////
+  if (runQAtask)
+  {
+    AliAnalysisTaskQAflow* taskQAflow = new AliAnalysisTaskQAflow("TaskQAflow");
+    taskQAflow->SetEventCuts(cutsEvent);
+    taskQAflow->SetTrackCuts(cutsRP);
+    taskQAflow->SetFillNTuple(FillQAntuple);
+    taskQAflow->SetDoCorrelations(DoQAcorrelations);
+    mgr->AddTask(taskQAflow);
     
-    AliAnalysisDataContainer *coutputMCEP = mgr->CreateContainer("cobjMCEP", TList::Class(),AliAnalysisManager::kOutputContainer,outputMCEP); 
-    mgr->ConnectInput(taskMCEP,0,coutputFE); 
-    mgr->ConnectOutput(taskMCEP,0,coutputMCEP); 
+    if((centrMin > 0)||(centrMax > 0))
+      Printf("centralityName %s",centralityName.Data());
+    TString taskQAoutputFileName(fileNameBase);
+    taskQAoutputFileName.Append("_QA.root");
+
+    TString taskQAname = "flowQA";
+    if((centrMin > 0)||(centrMax > 0)) {
+      taskQAname += "_"; taskQAname += centralityName.Data();}
+    AliAnalysisDataContainer* coutputQAtask = mgr->CreateContainer(taskQAname.Data(),
+                                              TObjArray::Class(),
+                                              AliAnalysisManager::kOutputContainer,
+                                              taskQAoutputFileName);
+
+    TString taskQATreeName = "flowQAntuple";
+    if((centrMin > 0)||(centrMax > 0)) {
+      taskQATreeName += "_"; taskQATreeName += centralityName.Data();}
+    AliAnalysisDataContainer* coutputQAtaskTree = mgr->CreateContainer(taskQATreeName.data(),
+                                              TNtuple::Class(),
+                                              AliAnalysisManager::kOutputContainer,
+                                              taskQAoutputFileName);
+    mgr->ConnectInput(taskQAflow,0,mgr->GetCommonInputContainer());
+    mgr->ConnectInput(taskQAflow,1,coutputFE);
+    mgr->ConnectOutput(taskQAflow,1,coutputQAtask);
+    if (FillQAntuple) mgr->ConnectOutput(taskQAflow,2,coutputQAtaskTree);
   }
-  
-
-  // Return analysis task
-  //===========================================================================
-  return taskFE;
-  
-
-
 }