[u/mrichter/AliRoot.git] / PWGHF / hfe / macros / configs / pPb / ConfigHFEnpepPb.C

Bool_t ReadContaminationFunctionsBeauty(TString filename, TF1 **functions, double sigma){
  TFile *in = TFile::Open(Form("$ALICE_ROOT/PWGHF/hfe/macros/configs/pPb/%s", filename.Data()));
  gROOT->cd();
  int isig = static_cast<int>(sigma * 100.);
  if (isig == -44) isig = -42;
  if (isig == 6) isig = 9;
  printf("Getting hadron background for the sigma cut: %d\n", isig);
  bool status = kTRUE;
  for(int icent = 0; icent < 12; icent++){
    functions[icent] = dynamic_cast<TF1 *>(in->Get(Form("hback_%d_%d", isig, icent)));
    if(functions[icent]) printf("Config for centrality class %d found\n", icent);
    else{
      printf("Config for the centrality class %d not found\n", icent);
      status = kFALSE;
    }
  }
  delete in;
  return status;
}

Bool_t ReadContaminationFunctions(TString filename, TF1 **functions, double sigma){
  TFile *in = TFile::Open(Form("$ALICE_ROOT/PWGHF/hfe/macros/configs/pPb/%s", filename.Data()));
  gROOT->cd();
  int isig = static_cast<int>(sigma * 100.);
  printf("Getting hadron background for the sigma cut: %d\n", isig);
  bool status = kTRUE;
  for(int icent = 0; icent < 12; icent++){
    functions[icent] = dynamic_cast<TF1 *>(in->Get(Form("hback_%d_%d", isig, icent)));
    if(functions[icent]) printf("Config for centrality class %d found\n", icent);
    else{
      printf("Config for the centrality class %d not found\n", icent);
      status = kFALSE;
    }
  }
  delete in;
  return status;
}


