[u/mrichter/AliRoot.git] / PWGCF / FLOW / macros / AddTaskPIDFlowSP.C


class AliAnalysisDataContainer;
class AliFlowTrackCuts;
class AliFlowTrackSimpleCuts;
class AliFlowEventCuts;
class AliFlowEventSimpleCuts;


void AddTaskPIDFlowSP(Int_t triggerSelectionString=AliVEvent::kMB,
                                   Float_t etamin=-0.8,
                                   Float_t etamax=0.8,
                                   Float_t EtaGap=0.2,
                                   TString fileNameBase="AnalysisResults",
                                   TString uniqueStr="Pion_02",
                                   TString Qvector ="Qa",
                                   Int_t AODfilterBit = 272,
                                   Int_t charge=0,
                                   Int_t MinTPCdedx = 10,
                                   Int_t ncentrality = 6,
                                   Bool_t doQA=kTRUE,
                                   Bool_t isPID = kTRUE,
                                   Bool_t VZERO = kFALSE, // use vzero sp method
                                   Bool_t is2011 = kFALSE,
                                   AliPID::EParticleType particleType=AliPID::kPion,
                                   AliFlowTrackCuts::PIDsource sourcePID=AliFlowTrackCuts::kTOFbayesian) {

// 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.;
   
//Define the range for eta subevents (for SP method) with TPC
Double_t minA = -0.8;//
Double_t maxA = -0.5*EtaGap;//
Double_t minB = +0.5*EtaGap;//
Double_t maxB = +0.8;//
    
int centrMin[9] = {0,5,10,20,30,40,50,60,70};
int centrMax[9] = {5,10,20,30,40,50,60,70,80};
const int ncentr = ncentrality;
    
    
//---------Data selection---------- ESD only!!!
//kMC, kGlobal, kESD_TPConly, kESD_SPDtracklet
AliFlowTrackCuts::trackParameterType rptype = AliFlowTrackCuts::kTPCstandalone;
AliFlowTrackCuts::trackParameterType poitype = AliFlowTrackCuts::kTPCstandalone;
    
//---------Parameter mixing--------
//kPure - no mixing, kTrackWithMCkine, kTrackWithMCPID, kTrackWithMCpt
AliFlowTrackCuts::trackParameterMix rpmix = AliFlowTrackCuts::kPure;
AliFlowTrackCuts::trackParameterMix poimix = AliFlowTrackCuts::kPure;
    
const char* rptypestr = AliFlowTrackCuts::GetParamTypeName(rptype);  //ESD
const char* poitypestr = AliFlowTrackCuts::GetParamTypeName(poitype); //ESD

    
//===========================================================================
// EVENTS CUTS:
AliFlowEventCuts* cutsEvent[ncentr];
AliFlowTrackCuts* cutsRP[ncentr];
AliFlowTrackCuts* cutsPOI[ncentr];
TString outputSlotName[ncentr][4];
TString suffixName[ncentr];
    
for(int icentr=0;icentr<ncentr;icentr++){
    cutsEvent[icentr] = new AliFlowEventCuts(Form("eventcuts_%d",icentr));
    cutsEvent[icentr]->SetLHC11h(is2011);
    cutsEvent[icentr]->SetCentralityPercentileRange(centrMin[icentr],centrMax[icentr]);
    cutsEvent[icentr]->SetCentralityPercentileMethod(AliFlowEventCuts::kV0);
    //  cutsEvent->SetRefMultMethod(AliFlowEventCuts::kVZERO);
    //cutsEvent->SetCentralityPercentileMethod(AliFlowEventCuts::kSPD1tracklets);
    //cutsEvent->SetNContributorsRange(2);
    cutsEvent[icentr]->SetPrimaryVertexZrange(-10.,10.);
    cutsEvent[icentr]->SetQA(doQA);
    cutsEvent[icentr]->SetCutTPCmultiplicityOutliers();
    
    
    // RP TRACK CUTS:
    if(!VZERO){
        cutsRP[icentr] = new AliFlowTrackCuts(Form("RP_%d",icentr));
        //cutsRP[icentr]->SetParamType(rptype);
        //cutsRP[icentr]->SetParamMix(rpmix);
        cutsRP[icentr]->SetPtRange(0.2,5.);
        cutsRP[icentr]->SetEtaRange(etamin,etamax);
        cutsRP[icentr]->SetMinNClustersTPC(70);
        cutsRP->SetMinChi2PerClusterTPC(0.1);
        cutsRP->SetMaxChi2PerClusterTPC(4.0);
        cutsRP[icentr]->SetMaxDCAToVertexXY(2.4);
        cutsRP[icentr]->SetMaxDCAToVertexZ(3.0);
        cutsRP[icentr]->SetAcceptKinkDaughters(kFALSE);
        cutsRP[icentr]->SetMinimalTPCdedx(MinTPCdedx);
        cutsRP[icentr]->SetAODfilterBit(AODfilterBit);
    }
    
    if(VZERO) { // use vzero sub analysis
        if(!is2011) cutsRP[icentr] = AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts2010(); // select vzero tracks
        if(is2011)  cutsRP[icentr] = AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts2011(); // select vzero tracks
        
        if(!cutsRP[icentr]) {
            cout << " Fatal error: no RP cuts found! " << endl;
            return 0x0;
        }

     }//vzero is not a tracking device it is just a scintillator. so pt range or DCAtoVertex are not set here.
    cutsRP[icentr]->SetQA(doQA);

    
    //POIs for SP and QC method
    //===========================================================================
    AliFlowTrackCuts  *SP_POI[ncentr];
    //half window for POIs
    //=======================SP POI Cuts
    SP_POI[icentr] = DefinePOIcuts(icentr);

    SP_POI[icentr]->GetBayesianResponse()->ForceOldDedx(); // for 2010 data to use old TPC PID Response instead of the official one
    //SP_POI[icentr]->SetParamType(poitype);
    //SP_POI[icentr]->SetParamMix(poimix);
    SP_POI[icentr]->SetPtRange(0.2,5.);//
    SP_POI[icentr]->SetMinNClustersTPC(70);
    SP_POI[icentr]->SetMinChi2PerClusterTPC(0.1);
    SP_POI[icentr]->SetMaxChi2PerClusterTPC(4.0);
    if(!VZERO && Qvector=="Qa"){
        SP_POI[icentr]->SetEtaRange( +0.5*EtaGap, etamax );
        printf(" > NOTE: Using half TPC (Qa) as POI selection u < \n");

    }
    if(!VZERO && Qvector=="Qb"){
        SP_POI[icentr]->SetEtaRange( etamin,-0.5*EtaGap );
        printf(" > NOTE: Using half TPC (Qb) as POI selection u < \n");

    }
    if(VZERO){
        SP_POI[icentr]->SetEtaRange( etamin,etamax );
        printf(" > NOTE: Using full TPC as POI selection u < \n");
    }
    //SP_POI->SetRequireITSRefit(kTRUE);
    //SP_POI->SetRequireTPCRefit(kTRUE);
    //SP_POI->SetMinNClustersITS(2);
    //SP_POI->SetMaxChi2PerClusterITS(1.e+09);
    SP_POI[icentr]->SetMaxDCAToVertexXY(2.4);
    SP_POI[icentr]->SetMaxDCAToVertexZ(3.0);
    //SP_POI->SetDCAToVertex2D(kTRUE);
    //SP_POI->SetMaxNsigmaToVertex(1.e+10);
    //SP_POI->SetRequireSigmaToVertex(kFALSE);
    SP_POI[icentr]->SetAcceptKinkDaughters(kFALSE);
    if(isPID) SP_POI[icentr]->SetPID(particleType, sourcePID);//particleType, sourcePID
    if (charge!=0) SP_POI[icentr]->SetCharge(charge);
    //SP_POI->SetAllowTOFmismatch(kFALSE);
    SP_POI[icentr]->SetRequireStrictTOFTPCagreement(kTRUE);
    SP_POI[icentr]->SetMinimalTPCdedx(MinTPCdedx);
    SP_POI[icentr]->SetAODfilterBit(AODfilterBit);
    SP_POI[icentr]->SetQA(doQA);
    SP_POI[icentr]->SetPriors((centrMin[icentr]+centrMax[icentr])*0.5);


    //=====================================================================
 
    if(!VZERO && Qvector=="Qa") suffixName[icentr] = "Qa";
    if(!VZERO && Qvector=="Qb") suffixName[icentr] = "Qb";
    if(VZERO) suffixName[icentr] = "vzero";
    suffixName[icentr] += "-highharmflow";
    suffixName[icentr] += Form("%i_", centrMin[icentr]);
    suffixName[icentr] += Form("%i_", centrMax[icentr]);
    suffixName[icentr] += Form("%.f_", EtaGap*10);

    if(isPID){
        suffixName[icentr]+=AliFlowTrackCuts::PIDsourceName(sourcePID);
        suffixName[icentr]+="_";
        suffixName[icentr]+=AliPID::ParticleName(particleType);//particleType
    }
    else{
        suffixName[icentr]+="AllCharged";
    }
    if (charge<0) suffixName[icentr]+="-";
    if (charge>0) suffixName[icentr]+="+";


    for(int harmonic=2;harmonic<6;harmonic++){  //for v2,v3,v4 and v5
        outputSlotName[icentr][harmonic-2] = "";
        outputSlotName[icentr][harmonic-2]+=uniqueStr;
        outputSlotName[icentr][harmonic-2]+=Form("_v%i_",harmonic);
        outputSlotName[icentr][harmonic-2]+=cutsRP[icentr]->GetName();
        outputSlotName[icentr][harmonic-2]+="_";
        outputSlotName[icentr][harmonic-2]+=SP_POI[icentr]->GetName();
        outputSlotName[icentr][harmonic-2]+=Form("_%i-",centrMin[icentr]);
        outputSlotName[icentr][harmonic-2]+=Form("%i_",centrMax[icentr]);
        
        
        if(isPID){
            outputSlotName[icentr][harmonic-2]+=AliFlowTrackCuts::PIDsourceName(sourcePID);//sourcePID
            outputSlotName[icentr][harmonic-2]+="_";
            outputSlotName[icentr][harmonic-2]+=AliPID::ParticleName(particleType);//particleType
        }
        else{
            outputSlotName[icentr][harmonic-2]+="AllCharged";
        }
        if (charge<0) outputSlotName[icentr][harmonic-2]+="-";
        if (charge>0) outputSlotName[icentr][harmonic-2]+="+";
    }
}


TString fileName(fileNameBase);
fileName.Append(".root");


//====================================FLOWPACKAGE TASKS=========================//
AliAnalysisDataContainer *cinput1[ncentr];
AliAnalysisDataContainer *coutputFE[ncentr];
AliAnalysisDataContainer* coutputFEQA[ncentr];
AliAnalysisTaskFlowEvent *taskFE[ncentr];

AliAnalysisDataContainer *flowEvent[ncentr][4];
AliAnalysisTaskFilterFE *tskFilter[ncentr][4];

AliAnalysisDataContainer *coutputSP[ncentr][4];
AliAnalysisTaskScalarProduct *taskSP[ncentr][4];

TString outputQA[ncentr];
TString myNameSP[ncentr][4];
TString slot[ncentr][4];

for (int icentr=0; icentr<ncentr; icentr++) {
        
    // 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 can also be checked here.
    //==============================================================================
    if (!mgr->GetInputEventHandler()) {
        ::Error("AddTaskFlowEvent", "This task requires an input event handler");
        return NULL;
    }

    taskFE[icentr] = new AliAnalysisTaskFlowEvent(Form("TaskFlowEvent_%s",suffixName[icentr].Data()),"",doQA);
    taskFE[icentr]->SelectCollisionCandidates(triggerSelectionString);
    
  //  if(taskFE[icentr]->SetVZEROSubEvents(EP3sub)) cout << " --> Setting up VZERO subevents method ... " << endl;
    
    if(!VZERO) taskFE[icentr]->SetSubeventEtaRange(minA, maxA, minB, maxB);
    if(VZERO)  taskFE[icentr]->SetSubeventEtaRange(-5,-1.5,+1.5,5);
    mgr->AddTask(taskFE[icentr]);

    // Pass cuts for RPs and POIs to the task:
    taskFE[icentr]->SetCutsEvent(cutsEvent[icentr]);
    taskFE[icentr]->SetCutsRP(cutsRP[icentr]);
    taskFE[icentr]->SetCutsPOI(SP_POI[icentr]);
    if (cutsRP[icentr]->GetParamType()==AliFlowTrackCuts::kVZERO)
    {
        //TODO: since this is set in a static object all analyses in an analysis train
        //will be affected.
        taskFE[icentr]->SetHistWeightvsPhiMin(0.);
        taskFE[icentr]->SetHistWeightvsPhiMax(200.);
    }
    
    AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
    
    cinput1[icentr] = mgr->GetCommonInputContainer();
    
    coutputFE[icentr] = mgr->CreateContainer(Form("FlowEvent_%s",suffixName[icentr].Data()),AliFlowEventSimple::Class(),AliAnalysisManager::kExchangeContainer);

    mgr->ConnectInput(taskFE[icentr],0,cinput1[icentr]);
    mgr->ConnectOutput(taskFE[icentr],1,coutputFE[icentr]);
    //==========================================================
    
    TString Species = "";
    if(isPID) Species += AliPID::ParticleName(particleType);
    else      Species += "Allcharged";
        

    for(int harm=2;harm<6;harm++){
        myNameSP[icentr][harm-2] = "SP_";
        myNameSP[icentr][harm-2] += Qvector;
        myNameSP[icentr][harm-2] += Form("_v%i_%s_%.f",harm,outputSlotName[icentr][harm-2].Data(),EtaGap*10);
        
        flowEvent[icentr][harm-2] = mgr->CreateContainer( Form("Filter_%s", myNameSP[icentr][harm-2].Data()),
                                                                   AliFlowEventSimple::Class(),
                                                                   AliAnalysisManager::kExchangeContainer );
        
        tskFilter[icentr][harm-2] = new AliAnalysisTaskFilterFE( Form("TaskFilter_%s",myNameSP[icentr][harm-2].Data()),cutsRP[icentr], NULL);
        if(!VZERO){
            tskFilter[icentr][harm-2]->SetSubeventEtaRange(etamin, -.5*EtaGap, +.5*EtaGap, etamax);
        }
        if(VZERO) tskFilter[icentr][harm-2]->SetSubeventEtaRange(-5,-1.5,+1.5,5);
        mgr->AddTask(tskFilter[icentr][harm-2]);
        mgr->ConnectInput( tskFilter[icentr][harm-2],0,coutputFE[icentr]);
        mgr->ConnectOutput(tskFilter[icentr][harm-2],1,flowEvent[icentr][harm-2]);
        
        
        taskSP[icentr][harm-2] = new AliAnalysisTaskScalarProduct(Form("TaskScalarProduct_%s",outputSlotName[icentr][harm-2].Data()),kFALSE);
        taskSP[icentr][harm-2]->SetHarmonic(harm);
        taskSP[icentr][harm-2]->SelectCollisionCandidates(triggerSelectionString);
        taskSP[icentr][harm-2]->SetRelDiffMsub(1.0);
        taskSP[icentr][harm-2]->SetTotalQvector(Qvector);
        taskSP[icentr][harm-2]->SetApplyCorrectionForNUA(kTRUE);
        
        TString outputSP = fileName;
        outputSP += ":outputSPanalysis";
        outputSP+= rptypestr;
        slot[icentr][harm-2] = "SP_";
        slot[icentr][harm-2] += outputSlotName[icentr][harm-2];
        slot[icentr][harm-2] += "_";
        slot[icentr][harm-2] += Qvector;
        coutputSP[icentr][harm-2] = mgr->CreateContainer(Form("%s_%.f",slot[icentr][harm-2].Data(),EtaGap*10),
                                         TList::Class(),AliAnalysisManager::kOutputContainer,outputSP);
        mgr->AddTask(taskSP[icentr][harm-2]);
        mgr->ConnectInput(taskSP[icentr][harm-2],0,flowEvent[icentr][harm-2]);
        mgr->ConnectInput(taskSP[icentr][harm-2],0,coutputFE[icentr]);
        mgr->ConnectOutput(taskSP[icentr][harm-2],1,coutputSP[icentr][harm-2]);
    }


    if (taskFE[icentr]->GetQAOn()) {
        outputQA[icentr] = fileName;
        outputQA[icentr] += ":QA";
        coutputFEQA[icentr] = mgr->CreateContainer(Form("QA_%s",suffixName[icentr].Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputQA[icentr]);
        mgr->ConnectOutput(taskFE[icentr],2,coutputFEQA[icentr]);
    }


}
}

//===========================================================================

AliFlowEventCuts* DefinecutsEvent(Int_t icentr){
    AliFlowEventCuts* cutsEvent = new AliFlowEventCuts(Form("eventcuts_%d",icentr));
    return cutsEvent;
}
AliFlowTrackCuts* DefineRPcuts(Int_t icentr){
    AliFlowTrackCuts* cutsRP = new AliFlowTrackCuts(Form("RP_%d",icentr));
    return cutsRP;
}
AliFlowTrackCuts* DefinePOIcuts(Int_t icentr){
    AliFlowTrackCuts* cutsPOI = new AliFlowTrackCuts(Form("POI_%d",icentr));
    return cutsPOI;
}
Commit	Line	Data
faf09db4	1
8ffcb38a	2	class AliAnalysisDataContainer;
8ffcb38a	3	class AliFlowTrackCuts;
faf09db4	4	class AliFlowTrackSimpleCuts;
8ffcb38a	5	class AliFlowEventCuts;
faf09db4	6	class AliFlowEventSimpleCuts;
8ffcb38a	7
8ffcb38a	8
db83c48d	9	void AddTaskPIDFlowSP(Int_t triggerSelectionString=AliVEvent::kMB,
8ffcb38a	10	Float_t etamin=-0.8,
	11	Float_t etamax=0.8,
	12	Float_t EtaGap=0.2,
	13	TString fileNameBase="AnalysisResults",
	14	TString uniqueStr="Pion_02",
	15	TString Qvector ="Qa",
02b20886	16	Int_t AODfilterBit = 272,
8ffcb38a	17	Int_t charge=0,
	18	Int_t MinTPCdedx = 10,
	19	Int_t ncentrality = 6,
	20	Bool_t doQA=kTRUE,
	21	Bool_t isPID = kTRUE,
db83c48d	22	Bool_t VZERO = kFALSE, // use vzero sp method
faf09db4	23	Bool_t is2011 = kFALSE,
8ffcb38a	24	AliPID::EParticleType particleType=AliPID::kPion,
	25	AliFlowTrackCuts::PIDsource sourcePID=AliFlowTrackCuts::kTOFbayesian) {
	26
	27	// Define a range of the detector to exclude
	28	Bool_t ExcludeRegion = kFALSE;
	29	Double_t excludeEtaMin = -0.;
	30	Double_t excludeEtaMax = 0.;
	31	Double_t excludePhiMin = 0.;
	32	Double_t excludePhiMax = 0.;
	33
faf09db4	34	//Define the range for eta subevents (for SP method) with TPC
8ffcb38a	35	Double_t minA = -0.8;//
	36	Double_t maxA = -0.5*EtaGap;//
	37	Double_t minB = +0.5*EtaGap;//
	38	Double_t maxB = +0.8;//
8ffcb38a	39
	40	int centrMin[9] = {0,5,10,20,30,40,50,60,70};
	41	int centrMax[9] = {5,10,20,30,40,50,60,70,80};
	42	const int ncentr = ncentrality;
	43
	44
db83c48d	45	//---------Data selection---------- ESD only!!!
8ffcb38a	46	//kMC, kGlobal, kESD_TPConly, kESD_SPDtracklet
	47	AliFlowTrackCuts::trackParameterType rptype = AliFlowTrackCuts::kTPCstandalone;
	48	AliFlowTrackCuts::trackParameterType poitype = AliFlowTrackCuts::kTPCstandalone;
	49
	50	//---------Parameter mixing--------
	51	//kPure - no mixing, kTrackWithMCkine, kTrackWithMCPID, kTrackWithMCpt
	52	AliFlowTrackCuts::trackParameterMix rpmix = AliFlowTrackCuts::kPure;
	53	AliFlowTrackCuts::trackParameterMix poimix = AliFlowTrackCuts::kPure;
	54
db83c48d	55	const char* rptypestr = AliFlowTrackCuts::GetParamTypeName(rptype); //ESD
db83c48d	56	const char* poitypestr = AliFlowTrackCuts::GetParamTypeName(poitype); //ESD
8ffcb38a	57
	58
	59	//===========================================================================
	60	// EVENTS CUTS:
	61	AliFlowEventCuts* cutsEvent[ncentr];
	62	AliFlowTrackCuts* cutsRP[ncentr];
	63	AliFlowTrackCuts* cutsPOI[ncentr];
	64	TString outputSlotName[ncentr][4];
	65	TString suffixName[ncentr];
	66
	67	for(int icentr=0;icentr<ncentr;icentr++){
	68	cutsEvent[icentr] = new AliFlowEventCuts(Form("eventcuts_%d",icentr));
02b20886	69	cutsEvent[icentr]->SetLHC11h(is2011);
8ffcb38a	70	cutsEvent[icentr]->SetCentralityPercentileRange(centrMin[icentr],centrMax[icentr]);
	71	cutsEvent[icentr]->SetCentralityPercentileMethod(AliFlowEventCuts::kV0);
	72	// cutsEvent->SetRefMultMethod(AliFlowEventCuts::kVZERO);
	73	//cutsEvent->SetCentralityPercentileMethod(AliFlowEventCuts::kSPD1tracklets);
	74	//cutsEvent->SetNContributorsRange(2);
	75	cutsEvent[icentr]->SetPrimaryVertexZrange(-10.,10.);
	76	cutsEvent[icentr]->SetQA(doQA);
	77	cutsEvent[icentr]->SetCutTPCmultiplicityOutliers();
	78
faf09db4	79
8ffcb38a	80	// RP TRACK CUTS:
db83c48d	81	if(!VZERO){
faf09db4	82	cutsRP[icentr] = new AliFlowTrackCuts(Form("RP_%d",icentr));
db83c48d	83	//cutsRP[icentr]->SetParamType(rptype);
02b20886	84	//cutsRP[icentr]->SetParamMix(rpmix);
db83c48d	85	cutsRP[icentr]->SetPtRange(0.2,5.);
	86	cutsRP[icentr]->SetEtaRange(etamin,etamax);
	87	cutsRP[icentr]->SetMinNClustersTPC(70);
02b20886	88	cutsRP->SetMinChi2PerClusterTPC(0.1);
	89	cutsRP->SetMaxChi2PerClusterTPC(4.0);
	90	cutsRP[icentr]->SetMaxDCAToVertexXY(2.4);
db83c48d	91	cutsRP[icentr]->SetMaxDCAToVertexZ(3.0);
	92	cutsRP[icentr]->SetAcceptKinkDaughters(kFALSE);
	93	cutsRP[icentr]->SetMinimalTPCdedx(MinTPCdedx);
02b20886	94	cutsRP[icentr]->SetAODfilterBit(AODfilterBit);
db83c48d	95	}
	96
	97	if(VZERO) { // use vzero sub analysis
faf09db4	98	if(!is2011) cutsRP[icentr] = AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts2010(); // select vzero tracks
	99	if(is2011) cutsRP[icentr] = AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts2011(); // select vzero tracks
	100
	101	if(!cutsRP[icentr]) {
	102	cout << " Fatal error: no RP cuts found! " << endl;
	103	return 0x0;
	104	}
	105
	106	}//vzero is not a tracking device it is just a scintillator. so pt range or DCAtoVertex are not set here.
8ffcb38a	107	cutsRP[icentr]->SetQA(doQA);
	108
	109
	110	//POIs for SP and QC method
	111	//===========================================================================
	112	AliFlowTrackCuts *SP_POI[ncentr];
	113	//half window for POIs
	114	//=======================SP POI Cuts
	115	SP_POI[icentr] = DefinePOIcuts(icentr);
	116
	117	SP_POI[icentr]->GetBayesianResponse()->ForceOldDedx(); // for 2010 data to use old TPC PID Response instead of the official one
b8256652	118	//SP_POI[icentr]->SetParamType(poitype);
b8256652	119	//SP_POI[icentr]->SetParamMix(poimix);
8ffcb38a	120	SP_POI[icentr]->SetPtRange(0.2,5.);//
8ffcb38a	121	SP_POI[icentr]->SetMinNClustersTPC(70);
02b20886	122	SP_POI[icentr]->SetMinChi2PerClusterTPC(0.1);
02b20886	123	SP_POI[icentr]->SetMaxChi2PerClusterTPC(4.0);
db83c48d	124	if(!VZERO && Qvector=="Qa"){
8ffcb38a	125	SP_POI[icentr]->SetEtaRange( +0.5*EtaGap, etamax );
db83c48d	126	printf(" > NOTE: Using half TPC (Qa) as POI selection u < \n");
db83c48d	127
8ffcb38a	128	}
db83c48d	129	if(!VZERO && Qvector=="Qb"){
8ffcb38a	130	SP_POI[icentr]->SetEtaRange( etamin,-0.5*EtaGap );
db83c48d	131	printf(" > NOTE: Using half TPC (Qb) as POI selection u < \n");
8ffcb38a	132
db83c48d	133	}
	134	if(VZERO){
	135	SP_POI[icentr]->SetEtaRange( etamin,etamax );
	136	printf(" > NOTE: Using full TPC as POI selection u < \n");
	137	}
02b20886	138	//SP_POI->SetRequireITSRefit(kTRUE);
	139	//SP_POI->SetRequireTPCRefit(kTRUE);
	140	//SP_POI->SetMinNClustersITS(2);
	141	//SP_POI->SetMaxChi2PerClusterITS(1.e+09);
	142	SP_POI[icentr]->SetMaxDCAToVertexXY(2.4);
8ffcb38a	143	SP_POI[icentr]->SetMaxDCAToVertexZ(3.0);
	144	//SP_POI->SetDCAToVertex2D(kTRUE);
	145	//SP_POI->SetMaxNsigmaToVertex(1.e+10);
	146	//SP_POI->SetRequireSigmaToVertex(kFALSE);
	147	SP_POI[icentr]->SetAcceptKinkDaughters(kFALSE);
	148	if(isPID) SP_POI[icentr]->SetPID(particleType, sourcePID);//particleType, sourcePID
	149	if (charge!=0) SP_POI[icentr]->SetCharge(charge);
	150	//SP_POI->SetAllowTOFmismatch(kFALSE);
	151	SP_POI[icentr]->SetRequireStrictTOFTPCagreement(kTRUE);
	152	SP_POI[icentr]->SetMinimalTPCdedx(MinTPCdedx);
02b20886	153	SP_POI[icentr]->SetAODfilterBit(AODfilterBit);
8ffcb38a	154	SP_POI[icentr]->SetQA(doQA);
	155	SP_POI[icentr]->SetPriors((centrMin[icentr]+centrMax[icentr])*0.5);
	156
	157
	158
	159
	160	//=====================================================================
	161
02b20886	162	if(!VZERO && Qvector=="Qa") suffixName[icentr] = "Qa";
	163	if(!VZERO && Qvector=="Qb") suffixName[icentr] = "Qb";
	164	if(VZERO) suffixName[icentr] = "vzero";
8ffcb38a	165	suffixName[icentr] += "-highharmflow";
	166	suffixName[icentr] += Form("%i_", centrMin[icentr]);
	167	suffixName[icentr] += Form("%i_", centrMax[icentr]);
	168	suffixName[icentr] += Form("%.f_", EtaGap*10);
	169
	170	if(isPID){
	171	suffixName[icentr]+=AliFlowTrackCuts::PIDsourceName(sourcePID);
	172	suffixName[icentr]+="_";
	173	suffixName[icentr]+=AliPID::ParticleName(particleType);//particleType
	174	}
	175	else{
	176	suffixName[icentr]+="AllCharged";
	177	}
	178	if (charge<0) suffixName[icentr]+="-";
	179	if (charge>0) suffixName[icentr]+="+";
	180
	181
	182	for(int harmonic=2;harmonic<6;harmonic++){ //for v2,v3,v4 and v5
	183	outputSlotName[icentr][harmonic-2] = "";
	184	outputSlotName[icentr][harmonic-2]+=uniqueStr;
	185	outputSlotName[icentr][harmonic-2]+=Form("_v%i_",harmonic);
	186	outputSlotName[icentr][harmonic-2]+=cutsRP[icentr]->GetName();
	187	outputSlotName[icentr][harmonic-2]+="_";
	188	outputSlotName[icentr][harmonic-2]+=SP_POI[icentr]->GetName();
	189	outputSlotName[icentr][harmonic-2]+=Form("_%i-",centrMin[icentr]);
	190	outputSlotName[icentr][harmonic-2]+=Form("%i_",centrMax[icentr]);
	191
	192
	193	if(isPID){
	194	outputSlotName[icentr][harmonic-2]+=AliFlowTrackCuts::PIDsourceName(sourcePID);//sourcePID
	195	outputSlotName[icentr][harmonic-2]+="_";
	196	outputSlotName[icentr][harmonic-2]+=AliPID::ParticleName(particleType);//particleType
	197	}
	198	else{
	199	outputSlotName[icentr][harmonic-2]+="AllCharged";
	200	}
	201	if (charge<0) outputSlotName[icentr][harmonic-2]+="-";
	202	if (charge>0) outputSlotName[icentr][harmonic-2]+="+";
	203	}
	204	}
	205
	206
	207	TString fileName(fileNameBase);
	208	fileName.Append(".root");
	209
	210
	211
	212	//====================================FLOWPACKAGE TASKS=========================//
	213	AliAnalysisDataContainer *cinput1[ncentr];
	214	AliAnalysisDataContainer *coutputFE[ncentr];
	215	AliAnalysisDataContainer* coutputFEQA[ncentr];
	216	AliAnalysisTaskFlowEvent *taskFE[ncentr];
	217
	218	AliAnalysisDataContainer *flowEvent[ncentr][4];
	219	AliAnalysisTaskFilterFE *tskFilter[ncentr][4];
	220
	221	AliAnalysisDataContainer *coutputSP[ncentr][4];
	222	AliAnalysisTaskScalarProduct *taskSP[ncentr][4];
	223
	224	TString outputQA[ncentr];
	225	TString myNameSP[ncentr][4];
	226	TString slot[ncentr][4];
	227
	228	for (int icentr=0; icentr<ncentr; icentr++) {
229
230	// Get the pointer to the existing analysis manager via the static access method.
231	//==============================================================================
232	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
233	if (!mgr) {
234	Error("AddTaskFlowEvent", "No analysis manager to connect to.");
235	return NULL;
236	}
237
238	// Check the analysis type using the event handlers connected to the analysis
239	// manager. The availability of MC handler can also be checked here.
240	//==============================================================================
241	if (!mgr->GetInputEventHandler()) {
242	::Error("AddTaskFlowEvent", "This task requires an input event handler");
243	return NULL;
244	}
245
246	taskFE[icentr] = new AliAnalysisTaskFlowEvent(Form("TaskFlowEvent_%s",suffixName[icentr].Data()),"",doQA);
247	taskFE[icentr]->SelectCollisionCandidates(triggerSelectionString);
faf09db4	248
	249	// if(taskFE[icentr]->SetVZEROSubEvents(EP3sub)) cout << " --> Setting up VZERO subevents method ... " << endl;
	250
	251	if(!VZERO) taskFE[icentr]->SetSubeventEtaRange(minA, maxA, minB, maxB);
	252	if(VZERO) taskFE[icentr]->SetSubeventEtaRange(-5,-1.5,+1.5,5);
8ffcb38a	253	mgr->AddTask(taskFE[icentr]);
	254
	255	// Pass cuts for RPs and POIs to the task:
	256	taskFE[icentr]->SetCutsEvent(cutsEvent[icentr]);
	257	taskFE[icentr]->SetCutsRP(cutsRP[icentr]);
	258	taskFE[icentr]->SetCutsPOI(SP_POI[icentr]);
	259	if (cutsRP[icentr]->GetParamType()==AliFlowTrackCuts::kVZERO)
	260	{
	261	//TODO: since this is set in a static object all analyses in an analysis train
	262	//will be affected.
	263	taskFE[icentr]->SetHistWeightvsPhiMin(0.);
	264	taskFE[icentr]->SetHistWeightvsPhiMax(200.);
	265	}
	266
	267	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
	268
	269	cinput1[icentr] = mgr->GetCommonInputContainer();
	270
	271	coutputFE[icentr] = mgr->CreateContainer(Form("FlowEvent_%s",suffixName[icentr].Data()),AliFlowEventSimple::Class(),AliAnalysisManager::kExchangeContainer);
	272
	273	mgr->ConnectInput(taskFE[icentr],0,cinput1[icentr]);
	274	mgr->ConnectOutput(taskFE[icentr],1,coutputFE[icentr]);
	275	//==========================================================
	276
	277	TString Species = "";
	278	if(isPID) Species += AliPID::ParticleName(particleType);
	279	else Species += "Allcharged";
	280
	281
	282	for(int harm=2;harm<6;harm++){
	283	myNameSP[icentr][harm-2] = "SP_";
	284	myNameSP[icentr][harm-2] += Qvector;
	285	myNameSP[icentr][harm-2] += Form("_v%i_%s_%.f",harm,outputSlotName[icentr][harm-2].Data(),EtaGap*10);
	286
	287	flowEvent[icentr][harm-2] = mgr->CreateContainer( Form("Filter_%s", myNameSP[icentr][harm-2].Data()),
	288	AliFlowEventSimple::Class(),
	289	AliAnalysisManager::kExchangeContainer );
	290
	291	tskFilter[icentr][harm-2] = new AliAnalysisTaskFilterFE( Form("TaskFilter_%s",myNameSP[icentr][harm-2].Data()),cutsRP[icentr], NULL);
db83c48d	292	if(!VZERO){
	293	tskFilter[icentr][harm-2]->SetSubeventEtaRange(etamin, -.5EtaGap, +.5EtaGap, etamax);
	294	}
faf09db4	295	if(VZERO) tskFilter[icentr][harm-2]->SetSubeventEtaRange(-5,-1.5,+1.5,5);
8ffcb38a	296	mgr->AddTask(tskFilter[icentr][harm-2]);
	297	mgr->ConnectInput( tskFilter[icentr][harm-2],0,coutputFE[icentr]);
	298	mgr->ConnectOutput(tskFilter[icentr][harm-2],1,flowEvent[icentr][harm-2]);
	299
	300
	301	taskSP[icentr][harm-2] = new AliAnalysisTaskScalarProduct(Form("TaskScalarProduct_%s",outputSlotName[icentr][harm-2].Data()),kFALSE);
	302	taskSP[icentr][harm-2]->SetHarmonic(harm);
	303	taskSP[icentr][harm-2]->SelectCollisionCandidates(triggerSelectionString);
	304	taskSP[icentr][harm-2]->SetRelDiffMsub(1.0);
	305	taskSP[icentr][harm-2]->SetTotalQvector(Qvector);
	306	taskSP[icentr][harm-2]->SetApplyCorrectionForNUA(kTRUE);
	307
	308	TString outputSP = fileName;
	309	outputSP += ":outputSPanalysis";
	310	outputSP+= rptypestr;
	311	slot[icentr][harm-2] = "SP_";
	312	slot[icentr][harm-2] += outputSlotName[icentr][harm-2];
	313	slot[icentr][harm-2] += "_";
	314	slot[icentr][harm-2] += Qvector;
	315	coutputSP[icentr][harm-2] = mgr->CreateContainer(Form("%s_%.f",slot[icentr][harm-2].Data(),EtaGap*10),
	316	TList::Class(),AliAnalysisManager::kOutputContainer,outputSP);
	317	mgr->AddTask(taskSP[icentr][harm-2]);
	318	mgr->ConnectInput(taskSP[icentr][harm-2],0,flowEvent[icentr][harm-2]);
faf09db4	319	mgr->ConnectInput(taskSP[icentr][harm-2],0,coutputFE[icentr]);
8ffcb38a	320	mgr->ConnectOutput(taskSP[icentr][harm-2],1,coutputSP[icentr][harm-2]);
	321	}
	322
	323
	324	if (taskFE[icentr]->GetQAOn()) {
	325	outputQA[icentr] = fileName;
	326	outputQA[icentr] += ":QA";
	327	coutputFEQA[icentr] = mgr->CreateContainer(Form("QA_%s",suffixName[icentr].Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputQA[icentr]);
	328	mgr->ConnectOutput(taskFE[icentr],2,coutputFEQA[icentr]);
	329	}
	330
	331
	332	}
	333	}
	334
	335	//===========================================================================
	336
	337	AliFlowEventCuts* DefinecutsEvent(Int_t icentr){
	338	AliFlowEventCuts* cutsEvent = new AliFlowEventCuts(Form("eventcuts_%d",icentr));
	339	return cutsEvent;
	340	}
	341	AliFlowTrackCuts* DefineRPcuts(Int_t icentr){
db83c48d	342	AliFlowTrackCuts* cutsRP = new AliFlowTrackCuts(Form("RP_%d",icentr));
8ffcb38a	343	return cutsRP;
	344	}
	345	AliFlowTrackCuts* DefinePOIcuts(Int_t icentr){
db83c48d	346	AliFlowTrackCuts* cutsPOI = new AliFlowTrackCuts(Form("POI_%d",icentr));
8ffcb38a	347	return cutsPOI;
	348	}
	349
	350
	351
	352