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