[u/mrichter/AliRoot.git] / PWGLF / RESONANCES / macros / train_old / RsnConfigPhi.C

//
// *** Configuration script for phi->KK analysis with 2010 runs ***
// 
// A configuration script for RSN package needs to define the followings:
//
// (1) decay tree of each resonance to be studied, which is needed to select
//     true pairs and to assign the right mass to all candidate daughters
// (2) cuts at all levels: single daughters, tracks, events
// (3) output objects: histograms or trees
//
Bool_t RsnConfigPhi
(
   Bool_t      isMC,
   const char *options,
   const char *path     = "$(ALICE_ROOT)/PWG2/RESONANCES/macros/train",
   const char *taskName = "RSNtask"
)
{
   // ==================================================================================================================
   // == OPTIONS =======================================================================================================
   // ==================================================================================================================
   
   // Instead or getting confused with plenty of arguments in the macro (with default values),
   // we use a unique string of options with a set of conventional strings to set up the job:
   // -- "MC"/"DATA" --> what kind of sample
   // -- "ITS"/"TPC" --> what tracks to use (ITS standalone and/or TPC+ITS)
   // -- "xxxPID"    --> add the PID cuts for the detector xxx.
   //
   // In this point, these options are converted into boolean variables.
   
   TString opt(options);
   opt.ToUpper();
   opt.ReplaceAll(" ", "");
   
   Bool_t addITS = opt.Contains("ITS");
   Bool_t addTPC = opt.Contains("TPC");
   Bool_t useITS = opt.Contains("ITSPID");
   Bool_t useTPC = opt.Contains("TPCPID");
   Bool_t useTOF = opt.Contains("TOFPID");
   
   // correct options when needed
   if (!addITS) useITS = kFALSE;
   if (!addTPC) useTPC = useTOF = kFALSE;
   
   // ==================================================================================================================
   // == DEFINITIONS ===================================================================================================
   // ==================================================================================================================
   
   // We put here the definitions of all objects which are needed in the following, in order to have then
   // a more readable code in the points where these objects are added to the analysis manager.
   
   // pair definitions --> decay channels:
   // in our case, unlike-charged KK pairs for the signal, and like-charged ones for background
   AliRsnPairDef *pairDef[3];
   pairDef[0] = new AliRsnPairDef(AliRsnDaughter::kKaon, '+', AliRsnDaughter::kKaon, '-', 333, 1.019455); // unlike
   pairDef[1] = new AliRsnPairDef(AliRsnDaughter::kKaon, '+', AliRsnDaughter::kKaon, '+', 333, 1.019455); // like ++
   pairDef[2] = new AliRsnPairDef(AliRsnDaughter::kKaon, '-', AliRsnDaughter::kKaon, '-', 333, 1.019455); // like --

   // computation objects:
   // these are the objects which drive the computations, whatever it is (histos or tree filling)
   // and all tracks passed to them will be given the mass according to the reference pair definition
   // we create two unlike-sign pair computators, one for all pairs and another for true pairs (useful in MC)
   AliRsnPairFunctions *pair[4];
   pair[0] = new AliRsnPairFunctions(Form("%s_kaonP_kaonM_phi", opt.Data()), pairDef[0]); // unlike - true
   pair[1] = new AliRsnPairFunctions(Form("%s_kaonP_kaonM_all", opt.Data()), pairDef[0]); // unlike - all
   pair[2] = new AliRsnPairFunctions(Form("%s_kaonP_kaonM_all", opt.Data()), pairDef[1]); // like ++
   pair[3] = new AliRsnPairFunctions(Form("%s_kaonP_kaonM_all", opt.Data()), pairDef[2]); // like --

   // set additional option for true pairs (slot [0])
   pair[0]->SetOnlyTrue(kTRUE);
   pair[0]->SetCheckDecay(kTRUE);
   
   // ==================================================================================================================
   // == CUTS AND AXES =================================================================================================
   // ==================================================================================================================
   
   //
   // Track quality for ITS standalone:
   // this cut is used to select tracks of good quality, irrespective of the PID.
   // When adding status flags, the second argument tells if each considered flag
   // must be active or not in the track status, since the ITS-SA tracks need that
   // some of them are OFF (e.g.: kTPCin)
   //
   AliRsnCutTrackQuality *cutQualityITS = new AliRsnCutTrackQuality("cutQualityITS");
   cutQualityITS->AddStatusFlag(AliESDtrack::kITSin    , kTRUE);
   cutQualityITS->AddStatusFlag(AliESDtrack::kTPCin    , kFALSE);
   cutQualityITS->AddStatusFlag(AliESDtrack::kITSrefit , kTRUE);
   cutQualityITS->AddStatusFlag(AliESDtrack::kTPCrefit , kFALSE);
   cutQualityITS->AddStatusFlag(AliESDtrack::kITSpureSA, kFALSE);
   cutQualityITS->AddStatusFlag(AliESDtrack::kITSpid   , kTRUE);
   cutQualityITS->SetPtRange(0.15, 1E+20);
   cutQualityITS->SetEtaRange(-0.8, 0.8);
   cutQualityITS->SetDCARPtFormula("0.0595+0.0182/pt^1.55");
   cutQualityITS->SetDCAZmax(2.0);
   cutQualityITS->SetSPDminNClusters(1);
   cutQualityITS->SetITSminNClusters(4);
   cutQualityITS->SetITSmaxChi2(2.0);
   cutQualityITS->SetTPCminNClusters(0);
   cutQualityITS->SetTPCmaxChi2(1E+10);
   cutQualityITS->SetRejectKinkDaughters();
      
   //
   // Track quality for TPC+ITS:
   // works exactly like the one above, but has settings for selecting TPC+ITS tracks
   // in this case, the flags required are all necessary, so here the procedure is simpler
   //
   AliRsnCutTrackQuality *cutQualityTPC = new AliRsnCutTrackQuality("cutQualityTPC");
   cutQualityTPC->AddStatusFlag(AliESDtrack::kTPCin   , kTRUE);
   cutQualityTPC->AddStatusFlag(AliESDtrack::kTPCrefit, kTRUE);
   cutQualityTPC->AddStatusFlag(AliESDtrack::kITSrefit, kTRUE);
   cutQualityTPC->SetPtRange(0.15, 1E+20);
   cutQualityTPC->SetEtaRange(-0.8, 0.8);
   cutQualityTPC->SetDCARPtFormula("0.0182+0.0350/pt^1.01");
   cutQualityTPC->SetDCAZmax(2.0);
   cutQualityTPC->SetSPDminNClusters(1);
   cutQualityTPC->SetITSminNClusters(0);
   cutQualityTPC->SetITSmaxChi2(1E+20);
   cutQualityTPC->SetTPCminNClusters(70);
   cutQualityTPC->SetTPCmaxChi2(4.0);
   cutQualityTPC->SetRejectKinkDaughters();
   
   //
   // ITS PID
   // In this implementation, it is a 3sigma cut around the Bethe-Bloch value.
   //
   // NOTE:
   // --> The initialization of the BB is different between data and MC.
   // --> The cut is the same for all momenta.
   //
   AliRsnCutPIDITS *cutPIDITSkaon = new AliRsnCutPIDITS("cutPIDITSkaon", AliPID::kKaon, -3.0, 3.0);
   cutPIDITSkaon->SetMC(isMC);
   cutPIDITSkaon->SetMomentumRange(0.0, 1E20);
   
   //
   // TPC PID
   // In this implementation, there are two instances:
   // - below 350 MeV --> 5 sigma cut
   // - above 350 MeV --> 3 sigma cut
   //
   // NOTE:
   // --> The initialization of the BB is different between data and MC.
   //
   AliRsnCutPIDTPC *cutPIDTPCkaonLow  = new AliRsnCutPIDTPC("cutPIDTPCkaonLow" , AliPID::kKaon, -5.0, 5.0);
   AliRsnCutPIDTPC *cutPIDTPCkaonHigh = new AliRsnCutPIDTPC("cutPIDTPCkaonHigh", AliPID::kKaon, -3.0, 3.0);
   
   // assign the momentum range and tell to reject tracks outside it
   cutPIDTPCkaonLow ->SetMomentumRange(0.00, 0.35);
   cutPIDTPCkaonHigh->SetMomentumRange(0.35, 1E20);
   cutPIDTPCkaonLow ->SetRejectOutside(kTRUE);
   cutPIDTPCkaonHigh->SetRejectOutside(kTRUE);
   
   // BB parameterization depends on data sample (MC, data)
   // the momentum range is passed and tracks outside it are rejected
   Double_t bbPar[5];
   if (isMC) {
      bbPar[0] = 2.15898 / 50.0;
      bbPar[1] = 1.75295E1;
      bbPar[2] = 3.40030E-9;
      bbPar[3] = 1.96178;
      bbPar[4] = 3.91720;
   } else {
      bbPar[0] = 1.41543 / 50.0;
      bbPar[1] = 2.63394E1;
      bbPar[2] = 5.0411E-11;
      bbPar[3] = 2.12543;
      bbPar[4] = 4.88663;
   }
   cutPIDTPCkaonLow ->SetBBParam(bbPar);
   cutPIDTPCkaonHigh->SetBBParam(bbPar);
   
   //
   // TOF PID
   // In this implementation it is a 3sigma cout aroung expected kaon time.
   //
   // NOTE:
   // --> It is important to choose if this cut must reject tracks not matched in TOF.
   //     Usually, if TPC pid is used, we can accept them, otherwise we must reject.
   //     (here we assume TPC is used)
   //
   AliRsnCutPIDTOF *cutPIDTOFkaon = new AliRsnCutPIDTOF("cutPIDTOFkaon", AliPID::kKaon, -3.0, 3.0);
   cutPIDTOFkaon->SetRejectUnmatched(!useTPC);
   
   //
   // Rapidity cut
   // Only thing to consider is that it needs a support object to define mass
   //
   AliRsnCutValue *cutRapidity = new AliRsnCutValue("cutY", AliRsnValue::kPairY, -0.5, 0.5);
   cutRapidity->GetValueObj()->SetSupportObject(pairDef[0]);
   
   //
   // Axes
   //
   // NOTE:
   // --> multiplicity has variable bins, defined by array below
   //
   
   Double_t mult[] = { 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12., 13.,  14.,  15.,  16.,  17.,  18.,  19., 
                      20., 21., 22., 23., 24., 25., 30., 35., 40., 50., 60., 70., 80., 90., 100., 120., 140., 160., 180., 200., 500.};
   Int_t    nmult  = sizeof(mult) / sizeof(mult[0]);
   
   AliRsnValue *axisIM      = new AliRsnValue("IM"  , AliRsnValue::kPairInvMass   ,  0.9, 1.4, 0.001);
   AliRsnValue *axisRes     = new AliRsnValue("Res" , AliRsnValue::kPairInvMassRes, -0.5, 0.5, 0.001);
   AliRsnValue *axisPt      = new AliRsnValue("PT"  , AliRsnValue::kPairPt        ,  0.0, 5.0, 0.1  );
   AliRsnValue *axisY       = new AliRsnValue("Y"   , AliRsnValue::kPairY         , -1.1, 1.1, 0.1  );
   AliRsnValue *axisMultESD = new AliRsnValue("MESD", AliRsnValue::kEventMultESDCuts, nmult, mult);
   AliRsnValue *axisMultSPD = new AliRsnValue("MSPD", AliRsnValue::kEventMultSPD    , nmult, mult);
   AliRsnValue *axisMultMC  = new AliRsnValue("MMC" , AliRsnValue::kEventMultMC     , nmult, mult);
   
   // ==================================================================================================================
   // == PRELIMINARY OPERATIONS ========================================================================================
   // ==================================================================================================================

   // retrieve task from manager, using its name
   AliAnalysisManager *mgr  = AliAnalysisManager::GetAnalysisManager();
   AliRsnAnalysisTask *task = (AliRsnAnalysisTask*)mgr->GetTask(taskName);
   if (!task) {
      Error("RsnConfigMonitor", "Task not found");
      return kFALSE;
   }
   
   // ==================================================================================================================
   // == SINGLE DAUGHTER CUTS ==========================================================================================
   // ==================================================================================================================
   
   // we need to define cuts for track quality and for PID
   // which one will be used, depend on the analysis type
   
   Char_t qualityITS[255], qualityTPC[255];
   Char_t pidITS[255], pidTPC[255], pidTOF[255];
   Char_t schemeITS[255], schemeTPC[255], scheme[255];
   
   sprintf(qualityITS, "%s"     , cutQualityITS->GetName());
   sprintf(qualityTPC, "%s"     , cutQualityTPC->GetName());
   sprintf(pidITS    , "%s"     , cutPIDITSkaon->GetName());
   sprintf(pidTPC    , "(%s|%s)", cutPIDTPCkaonHigh->GetName(), cutPIDTPCkaonLow->GetName());
   sprintf(pidTOF    , "%s"     , cutPIDTOFkaon->GetName());
   sprintf(schemeITS , "");
   sprintf(schemeTPC , "");
   
   // choose cuts to add depending on used tracks
   for (Int_t i = 0; i < 4; i++) {
      // assign cuts to common manager for daughters
      AliRsnCutSet *cutSet = pair[i]->GetCutManager()->GetCommonDaughterCuts();
      // add cuts and define schemes depending on options
      if (addITS) {
         // if adding ITS standalone,
         // by default we use the quality check
         // and if required we add the PID (option "ITSPID")
         cutSet->AddCut(cutQualityITS);
         sprintf(schemeITS, "%s", qualityITS);
         if (useITS) {
            sprintf(schemeITS, "%s & %s", qualityITS, pidITS);
            cutSet->AddCut(cutPIDITSkaon);
         }
      }
      if (addTPC) {
         // if adding TPC+ITS tracks,
         // by default we use the quality check
         // and then we can use no PID, or both, 
         // or only one among TPC and TOF
         // the TPC PID cut is always an 'OR' of 
         // the two ones, in order to cover the full momentum range
         cutSet->AddCut(cutQualityTPC);
         sprintf(schemeTPC, "%s", qualityTPC);
         if (useTPC && useTOF) {
            cutSet->AddCut(cutPIDTPCkaonLow);
            cutSet->AddCut(cutPIDTPCkaonHigh);
            cutSet->AddCut(cutPIDTOFkaon);
            sprintf(schemeTPC, "%s & %s & %s", qualityTPC, pidTPC, pidTOF);
         } else if (useTPC) {
            cutSet->AddCut(cutPIDTPCkaonLow);
            cutSet->AddCut(cutPIDTPCkaonHigh);
            sprintf(schemeTPC, "%s & %s", qualityTPC, pidTPC);
         } else if (useTOF) {
            cutSet->AddCut(cutPIDTOFkaon);
            sprintf(schemeTPC, "%s & %s", qualityTPC, pidTOF);
         }
      }
      
      // final scheme depends on what of the above were added
      // in case both ITS-SA and TPC tracks are added, the scheme
      // is the OR of the cuts for the first and those for the second
      // category of tracks
      if (strlen(schemeITS) > 0 && strlen(schemeTPC) > 0) {
         sprintf(scheme, "(%s) | (%s)", schemeITS, schemeTPC);
      } else if (strlen(schemeITS) > 0) {
         sprintf(scheme, "%s", schemeITS);
      } else if (strlen(schemeTPC) > 0) {
         sprintf(scheme, "%s", schemeTPC);
      } else {
         ::Error("Scheme is empty!");
         return kFALSE;
      }
      cutSet->SetCutScheme(scheme);
      ::Info("RsnConfigPhi", "Scheme for daughter cuts: %s", cutSet->GetCutScheme().Data());
   }
      
   // ==================================================================================================================
   // == PAIR CUTS =====================================================================================================
   // ==================================================================================================================

   // in this case, we add the cut to the specific cut sets of all pairs
   // and we must then loop over all pairs, to add cuts to the related sets
   for (Int_t ipair = 0; ipair < 4; ipair++) {
      pair[ipair]->GetCutManager()->GetMotherCuts()->AddCut(cutRapidity);
      pair[ipair]->GetCutManager()->GetMotherCuts()->SetCutScheme(cutRapidity->GetName());
      ::Info("RsnConfigPhi", "Scheme for pair cuts: %s", pair[ipair]->GetCutManager()->GetMotherCuts()->GetCutScheme().Data());
   }
   
   // ==================================================================================================================
   // == FUNCTIONS FOR HISTOGRAMS ======================================================================================
   // ==================================================================================================================
   
   // we define a histogram with inv mass and other support binning values
   // and, when possible, a resolution on inv. mass

   // create function for inv. mass and add axes
   AliRsnFunction *fcnIM = new AliRsnFunction;
   if (!fcnIM->AddAxis(axisIM)     ) return kFALSE;
   if (!fcnIM->AddAxis(axisPt)     ) return kFALSE;
   if (!fcnIM->AddAxis(axisMultSPD)) return kFALSE;
   fcnIM->UseSparse();
   
   // create function for inv. mass and add axes
   AliRsnFunction *fcnIMRes = new AliRsnFunction;
   if (!fcnIMRes->AddAxis(axisIM)     ) return kFALSE;
   if (!fcnIMRes->AddAxis(axisPt)     ) return kFALSE;
   if (!fcnIMRes->AddAxis(axisMultSPD)) return kFALSE;
   if (!fcnIMRes->AddAxis(axisRes)    ) return kFALSE;
   
   // add functions to pairs
   pair[0]->AddFunction(fcnIMRes);
   for (Int_t ipair = 1; ipair < 4; ipair++) pair[ipair]->AddFunction(fcnIM);
   
   // ==================================================================================================================
   // == CONCLUSION ====================================================================================================
   // ==================================================================================================================
   
   // add all created AliRsnPair objects to the AliRsnAnalysisManager in the task
   task->GetAnalysisManager()->Add(pair[1]);
   task->GetAnalysisManager()->Add(pair[2]);
   task->GetAnalysisManager()->Add(pair[3]);
   if (isMC) task->GetAnalysisManager()->Add(pair[0]);
   
   return kTRUE;
}
Commit	Line	Data
4873fa4a	1	//
	2	// * Configuration script for phi->KK analysis with 2010 runs *
	3	//
	4	// A configuration script for RSN package needs to define the followings:
	5	//
	6	// (1) decay tree of each resonance to be studied, which is needed to select
	7	// true pairs and to assign the right mass to all candidate daughters
	8	// (2) cuts at all levels: single daughters, tracks, events
	9	// (3) output objects: histograms or trees
	10	//
	11	Bool_t RsnConfigPhi
	12	(
	13	Bool_t isMC,
	14	const char *options,
	15	const char *path = "$(ALICE_ROOT)/PWG2/RESONANCES/macros/train",
	16	const char *taskName = "RSNtask"
	17	)
	18	{
4873fa4a	19	// ==================================================================================================================
	20	// == OPTIONS =======================================================================================================
	21	// ==================================================================================================================
	22
	23	// Instead or getting confused with plenty of arguments in the macro (with default values),
	24	// we use a unique string of options with a set of conventional strings to set up the job:
	25	// -- "MC"/"DATA" --> what kind of sample
	26	// -- "ITS"/"TPC" --> what tracks to use (ITS standalone and/or TPC+ITS)
	27	// -- "xxxPID" --> add the PID cuts for the detector xxx.
	28	//
	29	// In this point, these options are converted into boolean variables.
	30
	31	TString opt(options);
	32	opt.ToUpper();
	33	opt.ReplaceAll(" ", "");
	34
	35	Bool_t addITS = opt.Contains("ITS");
	36	Bool_t addTPC = opt.Contains("TPC");
	37	Bool_t useITS = opt.Contains("ITSPID");
	38	Bool_t useTPC = opt.Contains("TPCPID");
	39	Bool_t useTOF = opt.Contains("TOFPID");
	40
	41	// correct options when needed
	42	if (!addITS) useITS = kFALSE;
	43	if (!addTPC) useTPC = useTOF = kFALSE;
	44
	45	// ==================================================================================================================
	46	// == DEFINITIONS ===================================================================================================
	47	// ==================================================================================================================
	48
	49	// We put here the definitions of all objects which are needed in the following, in order to have then
	50	// a more readable code in the points where these objects are added to the analysis manager.
	51
	52	// pair definitions --> decay channels:
	53	// in our case, unlike-charged KK pairs for the signal, and like-charged ones for background
	54	AliRsnPairDef *pairDef[3];
	55	pairDef[0] = new AliRsnPairDef(AliRsnDaughter::kKaon, '+', AliRsnDaughter::kKaon, '-', 333, 1.019455); // unlike
	56	pairDef[1] = new AliRsnPairDef(AliRsnDaughter::kKaon, '+', AliRsnDaughter::kKaon, '+', 333, 1.019455); // like ++
	57	pairDef[2] = new AliRsnPairDef(AliRsnDaughter::kKaon, '-', AliRsnDaughter::kKaon, '-', 333, 1.019455); // like --
	58
	59	// computation objects:
	60	// these are the objects which drive the computations, whatever it is (histos or tree filling)
	61	// and all tracks passed to them will be given the mass according to the reference pair definition
	62	// we create two unlike-sign pair computators, one for all pairs and another for true pairs (useful in MC)
	63	AliRsnPairFunctions *pair[4];
	64	pair[0] = new AliRsnPairFunctions(Form("%s_kaonP_kaonM_phi", opt.Data()), pairDef[0]); // unlike - true
	65	pair[1] = new AliRsnPairFunctions(Form("%s_kaonP_kaonM_all", opt.Data()), pairDef[0]); // unlike - all
	66	pair[2] = new AliRsnPairFunctions(Form("%s_kaonP_kaonM_all", opt.Data()), pairDef[1]); // like ++
	67	pair[3] = new AliRsnPairFunctions(Form("%s_kaonP_kaonM_all", opt.Data()), pairDef[2]); // like --
	68
	69	// set additional option for true pairs (slot [0])
	70	pair[0]->SetOnlyTrue(kTRUE);
	71	pair[0]->SetCheckDecay(kTRUE);
	72
f4459279	73	// ==================================================================================================================
	74	// == CUTS AND AXES =================================================================================================
	75	// ==================================================================================================================
	76
	77	//
	78	// Track quality for ITS standalone:
	79	// this cut is used to select tracks of good quality, irrespective of the PID.
	80	// When adding status flags, the second argument tells if each considered flag
	81	// must be active or not in the track status, since the ITS-SA tracks need that
	82	// some of them are OFF (e.g.: kTPCin)
	83	//
	84	AliRsnCutTrackQuality *cutQualityITS = new AliRsnCutTrackQuality("cutQualityITS");
	85	cutQualityITS->AddStatusFlag(AliESDtrack::kITSin , kTRUE);
	86	cutQualityITS->AddStatusFlag(AliESDtrack::kTPCin , kFALSE);
	87	cutQualityITS->AddStatusFlag(AliESDtrack::kITSrefit , kTRUE);
	88	cutQualityITS->AddStatusFlag(AliESDtrack::kTPCrefit , kFALSE);
	89	cutQualityITS->AddStatusFlag(AliESDtrack::kITSpureSA, kFALSE);
	90	cutQualityITS->AddStatusFlag(AliESDtrack::kITSpid , kTRUE);
	91	cutQualityITS->SetPtRange(0.15, 1E+20);
	92	cutQualityITS->SetEtaRange(-0.8, 0.8);
	93	cutQualityITS->SetDCARPtFormula("0.0595+0.0182/pt^1.55");
	94	cutQualityITS->SetDCAZmax(2.0);
	95	cutQualityITS->SetSPDminNClusters(1);
	96	cutQualityITS->SetITSminNClusters(4);
	97	cutQualityITS->SetITSmaxChi2(2.0);
	98	cutQualityITS->SetTPCminNClusters(0);
	99	cutQualityITS->SetTPCmaxChi2(1E+10);
	100	cutQualityITS->SetRejectKinkDaughters();
	101
	102	//
	103	// Track quality for TPC+ITS:
	104	// works exactly like the one above, but has settings for selecting TPC+ITS tracks
	105	// in this case, the flags required are all necessary, so here the procedure is simpler
	106	//
	107	AliRsnCutTrackQuality *cutQualityTPC = new AliRsnCutTrackQuality("cutQualityTPC");
	108	cutQualityTPC->AddStatusFlag(AliESDtrack::kTPCin , kTRUE);
	109	cutQualityTPC->AddStatusFlag(AliESDtrack::kTPCrefit, kTRUE);
	110	cutQualityTPC->AddStatusFlag(AliESDtrack::kITSrefit, kTRUE);
	111	cutQualityTPC->SetPtRange(0.15, 1E+20);
	112	cutQualityTPC->SetEtaRange(-0.8, 0.8);
	113	cutQualityTPC->SetDCARPtFormula("0.0182+0.0350/pt^1.01");
	114	cutQualityTPC->SetDCAZmax(2.0);
	115	cutQualityTPC->SetSPDminNClusters(1);
	116	cutQualityTPC->SetITSminNClusters(0);
	117	cutQualityTPC->SetITSmaxChi2(1E+20);
	118	cutQualityTPC->SetTPCminNClusters(70);
	119	cutQualityTPC->SetTPCmaxChi2(4.0);
	120	cutQualityTPC->SetRejectKinkDaughters();
	121
	122	//
	123	// ITS PID
	124	// In this implementation, it is a 3sigma cut around the Bethe-Bloch value.
	125	//
	126	// NOTE:
	127	// --> The initialization of the BB is different between data and MC.
	128	// --> The cut is the same for all momenta.
	129	//
	130	AliRsnCutPIDITS *cutPIDITSkaon = new AliRsnCutPIDITS("cutPIDITSkaon", AliPID::kKaon, -3.0, 3.0);
	131	cutPIDITSkaon->SetMC(isMC);
	132	cutPIDITSkaon->SetMomentumRange(0.0, 1E20);
	133
	134	//
	135	// TPC PID
	136	// In this implementation, there are two instances:
137	// - below 350 MeV --> 5 sigma cut
138	// - above 350 MeV --> 3 sigma cut
139	//
140	// NOTE:
141	// --> The initialization of the BB is different between data and MC.
142	//
143	AliRsnCutPIDTPC *cutPIDTPCkaonLow = new AliRsnCutPIDTPC("cutPIDTPCkaonLow" , AliPID::kKaon, -5.0, 5.0);
144	AliRsnCutPIDTPC *cutPIDTPCkaonHigh = new AliRsnCutPIDTPC("cutPIDTPCkaonHigh", AliPID::kKaon, -3.0, 3.0);
145
146	// assign the momentum range and tell to reject tracks outside it
147	cutPIDTPCkaonLow ->SetMomentumRange(0.00, 0.35);
148	cutPIDTPCkaonHigh->SetMomentumRange(0.35, 1E20);
149	cutPIDTPCkaonLow ->SetRejectOutside(kTRUE);
150	cutPIDTPCkaonHigh->SetRejectOutside(kTRUE);
151
152	// BB parameterization depends on data sample (MC, data)
153	// the momentum range is passed and tracks outside it are rejected
154	Double_t bbPar[5];
155	if (isMC) {
156	bbPar[0] = 2.15898 / 50.0;
157	bbPar[1] = 1.75295E1;
158	bbPar[2] = 3.40030E-9;
159	bbPar[3] = 1.96178;
160	bbPar[4] = 3.91720;
161	} else {
162	bbPar[0] = 1.41543 / 50.0;
163	bbPar[1] = 2.63394E1;
164	bbPar[2] = 5.0411E-11;
165	bbPar[3] = 2.12543;
166	bbPar[4] = 4.88663;
167	}
168	cutPIDTPCkaonLow ->SetBBParam(bbPar);
169	cutPIDTPCkaonHigh->SetBBParam(bbPar);
170
171	//
172	// TOF PID
173	// In this implementation it is a 3sigma cout aroung expected kaon time.
174	//
175	// NOTE:
176	// --> It is important to choose if this cut must reject tracks not matched in TOF.
177	// Usually, if TPC pid is used, we can accept them, otherwise we must reject.
178	// (here we assume TPC is used)
179	//
180	AliRsnCutPIDTOF *cutPIDTOFkaon = new AliRsnCutPIDTOF("cutPIDTOFkaon", AliPID::kKaon, -3.0, 3.0);
181	cutPIDTOFkaon->SetRejectUnmatched(!useTPC);
182
183	//
184	// Rapidity cut
185	// Only thing to consider is that it needs a support object to define mass
186	//
187	AliRsnCutValue *cutRapidity = new AliRsnCutValue("cutY", AliRsnValue::kPairY, -0.5, 0.5);
188	cutRapidity->GetValueObj()->SetSupportObject(pairDef[0]);
189
190	//
191	// Axes
192	//
193	// NOTE:
194	// --> multiplicity has variable bins, defined by array below
195	//
196
197	Double_t mult[] = { 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19.,
198	20., 21., 22., 23., 24., 25., 30., 35., 40., 50., 60., 70., 80., 90., 100., 120., 140., 160., 180., 200., 500.};
199	Int_t nmult = sizeof(mult) / sizeof(mult[0]);
200
201	AliRsnValue *axisIM = new AliRsnValue("IM" , AliRsnValue::kPairInvMass , 0.9, 1.4, 0.001);
202	AliRsnValue *axisRes = new AliRsnValue("Res" , AliRsnValue::kPairInvMassRes, -0.5, 0.5, 0.001);
203	AliRsnValue *axisPt = new AliRsnValue("PT" , AliRsnValue::kPairPt , 0.0, 5.0, 0.1 );
204	AliRsnValue *axisY = new AliRsnValue("Y" , AliRsnValue::kPairY , -1.1, 1.1, 0.1 );
30447d77	205	AliRsnValue *axisMultESD = new AliRsnValue("MESD", AliRsnValue::kEventMultESDCuts, nmult, mult);
	206	AliRsnValue *axisMultSPD = new AliRsnValue("MSPD", AliRsnValue::kEventMultSPD , nmult, mult);
	207	AliRsnValue *axisMultMC = new AliRsnValue("MMC" , AliRsnValue::kEventMultMC , nmult, mult);
f4459279	208
	209	// ==================================================================================================================
	210	// == PRELIMINARY OPERATIONS ========================================================================================
	211	// ==================================================================================================================
	212
	213	// retrieve task from manager, using its name
	214	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
	215	AliRsnAnalysisTask task = (AliRsnAnalysisTask)mgr->GetTask(taskName);
	216	if (!task) {
	217	Error("RsnConfigMonitor", "Task not found");
	218	return kFALSE;
	219	}
4873fa4a	220
	221	// ==================================================================================================================
	222	// == SINGLE DAUGHTER CUTS ==========================================================================================
	223	// ==================================================================================================================
	224
	225	// we need to define cuts for track quality and for PID
	226	// which one will be used, depend on the analysis type
	227
4873fa4a	228	Char_t qualityITS[255], qualityTPC[255];
	229	Char_t pidITS[255], pidTPC[255], pidTOF[255];
	230	Char_t schemeITS[255], schemeTPC[255], scheme[255];
	231
	232	sprintf(qualityITS, "%s" , cutQualityITS->GetName());
	233	sprintf(qualityTPC, "%s" , cutQualityTPC->GetName());
	234	sprintf(pidITS , "%s" , cutPIDITSkaon->GetName());
	235	sprintf(pidTPC , "(%s\|%s)", cutPIDTPCkaonHigh->GetName(), cutPIDTPCkaonLow->GetName());
	236	sprintf(pidTOF , "%s" , cutPIDTOFkaon->GetName());
	237	sprintf(schemeITS , "");
	238	sprintf(schemeTPC , "");
	239
	240	// choose cuts to add depending on used tracks
	241	for (Int_t i = 0; i < 4; i++) {
	242	// assign cuts to common manager for daughters
	243	AliRsnCutSet *cutSet = pair[i]->GetCutManager()->GetCommonDaughterCuts();
	244	// add cuts and define schemes depending on options
	245	if (addITS) {
	246	// if adding ITS standalone,
	247	// by default we use the quality check
	248	// and if required we add the PID (option "ITSPID")
	249	cutSet->AddCut(cutQualityITS);
	250	sprintf(schemeITS, "%s", qualityITS);
	251	if (useITS) {
	252	sprintf(schemeITS, "%s & %s", qualityITS, pidITS);
	253	cutSet->AddCut(cutPIDITSkaon);
	254	}
	255	}
	256	if (addTPC) {
	257	// if adding TPC+ITS tracks,
	258	// by default we use the quality check
	259	// and then we can use no PID, or both,
	260	// or only one among TPC and TOF
	261	// the TPC PID cut is always an 'OR' of
	262	// the two ones, in order to cover the full momentum range
	263	cutSet->AddCut(cutQualityTPC);
	264	sprintf(schemeTPC, "%s", qualityTPC);
	265	if (useTPC && useTOF) {
	266	cutSet->AddCut(cutPIDTPCkaonLow);
	267	cutSet->AddCut(cutPIDTPCkaonHigh);
	268	cutSet->AddCut(cutPIDTOFkaon);
	269	sprintf(schemeTPC, "%s & %s & %s", qualityTPC, pidTPC, pidTOF);
	270	} else if (useTPC) {
	271	cutSet->AddCut(cutPIDTPCkaonLow);
	272	cutSet->AddCut(cutPIDTPCkaonHigh);
	273	sprintf(schemeTPC, "%s & %s", qualityTPC, pidTPC);
	274	} else if (useTOF) {
	275	cutSet->AddCut(cutPIDTOFkaon);
	276	sprintf(schemeTPC, "%s & %s", qualityTPC, pidTOF);
	277	}
	278	}
	279
	280	// final scheme depends on what of the above were added
	281	// in case both ITS-SA and TPC tracks are added, the scheme
	282	// is the OR of the cuts for the first and those for the second
	283	// category of tracks
	284	if (strlen(schemeITS) > 0 && strlen(schemeTPC) > 0) {
	285	sprintf(scheme, "(%s) \| (%s)", schemeITS, schemeTPC);
	286	} else if (strlen(schemeITS) > 0) {
	287	sprintf(scheme, "%s", schemeITS);
	288	} else if (strlen(schemeTPC) > 0) {
	289	sprintf(scheme, "%s", schemeTPC);
	290	} else {
	291	::Error("Scheme is empty!");
292	return kFALSE;
293	}
294	cutSet->SetCutScheme(scheme);
295	::Info("RsnConfigPhi", "Scheme for daughter cuts: %s", cutSet->GetCutScheme().Data());
296	}
297
298	// ==================================================================================================================
299	// == PAIR CUTS =====================================================================================================
300	// ==================================================================================================================
4873fa4a	301
	302	// in this case, we add the cut to the specific cut sets of all pairs
	303	// and we must then loop over all pairs, to add cuts to the related sets
	304	for (Int_t ipair = 0; ipair < 4; ipair++) {
	305	pair[ipair]->GetCutManager()->GetMotherCuts()->AddCut(cutRapidity);
	306	pair[ipair]->GetCutManager()->GetMotherCuts()->SetCutScheme(cutRapidity->GetName());
	307	::Info("RsnConfigPhi", "Scheme for pair cuts: %s", pair[ipair]->GetCutManager()->GetMotherCuts()->GetCutScheme().Data());
	308	}
	309
	310	// ==================================================================================================================
	311	// == FUNCTIONS FOR HISTOGRAMS ======================================================================================
	312	// ==================================================================================================================
	313
	314	// we define a histogram with inv mass and other support binning values
	315	// and, when possible, a resolution on inv. mass
4873fa4a	316
	317	// create function for inv. mass and add axes
	318	AliRsnFunction *fcnIM = new AliRsnFunction;
	319	if (!fcnIM->AddAxis(axisIM) ) return kFALSE;
	320	if (!fcnIM->AddAxis(axisPt) ) return kFALSE;
	321	if (!fcnIM->AddAxis(axisMultSPD)) return kFALSE;
	322	fcnIM->UseSparse();
	323
	324	// create function for inv. mass and add axes
	325	AliRsnFunction *fcnIMRes = new AliRsnFunction;
	326	if (!fcnIMRes->AddAxis(axisIM) ) return kFALSE;
	327	if (!fcnIMRes->AddAxis(axisPt) ) return kFALSE;
	328	if (!fcnIMRes->AddAxis(axisMultSPD)) return kFALSE;
	329	if (!fcnIMRes->AddAxis(axisRes) ) return kFALSE;
	330
	331	// add functions to pairs
	332	pair[0]->AddFunction(fcnIMRes);
	333	for (Int_t ipair = 1; ipair < 4; ipair++) pair[ipair]->AddFunction(fcnIM);
	334
	335	// ==================================================================================================================
	336	// == CONCLUSION ====================================================================================================
	337	// ==================================================================================================================
	338
	339	// add all created AliRsnPair objects to the AliRsnAnalysisManager in the task
	340	task->GetAnalysisManager()->Add(pair[1]);
	341	task->GetAnalysisManager()->Add(pair[2]);
	342	task->GetAnalysisManager()->Add(pair[3]);
	343	if (isMC) task->GetAnalysisManager()->Add(pair[0]);
	344
	345	return kTRUE;
	346	}