AliAnalysisTaskHFE* ConfigHFEnpepPb(Bool_t useMC, Bool_t isAOD, Bool_t isBeauty, TString appendix,
                UChar_t TPCcl=70, UChar_t TPCclPID = 80, 
                UChar_t ITScl=3, Double_t DCAxy=1000., Double_t DCAz=1000., 
                Double_t* tpcdEdxcutlow=NULL, Double_t* tpcdEdxcuthigh=NULL, 
                Double_t TOFs=3., Int_t TOFmis=0, 
                Int_t itshitpixel = 0, Int_t icent, 
		Double_t etami=-0.8, Double_t etama=0.8,
		Double_t phimi=-1., Double_t phima=-1.,
		Double_t assETAm=-0.8, Double_t assETAp=0.8,
		Double_t assMinPt=0.2, Int_t assITS=2, 
                Int_t assTPCcl=100, Int_t assTPCPIDcl=80, 
                Double_t assDCAr=1.0, Double_t assDCAz=2.0, 
                Double_t *assTPCSminus=NULL, Double_t *assTPCSplus=NULL, 
                Double_t assITSpid=-3., 
		Double_t assTOFs=3.,
                Bool_t useCat1Tracks = kTRUE, Bool_t useCat2Tracks = kTRUE, Int_t weightlevelback = -1, 
	        Bool_t nonPhotonicElectronBeauty = kFALSE, Bool_t ipCharge = kFALSE, Bool_t ipOpp = kFALSE, Int_t ipsys = 0)
{
  Bool_t kAnalyseTaggedTracks = kTRUE;
  Bool_t kApplyPreselection = kFALSE;

  //***************************************//
  //        Setting up the HFE cuts        //
  //***************************************//

  AliHFEcuts *hfecuts = new AliHFEcuts(appendix,"HFE cuts for pPb");
  //hfecuts->SetQAOn();
  hfecuts->CreateStandardCuts();
  hfecuts->SetMinNClustersTPC(TPCcl);
  hfecuts->SetMinNClustersTPCPID(TPCclPID);
  hfecuts->SetMinNClustersITS(ITScl);
  hfecuts->SetMinRatioTPCclusters(0.6);
  hfecuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
  hfecuts->SetCutITSpixel(itshitpixel);
  hfecuts->SetCheckITSLayerStatus(kFALSE);
  hfecuts->SetEtaRange(etami,etama);
  if(phimi >= 0. && phima >= 0) hfecuts->SetPhiRange(phimi,phima);
  hfecuts->SetRejectKinkDaughters();
  hfecuts->SetAcceptKinkMothers();
  if(isAOD) hfecuts->SetAODFilterBit(4);
  
  //if((iPixelAny==AliHFEextraCuts::kAny) || (iPixelAny==AliHFEextraCuts::kSecond))     
 
  hfecuts->SetMaxImpactParam(DCAxy,DCAz);
  hfecuts->SetUseMixedVertex(kTRUE);
  hfecuts->SetVertexRange(10.);
  // New pPb cuts (February 2013)
  hfecuts->SetUseCorrelationVertex();
  hfecuts->SetSPDVtxResolutionCut();
  hfecuts->SetpApileupCut();

  Bool_t ipSig = kFALSE;

  hfecuts->SetIPcutParam(0.0054,0.078,-0.56,0,ipSig,ipCharge,ipOpp);
  if(ipsys==1) hfecuts->SetIPcutParam(0.0054,0.057,-0.66,0,ipSig,ipCharge,ipOpp);
  if(ipsys==2) hfecuts->SetIPcutParam(0.012,0.088,-0.65,0,ipSig,ipCharge,ipOpp);

  if(isBeauty) hfecuts->SetProductionVertex(0,100,0,100);

  // TOF settings:
  Int_t usetof=0;
  Bool_t kTOFmis=kFALSE;
  if (TOFs>0.){
    usetof = 1;
    printf("CONFIGURATION FILE: TOF is used \n");
    hfecuts->SetTOFPIDStep(kTRUE);
    printf("CONFIGURATION FILE: TOF PID step is requested !!!! \n");
    if (TOFmis>0){
      kTOFmis = kTRUE;
      printf("CONFIGURATION FILE: TOF mismatch rejection is set ON \n");
    }
  }

  //***************************************//
  //        Setting up the task            //
  //***************************************//

  AliAnalysisTaskHFE *task = new AliAnalysisTaskHFE(Form("HFEtask%s",appendix.Data()));
  printf("task %p\n", task);
  task->SetpPbAnalysis();
  if(!isAOD) task->SetRemoveFirstEventInChunk();
  task->SetRemovePileUp(kFALSE);
  task->SetHFECuts(hfecuts);
  task->GetPIDQAManager()->SetHighResolutionHistos();
  task->SetRejectKinkMother(kFALSE);

  // Determine the centrality estimator
  task->SetCentralityEstimator("V0A");
  if (icent == 2) task->SetCentralityEstimator("V0M");
  else if (icent == 3) task->SetCentralityEstimator("CL1");
  else if (icent == 4) task->SetCentralityEstimator("ZNA");

  //***************************************//
  //        Prepare preselection           //
  // This mimics the ESD->AOD filter in    //
  // case of the ESD analysis and selects  //
  // only tracks which will be selected in //
  // the AOD analysis with the given filter//
  // bit. Not to be applied for AODS.      //
  // For pPb the cuts used are (bit 4)     //
  // esdTrackCutsHG0 from file $ALICE_ROOT///
  // ANALYSIS/macros/AddTaskESDFilter.C    //
  //***************************************//

  if(kApplyPreselection){    
    AliESDtrackCuts* esdfilter = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
    esdfilter->SetMaxDCAToVertexXY(2.4);
    esdfilter->SetMaxDCAToVertexZ(3.2);
    esdfilter->SetDCAToVertex2D(kTRUE);

    task->SetHFECutsPreselect(esdfilter);
    printf("Put a preselection cut\n");
    task->SetFillNoCuts(kTRUE);
  }

  //***************************************//
  //          Variable manager             //
  //***************************************//
  // Define Variables
  Double_t ptbinning[36] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
  Double_t etabinning[17] = {-0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8};

  Int_t sizept=(sizeof(ptbinning)/sizeof(double))-1;
  Int_t sizeeta=(sizeof(etabinning)/sizeof(double))-1;

  AliHFEvarManager *vm = task->GetVarManager();
  vm->AddVariable("pt", sizept, ptbinning);
  vm->AddVariable("eta", sizeeta, -0.8,0.8);
  vm->AddVariable("phi",18, -0, 2*TMath::Pi());
  vm->AddVariable("charge");
  vm->AddVariable("source");
  vm->AddVariable("centrality");

  // For the moment, remove the part dedicated to the background subtraction.
  // It will be implemented in a different way, reading it from a root file.

  //***************************************//
  //          Configure the PID            //
  //***************************************//

  // Define PID
  AliHFEpid *pid = task->GetPID();
  if(useMC) pid->SetHasMCData(kTRUE);

  if (usetof){
    pid->AddDetector("TOF", 0);
    pid->AddDetector("TPC", 1);
  } else {
    pid->AddDetector("TPC", 0);
  }
  
  // Configure TPC PID
  // do the identical thing in data and MC
  Double_t paramsTPCdEdxcutlow[12] ={0.0, 0.0, 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
  if(tpcdEdxcutlow) memcpy(paramsTPCdEdxcutlow,tpcdEdxcutlow,sizeof(paramsTPCdEdxcutlow));

  Double_t paramsTPCdEdxcuthigh[12] ={3.0, 3.0, 3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0};
  if(tpcdEdxcuthigh) memcpy(paramsTPCdEdxcuthigh,tpcdEdxcuthigh,sizeof(paramsTPCdEdxcuthigh));

  char *cutmodel;

  if(useMC){ // constant (default) cut for MC
      cutmodel="pol0(0)";
      Double_t params[1];
      params[0]=paramsTPCdEdxcutlow[0];
      pid->ConfigureTPCdefaultCut(cutmodel, params,tpcdEdxcuthigh[0]);
  } else { // correct for mean shift in data
      cutmodel="min(pol1(0),pol0(2))";
      Double_t params[3];
      //params[0]=-0.12; params[1]=0.14; params[2]=0.09;
      params[0]=-0.21 + paramsTPCdEdxcutlow[0];
      params[1]=0.14;
      params[2]=paramsTPCdEdxcutlow[0];
      pid->ConfigureTPCdefaultCut(cutmodel, params,tpcdEdxcuthigh[0]);
  }
  /*
  char *cutmodel;
  cutmodel="pol0";

  for(Int_t a=0;a<11;a++){
    // Not necessary anymore, since the pPb case is handled similarly to the pp case
    //   cout << a << " " << paramsTPCdEdxcut[a] << endl;
    Double_t tpcparamlow[1]={paramsTPCdEdxcutlow[a]};
    Float_t tpcparamhigh=paramsTPCdEdxcuthigh[a];
    pid->ConfigureTPCcentralityCut(a,cutmodel,tpcparamlow,tpcparamhigh);
  }
  pid->ConfigureTPCdefaultCut(cutmodel,paramsTPCdEdxcutlow,paramsTPCdEdxcuthigh[0]); // After introducing the pPb flag, pPb is merged with pp and this line defines the cut
  */

  // Configure TOF PID
  if (usetof){
    pid->ConfigureTOF(TOFs);
    AliHFEpidTOF *tofpid = pid->GetDetPID(AliHFEpid::kTOFpid);
    if (kTOFmis){
      tofpid->SetRejectTOFmismatch();
    }
  }

  // To make different upper TOF cut to see contamination effect
  // The below two lines should be removed after this check
  //AliHFEpidTOF *tofpid = pid->GetDetPID(AliHFEpid::kTOFpid);
  //if(TOFs<3.) tofpid->SetTOFnSigmaBand(-3,TOFs); //only to check the assymmetric tof cut

  // Load hadron background
  if(!useMC){
    Bool_t status = kTRUE;
    TF1 *hBackground[12];
    if(isAOD==1) {
	if (usetof)  status = ReadContaminationFunctions("hadroncontamination_AOD139_TOFPID_pPb_eta06.root", hBackground, tpcdEdxcutlow[0]);
	else status = ReadContaminationFunctions("hadroncontamination_AOD139_noTOFPID_pPb_eta06.root", hBackground, tpcdEdxcutlow[0]);
    }
    else if (isBeauty==1) {
	status = ReadContaminationFunctionsBeauty("hadroncontamination_ESD_Beauty_TOFPID_pPb_eta06.root", hBackground, tpcdEdxcutlow[0]);
    }
  else  status = ReadContaminationFunctions("hadroncontamination_TOFTPC_pPb_eta06_newsplines_try3.root", hBackground, tpcdEdxcutlow[0]);
    for(Int_t a=0;a<12;a++) {
      //printf("back %f \n",hBackground[a]);
      if(status) task->SetBackGroundFactorsFunction(hBackground[a],a);
      else printf("not all background functions found\n");
    }
  }

  //***************************************//
  //       Configure NPE plugin            //
  //***************************************//

  AliHFENonPhotonicElectron *backe = new AliHFENonPhotonicElectron(Form("HFEBackGroundSubtractionPID2%s",appendix.Data()),"Background subtraction");  //appendix
    //Setting the Cuts for the Associated electron-pool
  AliHFEcuts *hfeBackgroundCuts = new AliHFEcuts(Form("HFEBackSub%s",appendix.Data()),"Background sub Cuts");
  //  hfeBackgroundCuts->SetEtaRange(assETA);
  hfeBackgroundCuts->SetEtaRange(assETAm,assETAp);
  hfeBackgroundCuts->SetPtRange(assMinPt,20.);

  hfeBackgroundCuts->SetMaxChi2perClusterTPC(4);
  hfeBackgroundCuts->SetMinNClustersITS(assITS);
  hfeBackgroundCuts->SetMinNClustersTPC(assTPCcl);
  hfeBackgroundCuts->SetMinNClustersTPCPID(assTPCPIDcl);
  hfeBackgroundCuts->SetMaxImpactParam(assDCAr,assDCAz);
  if(isAOD) hfeBackgroundCuts->SetAODFilterBit(4);
  hfeBackgroundCuts->SetQAOn();			        // QA

  AliHFEpid *pidbackground = backe->GetPIDBackground();
  if(useMC) pidbackground->SetHasMCData(kTRUE);

  if (assTOFs>0.){
    pidbackground->AddDetector("TOF", 0);
    pidbackground->AddDetector("TPC", 1);
  } else {
    pidbackground->AddDetector("TPC", 0);
  }

  Double_t paramsTPCdEdxcutlowAssoc[12] ={-3.0,-3.0,-3.0,-3.0,-3.0,-3.0,-3.0,-3.0,-3.0,-3.0,-3.0,-3.0};
  if(assTPCSminus) memcpy(paramsTPCdEdxcutlowAssoc,assTPCSminus,sizeof(paramsTPCdEdxcutlowAssoc));

  Double_t paramsTPCdEdxcuthighAssoc[12] ={3.0, 3.0, 3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0};
  if(assTPCSplus) memcpy(paramsTPCdEdxcuthighAssoc,assTPCSplus,sizeof(paramsTPCdEdxcuthighAssoc));
    
  char *cutmodelAssoc;
  cutmodelAssoc="pol0";
  for(Int_t a=0;a<11;a++){
    // Not necessary anymore, since the pPb case is handled similarly to the pp case
    //   cout << a << " " << paramsTPCdEdxcut[a] << endl;
    Double_t tpcparamlow[1]={paramsTPCdEdxcutlowAssoc[a]};
    Float_t tpcparamhigh=paramsTPCdEdxcuthighAssoc[a];
    pidbackground->ConfigureTPCcentralityCut(a,cutmodelAssoc,tpcparamlow,tpcparamhigh);
  }
  pidbackground->ConfigureTPCdefaultCut(cutmodelAssoc,paramsTPCdEdxcutlowAssoc,paramsTPCdEdxcuthighAssoc[0]); // After introducing the pPb flag, pPb is merged with pp and this line defines the cut
  //backe->GetPIDBackgroundQAManager()->SetHighResolutionHistos();

  if (assTOFs>0.){
    pidbackground->ConfigureTOF(TOFs);
  }

  backe->SetHFEBackgroundCuts(hfeBackgroundCuts);

  // Selection of associated tracks for the pool
  if(useCat1Tracks) backe->SelectCategory1Tracks(kTRUE);
  if(useCat2Tracks){
    backe->SelectCategory2Tracks(kTRUE);
    backe->SetITSMeanShift(-0.5);
    backe->SetITSnSigmaHigh(assITSpid);
    Double_t assITSminus = -1.0 * assITSpid;
    backe->SetITSnSigmaLow(assITSminus);
    //backe->SetminPt(assMinPt);
  }

  // apply opening angle cut to reduce file size
  backe->SetMaxInvMass(0.3);
  backe->SetPtBinning(sizept, ptbinning);
  backe->SetEtaBinning(sizeeta, etabinning);
  // MC weight
  if(useMC) {
    //printf("test put weight %d\n",weightlevelback);
    if((weightlevelback >=0) && (weightlevelback < 3)) backe->SetWithWeights(weightlevelback);
  }
  task->SetHFEBackgroundSubtraction(backe);

  //task->SetWeightHist(); 

  //***************************************//
  //          V0 tagged tracks             //
  //***************************************//

  if(kAnalyseTaggedTracks){
    AliHFEcuts *v0trackCuts = new AliHFEcuts("V0trackCuts", "Track Cuts for tagged track Analysis");
    v0trackCuts->CreateStandardCuts();
    v0trackCuts->SetMinNClustersTPC(TPCcl);
    v0trackCuts->SetMinNClustersTPCPID(TPCclPID);
    v0trackCuts->SetMinRatioTPCclusters(0.6);
    v0trackCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
    v0trackCuts->SetMinNClustersITS(1);
    v0trackCuts->SetCutITSpixel(AliHFEextraCuts::kFirst);
    v0trackCuts->SetCheckITSLayerStatus(kFALSE);
    v0trackCuts->UnsetVertexRequirement();
    //hfecuts->SetSigmaToVertex(10);
    if(usetof) v0trackCuts->SetTOFPIDStep(kTRUE);
    v0trackCuts->SetQAOn();

    task->SwitchOnPlugin(AliAnalysisTaskHFE::kTaggedTrackAnalysis);
    task->SetTaggedTrackCuts(v0trackCuts);
    task->SetCleanTaggedTrack(kTRUE);
  }

  // QA
  printf("task %p\n", task);
  task->SetQAOn(AliAnalysisTaskHFE::kPIDqa);
  task->SetQAOn(AliAnalysisTaskHFE::kMCqa);
  task->SwitchOnPlugin(AliAnalysisTaskHFE::kDEstep);
  if(nonPhotonicElectronBeauty) task->SwitchOnPlugin(AliAnalysisTaskHFE::kNonPhotonicElectronBeauty);
  else task->SwitchOnPlugin(AliAnalysisTaskHFE::kNonPhotonicElectron);

  printf("*************************************\n");
  printf("Configuring standard Task:\n");
  task->PrintStatus();
  pid->PrintStatus();
  printf("*************************************\n");
  return task;
}
Commit	Line	Data
39ebaf78	1	Bool_t ReadContaminationFunctionsBeauty(TString filename, TF1 **functions, double sigma){
	2	TFile *in = TFile::Open(Form("$ALICE_ROOT/PWGHF/hfe/macros/configs/pPb/%s", filename.Data()));
	3	gROOT->cd();
	4	int isig = static_cast<int>(sigma * 100.);
	5	if (isig == -44) isig = -42;
	6	if (isig == 6) isig = 9;
	7	printf("Getting hadron background for the sigma cut: %d\n", isig);
	8	bool status = kTRUE;
	9	for(int icent = 0; icent < 12; icent++){
	10	functions[icent] = dynamic_cast<TF1 *>(in->Get(Form("hback_%d_%d", isig, icent)));
	11	if(functions[icent]) printf("Config for centrality class %d found\n", icent);
	12	else{
	13	printf("Config for the centrality class %d not found\n", icent);
	14	status = kFALSE;
	15	}
	16	}
	17	delete in;
	18	return status;
	19	}
	20
8dc7fa14	21	Bool_t ReadContaminationFunctions(TString filename, TF1 **functions, double sigma){
420bd2ea	22	TFile *in = TFile::Open(Form("$ALICE_ROOT/PWGHF/hfe/macros/configs/pPb/%s", filename.Data()));
8dc7fa14	23	gROOT->cd();
	24	int isig = static_cast<int>(sigma * 100.);
	25	printf("Getting hadron background for the sigma cut: %d\n", isig);
	26	bool status = kTRUE;
	27	for(int icent = 0; icent < 12; icent++){
	28	functions[icent] = dynamic_cast<TF1 *>(in->Get(Form("hback_%d_%d", isig, icent)));
	29	if(functions[icent]) printf("Config for centrality class %d found\n", icent);
	30	else{
	31	printf("Config for the centrality class %d not found\n", icent);
	32	status = kFALSE;
	33	}
	34	}
	35	delete in;
	36	return status;
	37	}
	38
39ebaf78	39
39ebaf78	40	AliAnalysisTaskHFE* ConfigHFEnpepPb(Bool_t useMC, Bool_t isAOD, Bool_t isBeauty, TString appendix,
8dc7fa14	41	UChar_t TPCcl=70, UChar_t TPCclPID = 80,
	42	UChar_t ITScl=3, Double_t DCAxy=1000., Double_t DCAz=1000.,
	43	Double_t* tpcdEdxcutlow=NULL, Double_t* tpcdEdxcuthigh=NULL,
	44	Double_t TOFs=3., Int_t TOFmis=0,
	45	Int_t itshitpixel = 0, Int_t icent,
b74766cf	46	Double_t etami=-0.8, Double_t etama=0.8,
39ebaf78	47	Double_t phimi=-1., Double_t phima=-1.,
	48	Double_t assETAm=-0.8, Double_t assETAp=0.8,
	49	Double_t assMinPt=0.2, Int_t assITS=2,
8dc7fa14	50	Int_t assTPCcl=100, Int_t assTPCPIDcl=80,
	51	Double_t assDCAr=1.0, Double_t assDCAz=2.0,
	52	Double_t assTPCSminus=NULL, Double_t assTPCSplus=NULL,
39ebaf78	53	Double_t assITSpid=-3.,
	54	Double_t assTOFs=3.,
	55	Bool_t useCat1Tracks = kTRUE, Bool_t useCat2Tracks = kTRUE, Int_t weightlevelback = -1,
	56	Bool_t nonPhotonicElectronBeauty = kFALSE, Bool_t ipCharge = kFALSE, Bool_t ipOpp = kFALSE, Int_t ipsys = 0)
8dc7fa14	57	{
39ebaf78	58	Bool_t kAnalyseTaggedTracks = kTRUE;
b74766cf	59	Bool_t kApplyPreselection = kFALSE;
420bd2ea	60
8dc7fa14	61	//***************************************//
	62	// Setting up the HFE cuts //
	63	//***************************************//
420bd2ea	64
420bd2ea	65	AliHFEcuts *hfecuts = new AliHFEcuts(appendix,"HFE cuts for pPb");
8dc7fa14	66	//hfecuts->SetQAOn();
	67	hfecuts->CreateStandardCuts();
	68	hfecuts->SetMinNClustersTPC(TPCcl);
	69	hfecuts->SetMinNClustersTPCPID(TPCclPID);
	70	hfecuts->SetMinNClustersITS(ITScl);
	71	hfecuts->SetMinRatioTPCclusters(0.6);
	72	hfecuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
	73	hfecuts->SetCutITSpixel(itshitpixel);
	74	hfecuts->SetCheckITSLayerStatus(kFALSE);
	75	hfecuts->SetEtaRange(etami,etama);
39ebaf78	76	if(phimi >= 0. && phima >= 0) hfecuts->SetPhiRange(phimi,phima);
1452614c	77	hfecuts->SetRejectKinkDaughters();
1452614c	78	hfecuts->SetAcceptKinkMothers();
8dc7fa14	79	if(isAOD) hfecuts->SetAODFilterBit(4);
	80
	81	//if((iPixelAny==AliHFEextraCuts::kAny) \|\| (iPixelAny==AliHFEextraCuts::kSecond))
8dc7fa14	82
	83	hfecuts->SetMaxImpactParam(DCAxy,DCAz);
	84	hfecuts->SetUseMixedVertex(kTRUE);
	85	hfecuts->SetVertexRange(10.);
	86	// New pPb cuts (February 2013)
	87	hfecuts->SetUseCorrelationVertex();
	88	hfecuts->SetSPDVtxResolutionCut();
	89	hfecuts->SetpApileupCut();
	90
39ebaf78	91	Bool_t ipSig = kFALSE;
	92
	93	hfecuts->SetIPcutParam(0.0054,0.078,-0.56,0,ipSig,ipCharge,ipOpp);
	94	if(ipsys==1) hfecuts->SetIPcutParam(0.0054,0.057,-0.66,0,ipSig,ipCharge,ipOpp);
	95	if(ipsys==2) hfecuts->SetIPcutParam(0.012,0.088,-0.65,0,ipSig,ipCharge,ipOpp);
	96
	97	if(isBeauty) hfecuts->SetProductionVertex(0,100,0,100);
	98
8dc7fa14	99	// TOF settings:
	100	Int_t usetof=0;
	101	Bool_t kTOFmis=kFALSE;
	102	if (TOFs>0.){
	103	usetof = 1;
	104	printf("CONFIGURATION FILE: TOF is used \n");
	105	hfecuts->SetTOFPIDStep(kTRUE);
	106	printf("CONFIGURATION FILE: TOF PID step is requested !!!! \n");
	107	if (TOFmis>0){
	108	kTOFmis = kTRUE;
	109	printf("CONFIGURATION FILE: TOF mismatch rejection is set ON \n");
	110	}
	111	}
	112
	113	//***************************************//
	114	// Setting up the task //
	115	//***************************************//
	116
	117	AliAnalysisTaskHFE *task = new AliAnalysisTaskHFE(Form("HFEtask%s",appendix.Data()));
	118	printf("task %p\n", task);
	119	task->SetpPbAnalysis();
8dc7fa14	120	if(!isAOD) task->SetRemoveFirstEventInChunk();
420bd2ea	121	task->SetRemovePileUp(kFALSE);
420bd2ea	122	task->SetHFECuts(hfecuts);
8dc7fa14	123	task->GetPIDQAManager()->SetHighResolutionHistos();
1452614c	124	task->SetRejectKinkMother(kFALSE);
8dc7fa14	125
420bd2ea	126	// Determine the centrality estimator
	127	task->SetCentralityEstimator("V0A");
	128	if (icent == 2) task->SetCentralityEstimator("V0M");
	129	else if (icent == 3) task->SetCentralityEstimator("CL1");
	130	else if (icent == 4) task->SetCentralityEstimator("ZNA");
8dc7fa14	131
1452614c	132	//***************************************//
	133	// Prepare preselection //
	134	// This mimics the ESD->AOD filter in //
	135	// case of the ESD analysis and selects //
	136	// only tracks which will be selected in //
	137	// the AOD analysis with the given filter//
	138	// bit. Not to be applied for AODS. //
	139	// For pPb the cuts used are (bit 4) //
	140	// esdTrackCutsHG0 from file $ALICE_ROOT///
	141	// ANALYSIS/macros/AddTaskESDFilter.C //
	142	//***************************************//
	143
	144	if(kApplyPreselection){
	145	AliESDtrackCuts* esdfilter = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
	146	esdfilter->SetMaxDCAToVertexXY(2.4);
	147	esdfilter->SetMaxDCAToVertexZ(3.2);
	148	esdfilter->SetDCAToVertex2D(kTRUE);
	149
	150	task->SetHFECutsPreselect(esdfilter);
	151	printf("Put a preselection cut\n");
	152	task->SetFillNoCuts(kTRUE);
	153	}
	154
8dc7fa14	155	//***************************************//
	156	// Variable manager //
	157	//***************************************//
8dc7fa14	158	// Define Variables
	159	Double_t ptbinning[36] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
	160	Double_t etabinning[17] = {-0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8};
	161
	162	Int_t sizept=(sizeof(ptbinning)/sizeof(double))-1;
	163	Int_t sizeeta=(sizeof(etabinning)/sizeof(double))-1;
	164
	165	AliHFEvarManager *vm = task->GetVarManager();
	166	vm->AddVariable("pt", sizept, ptbinning);
	167	vm->AddVariable("eta", sizeeta, -0.8,0.8);
39ebaf78	168	vm->AddVariable("phi",18, -0, 2*TMath::Pi());
8dc7fa14	169	vm->AddVariable("charge");
	170	vm->AddVariable("source");
	171	vm->AddVariable("centrality");
	172
8dc7fa14	173	// For the moment, remove the part dedicated to the background subtraction.
8dc7fa14	174	// It will be implemented in a different way, reading it from a root file.
8dc7fa14	175
	176	//***************************************//
	177	// Configure the PID //
	178	//***************************************//
	179
	180	// Define PID
	181	AliHFEpid *pid = task->GetPID();
	182	if(useMC) pid->SetHasMCData(kTRUE);
	183
	184	if (usetof){
	185	pid->AddDetector("TOF", 0);
	186	pid->AddDetector("TPC", 1);
	187	} else {
	188	pid->AddDetector("TPC", 0);
	189	}
	190
	191	// Configure TPC PID
420bd2ea	192	// do the identical thing in data and MC
	193	Double_t paramsTPCdEdxcutlow[12] ={0.0, 0.0, 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
	194	if(tpcdEdxcutlow) memcpy(paramsTPCdEdxcutlow,tpcdEdxcutlow,sizeof(paramsTPCdEdxcutlow));
8dc7fa14	195
420bd2ea	196	Double_t paramsTPCdEdxcuthigh[12] ={3.0, 3.0, 3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0};
	197	if(tpcdEdxcuthigh) memcpy(paramsTPCdEdxcuthigh,tpcdEdxcuthigh,sizeof(paramsTPCdEdxcuthigh));
	198
39ebaf78	199	char *cutmodel;
eba5ec88	200
39ebaf78	201	if(useMC){ // constant (default) cut for MC
	202	cutmodel="pol0(0)";
	203	Double_t params[1];
	204	params[0]=paramsTPCdEdxcutlow[0];
	205	pid->ConfigureTPCdefaultCut(cutmodel, params,tpcdEdxcuthigh[0]);
	206	} else { // correct for mean shift in data
	207	cutmodel="min(pol1(0),pol0(2))";
	208	Double_t params[3];
	209	//params[0]=-0.12; params[1]=0.14; params[2]=0.09;
	210	params[0]=-0.21 + paramsTPCdEdxcutlow[0];
	211	params[1]=0.14;
	212	params[2]=paramsTPCdEdxcutlow[0];
	213	pid->ConfigureTPCdefaultCut(cutmodel, params,tpcdEdxcuthigh[0]);
	214	}
eba5ec88	215	/*
420bd2ea	216	char *cutmodel;
	217	cutmodel="pol0";
	218
	219	for(Int_t a=0;a<11;a++){
	220	// Not necessary anymore, since the pPb case is handled similarly to the pp case
	221	// cout << a << " " << paramsTPCdEdxcut[a] << endl;
	222	Double_t tpcparamlow[1]={paramsTPCdEdxcutlow[a]};
	223	Float_t tpcparamhigh=paramsTPCdEdxcuthigh[a];
	224	pid->ConfigureTPCcentralityCut(a,cutmodel,tpcparamlow,tpcparamhigh);
8dc7fa14	225	}
8dc7fa14	226	pid->ConfigureTPCdefaultCut(cutmodel,paramsTPCdEdxcutlow,paramsTPCdEdxcuthigh[0]); // After introducing the pPb flag, pPb is merged with pp and this line defines the cut
eba5ec88	227	*/
8dc7fa14	228
	229	// Configure TOF PID
	230	if (usetof){
	231	pid->ConfigureTOF(TOFs);
	232	AliHFEpidTOF *tofpid = pid->GetDetPID(AliHFEpid::kTOFpid);
	233	if (kTOFmis){
	234	tofpid->SetRejectTOFmismatch();
	235	}
	236	}
	237
	238	// To make different upper TOF cut to see contamination effect
	239	// The below two lines should be removed after this check
	240	//AliHFEpidTOF *tofpid = pid->GetDetPID(AliHFEpid::kTOFpid);
	241	//if(TOFs<3.) tofpid->SetTOFnSigmaBand(-3,TOFs); //only to check the assymmetric tof cut
420bd2ea	242
8dc7fa14	243	// Load hadron background
	244	if(!useMC){
	245	Bool_t status = kTRUE;
	246	TF1 *hBackground[12];
eba5ec88	247	if(isAOD==1) {
	248	if (usetof) status = ReadContaminationFunctions("hadroncontamination_AOD139_TOFPID_pPb_eta06.root", hBackground, tpcdEdxcutlow[0]);
	249	else status = ReadContaminationFunctions("hadroncontamination_AOD139_noTOFPID_pPb_eta06.root", hBackground, tpcdEdxcutlow[0]);
	250	}
39ebaf78	251	else if (isBeauty==1) {
	252	status = ReadContaminationFunctionsBeauty("hadroncontamination_ESD_Beauty_TOFPID_pPb_eta06.root", hBackground, tpcdEdxcutlow[0]);
	253	}
eba5ec88	254	else status = ReadContaminationFunctions("hadroncontamination_TOFTPC_pPb_eta06_newsplines_try3.root", hBackground, tpcdEdxcutlow[0]);
8dc7fa14	255	for(Int_t a=0;a<12;a++) {
420bd2ea	256	//printf("back %f \n",hBackground[a]);
8dc7fa14	257	if(status) task->SetBackGroundFactorsFunction(hBackground[a],a);
	258	else printf("not all background functions found\n");
	259	}
	260	}
	261
	262	//***************************************//
	263	// Configure NPE plugin //
	264	//***************************************//
420bd2ea	265
8dc7fa14	266	AliHFENonPhotonicElectron *backe = new AliHFENonPhotonicElectron(Form("HFEBackGroundSubtractionPID2%s",appendix.Data()),"Background subtraction"); //appendix
	267	//Setting the Cuts for the Associated electron-pool
	268	AliHFEcuts *hfeBackgroundCuts = new AliHFEcuts(Form("HFEBackSub%s",appendix.Data()),"Background sub Cuts");
	269	// hfeBackgroundCuts->SetEtaRange(assETA);
	270	hfeBackgroundCuts->SetEtaRange(assETAm,assETAp);
39ebaf78	271	hfeBackgroundCuts->SetPtRange(assMinPt,20.);
8dc7fa14	272
	273	hfeBackgroundCuts->SetMaxChi2perClusterTPC(4);
	274	hfeBackgroundCuts->SetMinNClustersITS(assITS);
	275	hfeBackgroundCuts->SetMinNClustersTPC(assTPCcl);
	276	hfeBackgroundCuts->SetMinNClustersTPCPID(assTPCPIDcl);
	277	hfeBackgroundCuts->SetMaxImpactParam(assDCAr,assDCAz);
	278	if(isAOD) hfeBackgroundCuts->SetAODFilterBit(4);
	279	hfeBackgroundCuts->SetQAOn(); // QA
	280
	281	AliHFEpid *pidbackground = backe->GetPIDBackground();
	282	if(useMC) pidbackground->SetHasMCData(kTRUE);
39ebaf78	283
	284	if (assTOFs>0.){
	285	pidbackground->AddDetector("TOF", 0);
	286	pidbackground->AddDetector("TPC", 1);
	287	} else {
	288	pidbackground->AddDetector("TPC", 0);
	289	}
	290
8dc7fa14	291	Double_t paramsTPCdEdxcutlowAssoc[12] ={-3.0,-3.0,-3.0,-3.0,-3.0,-3.0,-3.0,-3.0,-3.0,-3.0,-3.0,-3.0};
	292	if(assTPCSminus) memcpy(paramsTPCdEdxcutlowAssoc,assTPCSminus,sizeof(paramsTPCdEdxcutlowAssoc));
	293
	294	Double_t paramsTPCdEdxcuthighAssoc[12] ={3.0, 3.0, 3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0};
	295	if(assTPCSplus) memcpy(paramsTPCdEdxcuthighAssoc,assTPCSplus,sizeof(paramsTPCdEdxcuthighAssoc));
	296
	297	char *cutmodelAssoc;
	298	cutmodelAssoc="pol0";
420bd2ea	299	for(Int_t a=0;a<11;a++){
8dc7fa14	300	// Not necessary anymore, since the pPb case is handled similarly to the pp case
	301	// cout << a << " " << paramsTPCdEdxcut[a] << endl;
	302	Double_t tpcparamlow[1]={paramsTPCdEdxcutlowAssoc[a]};
	303	Float_t tpcparamhigh=paramsTPCdEdxcuthighAssoc[a];
	304	pidbackground->ConfigureTPCcentralityCut(a,cutmodelAssoc,tpcparamlow,tpcparamhigh);
	305	}
	306	pidbackground->ConfigureTPCdefaultCut(cutmodelAssoc,paramsTPCdEdxcutlowAssoc,paramsTPCdEdxcuthighAssoc[0]); // After introducing the pPb flag, pPb is merged with pp and this line defines the cut
b74766cf	307	//backe->GetPIDBackgroundQAManager()->SetHighResolutionHistos();
39ebaf78	308
	309	if (assTOFs>0.){
	310	pidbackground->ConfigureTOF(TOFs);
	311	}
	312
8dc7fa14	313	backe->SetHFEBackgroundCuts(hfeBackgroundCuts);
	314
	315	// Selection of associated tracks for the pool
	316	if(useCat1Tracks) backe->SelectCategory1Tracks(kTRUE);
1452614c	317	if(useCat2Tracks){
1452614c	318	backe->SelectCategory2Tracks(kTRUE);
39ebaf78	319	backe->SetITSMeanShift(-0.5);
	320	backe->SetITSnSigmaHigh(assITSpid);
	321	Double_t assITSminus = -1.0 * assITSpid;
	322	backe->SetITSnSigmaLow(assITSminus);
	323	//backe->SetminPt(assMinPt);
1452614c	324	}
8dc7fa14	325
	326	// apply opening angle cut to reduce file size
	327	backe->SetMaxInvMass(0.3);
1452614c	328	backe->SetPtBinning(sizept, ptbinning);
1452614c	329	backe->SetEtaBinning(sizeeta, etabinning);
39ebaf78	330	// MC weight
	331	if(useMC) {
	332	//printf("test put weight %d\n",weightlevelback);
	333	if((weightlevelback >=0) && (weightlevelback < 3)) backe->SetWithWeights(weightlevelback);
	334	}
8dc7fa14	335	task->SetHFEBackgroundSubtraction(backe);
8dc7fa14	336
39ebaf78	337	//task->SetWeightHist();
39ebaf78	338
8dc7fa14	339	//***************************************//
	340	// V0 tagged tracks //
	341	//***************************************//
	342
	343	if(kAnalyseTaggedTracks){
	344	AliHFEcuts *v0trackCuts = new AliHFEcuts("V0trackCuts", "Track Cuts for tagged track Analysis");
	345	v0trackCuts->CreateStandardCuts();
	346	v0trackCuts->SetMinNClustersTPC(TPCcl);
420bd2ea	347	v0trackCuts->SetMinNClustersTPCPID(TPCclPID);
8dc7fa14	348	v0trackCuts->SetMinRatioTPCclusters(0.6);
	349	v0trackCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
	350	v0trackCuts->SetMinNClustersITS(1);
	351	v0trackCuts->SetCutITSpixel(AliHFEextraCuts::kFirst);
	352	v0trackCuts->SetCheckITSLayerStatus(kFALSE);
	353	v0trackCuts->UnsetVertexRequirement();
	354	//hfecuts->SetSigmaToVertex(10);
420bd2ea	355	if(usetof) v0trackCuts->SetTOFPIDStep(kTRUE);
8dc7fa14	356	v0trackCuts->SetQAOn();
	357
	358	task->SwitchOnPlugin(AliAnalysisTaskHFE::kTaggedTrackAnalysis);
	359	task->SetTaggedTrackCuts(v0trackCuts);
	360	task->SetCleanTaggedTrack(kTRUE);
	361	}
	362
	363	// QA
	364	printf("task %p\n", task);
	365	task->SetQAOn(AliAnalysisTaskHFE::kPIDqa);
	366	task->SetQAOn(AliAnalysisTaskHFE::kMCqa);
8dc7fa14	367	task->SwitchOnPlugin(AliAnalysisTaskHFE::kDEstep);
39ebaf78	368	if(nonPhotonicElectronBeauty) task->SwitchOnPlugin(AliAnalysisTaskHFE::kNonPhotonicElectronBeauty);
39ebaf78	369	else task->SwitchOnPlugin(AliAnalysisTaskHFE::kNonPhotonicElectron);
8dc7fa14	370
	371	printf("*************************************\n");
	372	printf("Configuring standard Task:\n");
	373	task->PrintStatus();
	374	pid->PrintStatus();
	375	printf("*************************************\n");
	376	return task;
	377	}