//
// *** 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;
}