]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliCFTaskVertexingHF.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coveriry
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliCFTaskVertexingHF.cxx
CommitLineData
379592af 1/**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//-----------------------------------------------------------------------
17// Class for HF corrections as a function of many variables
18// 6 Steps introduced: MC, MC Acc, Reco, Reco Acc, Reco Acc + ITS Cl,
19// Reco Acc + ITS Cl + PPR cuts
20// 12 variables used: pt, y, cosThetaStar, ptPi, ptK, ct,
21// dca, d0Pi, d0K, d0Pixd0K, cosPointingAngle, phi
22//
23//-----------------------------------------------------------------------
24// Author : C. Zampolli, CERN
25// D. Caffarri, Univ & INFN Padova caffarri@pd.infn.it
26//-----------------------------------------------------------------------
27//-----------------------------------------------------------------------
28// Base class for HF Unfolding (pt and eta)
29// correlation matrix filled at Acceptance and PPR level
30// Author: A.Grelli , Utrecht - agrelli@uu.nl
31//-----------------------------------------------------------------------
32#include <TCanvas.h>
33#include <TParticle.h>
34#include <TDatabasePDG.h>
35#include <TH1I.h>
36#include <TStyle.h>
37#include <TFile.h>
38
39#include "AliCFTaskVertexingHF.h"
40#include "AliStack.h"
41#include "AliMCEvent.h"
42#include "AliCFManager.h"
43#include "AliCFContainer.h"
44#include "AliLog.h"
45#include "AliAnalysisManager.h"
46#include "AliAODHandler.h"
47#include "AliAODEvent.h"
48#include "AliAODRecoDecay.h"
49#include "AliAODRecoDecayHF.h"
50#include "AliAODRecoDecayHF2Prong.h"
f2703bd2 51#include "AliAODRecoDecayHF3Prong.h"
52#include "AliAODRecoDecayHF4Prong.h"
53#include "AliAODRecoCascadeHF.h"
379592af 54#include "AliAODMCParticle.h"
55#include "AliAODMCHeader.h"
56#include "AliESDtrack.h"
57#include "TChain.h"
58#include "THnSparse.h"
59#include "TH2D.h"
60#include "AliESDtrackCuts.h"
f2703bd2 61#include "AliRDHFCuts.h"
379592af 62#include "AliRDHFCutsD0toKpi.h"
f2703bd2 63#include "AliRDHFCutsDplustoKpipi.h"
64#include "AliRDHFCutsDStartoKpipi.h"
65#include "AliRDHFCutsDstoKKpi.h"
66#include "AliRDHFCutsLctopKpi.h"
67#include "AliRDHFCutsD0toKpipipi.h"
379592af 68#include "AliCFVertexingHF2Prong.h"
4e600f58 69#include "AliCFVertexingHF3Prong.h"
367e9aa3 70#include "AliCFVertexingHFCascade.h"
379592af 71#include "AliCFVertexingHF.h"
f2703bd2 72#include "AliAnalysisDataSlot.h"
73#include "AliAnalysisDataContainer.h"
379592af 74
75//__________________________________________________________________________
76AliCFTaskVertexingHF::AliCFTaskVertexingHF() :
f2703bd2 77 AliAnalysisTaskSE(),
f2703bd2 78 fCFManager(0x0),
79 fHistEventsProcessed(0x0),
80 fCorrelation(0x0),
81 fCountMC(0),
82 fCountAcc(0),
83 fCountVertex(0),
84 fCountRefit(0),
85 fCountReco(0),
86 fCountRecoAcc(0),
87 fCountRecoITSClusters(0),
88 fCountRecoPPR(0),
89 fCountRecoPID(0),
90 fEvents(0),
91 fDecayChannel(0),
92 fFillFromGenerated(kFALSE),
93 fOriginDselection(0),
94 fAcceptanceUnf(kTRUE),
95 fCuts(0),
96 fUseWeight(kFALSE),
97 fWeight(1.),
31da6b05 98 fNvar(0),
99 fPartName(""),
3ee5eb83 100 fDauNames(""),
b7af2639 101 fSign(2),
fbec9fa9 102 fCentralitySelection(kTRUE),
0ada222f 103 fFakeSelection(0),
104 fRejectIfNoQuark(kTRUE),
6c62cb59 105 fUseMCVertex(kFALSE),
ec5288c3 106 fDsOption(1),
107 fConfiguration(kCheetah) // by default, setting the fast configuration
379592af 108{
109 //
110 //Default ctor
111 //
112}
113//___________________________________________________________________________
f2703bd2 114AliCFTaskVertexingHF::AliCFTaskVertexingHF(const Char_t* name, AliRDHFCuts* cuts) :
115 AliAnalysisTaskSE(name),
f2703bd2 116 fCFManager(0x0),
117 fHistEventsProcessed(0x0),
118 fCorrelation(0x0),
119 fCountMC(0),
120 fCountAcc(0),
121 fCountVertex(0),
122 fCountRefit(0),
123 fCountReco(0),
124 fCountRecoAcc(0),
125 fCountRecoITSClusters(0),
126 fCountRecoPPR(0),
127 fCountRecoPID(0),
128 fEvents(0),
129 fDecayChannel(0),
130 fFillFromGenerated(kFALSE),
131 fOriginDselection(0),
132 fAcceptanceUnf(kTRUE),
133 fCuts(cuts),
134 fUseWeight(kFALSE),
135 fWeight(1.),
31da6b05 136 fNvar(0),
137 fPartName(""),
3ee5eb83 138 fDauNames(""),
b7af2639 139 fSign(2),
fbec9fa9 140 fCentralitySelection(kTRUE),
0ada222f 141 fFakeSelection(0),
142 fRejectIfNoQuark(kTRUE),
6c62cb59 143 fUseMCVertex(kFALSE),
ec5288c3 144 fDsOption(1),
145 fConfiguration(kCheetah) // by default, setting the fast configuration
379592af 146{
147 //
148 // Constructor. Initialization of Inputs and Outputs
149 //
379592af 150 /*
f2703bd2 151 DefineInput(0) and DefineOutput(0)
152 are taken care of by AliAnalysisTaskSE constructor
153 */
379592af 154 DefineOutput(1,TH1I::Class());
155 DefineOutput(2,AliCFContainer::Class());
156 DefineOutput(3,THnSparseD::Class());
f2703bd2 157 DefineOutput(4,AliRDHFCuts::Class());
379592af 158
159 fCuts->PrintAll();
160}
161
162//___________________________________________________________________________
163AliCFTaskVertexingHF& AliCFTaskVertexingHF::operator=(const AliCFTaskVertexingHF& c)
164{
165 //
166 // Assignment operator
167 //
168 if (this!=&c) {
169 AliAnalysisTaskSE::operator=(c) ;
379592af 170 fCFManager = c.fCFManager;
171 fHistEventsProcessed = c.fHistEventsProcessed;
172 fCuts = c.fCuts;
173 }
174 return *this;
175}
176
177//___________________________________________________________________________
178AliCFTaskVertexingHF::AliCFTaskVertexingHF(const AliCFTaskVertexingHF& c) :
f2703bd2 179 AliAnalysisTaskSE(c),
f2703bd2 180 fCFManager(c.fCFManager),
181 fHistEventsProcessed(c.fHistEventsProcessed),
182 fCorrelation(c.fCorrelation),
183 fCountMC(c.fCountMC),
184 fCountAcc(c.fCountAcc),
185 fCountVertex(c.fCountVertex),
186 fCountRefit(c.fCountRefit),
187 fCountReco(c.fCountReco),
188 fCountRecoAcc(c.fCountRecoAcc),
189 fCountRecoITSClusters(c.fCountRecoITSClusters),
190 fCountRecoPPR(c.fCountRecoPPR),
191 fCountRecoPID(c.fCountRecoPID),
192 fEvents(c.fEvents),
193 fDecayChannel(c.fDecayChannel),
194 fFillFromGenerated(c.fFillFromGenerated),
195 fOriginDselection(c.fOriginDselection),
196 fAcceptanceUnf(c.fAcceptanceUnf),
197 fCuts(c.fCuts),
198 fUseWeight(c.fUseWeight),
199 fWeight(c.fWeight),
31da6b05 200 fNvar(c.fNvar),
201 fPartName(c.fPartName),
3ee5eb83 202 fDauNames(c.fDauNames),
b7af2639 203 fSign(c.fSign),
fbec9fa9 204 fCentralitySelection(c.fCentralitySelection),
0ada222f 205 fFakeSelection(c.fFakeSelection),
206 fRejectIfNoQuark(c.fRejectIfNoQuark),
6c62cb59 207 fUseMCVertex(c.fUseMCVertex),
ec5288c3 208 fDsOption(c.fDsOption),
209 fConfiguration(c.fConfiguration)
379592af 210{
211 //
212 // Copy Constructor
213 //
214}
215
216//___________________________________________________________________________
f2703bd2 217AliCFTaskVertexingHF::~AliCFTaskVertexingHF()
218{
379592af 219 //
220 //destructor
221 //
222 if (fCFManager) delete fCFManager ;
223 if (fHistEventsProcessed) delete fHistEventsProcessed ;
f2703bd2 224 if (fCorrelation) delete fCorrelation ;
379592af 225 if (fCuts) delete fCuts;
226}
227
f2703bd2 228//_________________________________________________________________________-
229void AliCFTaskVertexingHF::Init()
230{
231 //
232 // Initialization
233 //
234
235 if (fDebug>1) printf("AliCFTaskVertexingHF::Init()");
236 AliRDHFCuts *copyfCuts = 0x0;
8a983038 237 if (!fCuts){
238 AliFatal("No cuts defined - Exiting...");
239 return;
240 }
241
f2703bd2 242 switch (fDecayChannel){
243 case 2:{
314ec0fe 244 copyfCuts = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
ec5288c3 245 switch (fConfiguration) {
246 case kSnail: // slow configuration: all variables in
247 fNvar = 16;
248 break;
249 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
250 fNvar = 8;
251 break;
252 }
31da6b05 253 fPartName="D0";
254 fDauNames="K+pi";
f2703bd2 255 break;
256 }
257 case 21:{
314ec0fe 258 copyfCuts = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
ec5288c3 259 switch (fConfiguration) {
260 case kSnail: // slow configuration: all variables in
261 fNvar = 16;
262 break;
263 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
264 fNvar = 8;
265 break;
266 }
31da6b05 267 fPartName="Dstar";
268 fDauNames="K+pi+pi";
f2703bd2 269 break;
270 }
271 case 31:{
314ec0fe 272 copyfCuts = new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fCuts)));
ec5288c3 273 switch (fConfiguration) {
274 case kSnail: // slow configuration: all variables in
275 fNvar = 14;
276 break;
277 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
278 fNvar = 8;
279 break;
280 }
31da6b05 281 fPartName="Dplus";
282 fDauNames="K+pi+pi";
f2703bd2 283 break;
284 }
285 case 32:{
314ec0fe 286 copyfCuts = new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fCuts)));
ec5288c3 287 switch (fConfiguration) {
288 case kSnail: // slow configuration: all variables in
289 fNvar = 18;
290 break;
291 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
292 fNvar = 8;
293 break;
294 }
31da6b05 295 fPartName="Lambdac";
296 fDauNames="p+K+pi";
f2703bd2 297 break;
298 }
299 case 33:{
314ec0fe 300 copyfCuts = new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fCuts)));
ec5288c3 301 switch (fConfiguration) {
302 case kSnail: // slow configuration: all variables in
303 fNvar = 14;
304 break;
305 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
306 fNvar = 8;
307 break;
308 }
31da6b05 309 fPartName="Ds";
310 fDauNames="K+K+pi";
f2703bd2 311 break;
312 }
313 case 4:{
314ec0fe 314 copyfCuts = new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fCuts)));
ec5288c3 315 switch (fConfiguration) {
316 case kSnail: // slow configuration: all variables in
317 fNvar = 16;
318 break;
319 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
320 fNvar = 8;
321 break;
322 }
31da6b05 323 fPartName="D0";
324 fDauNames="K+pi+pi+pi";
f2703bd2 325 break;
326 }
327 default:
328 AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exiting...");
329 break;
330 }
331
332 const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
8a983038 333 if (copyfCuts){
334 copyfCuts->SetName(nameoutput);
335
336 //Post the data
337 PostData(4, copyfCuts);
338 }
339 else{
340 AliFatal("Failing initializing AliRDHFCuts object - Exiting...");
341 }
f2703bd2 342
343 return;
344}
345
379592af 346//_________________________________________________
347void AliCFTaskVertexingHF::UserExec(Option_t *)
348{
379592af 349 //
350 // Main loop function
351 //
352
353 PostData(1,fHistEventsProcessed) ;
354 PostData(2,fCFManager->GetParticleContainer()) ;
0e355b7f 355 PostData(3,fCorrelation) ;
b7af2639 356
379592af 357 if (fFillFromGenerated){
358 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
359 }
360
361 if (!fInputEvent) {
362 Error("UserExec","NO EVENT FOUND!");
363 return;
364 }
365
366 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
367
f2703bd2 368 TClonesArray *arrayBranch=0;
379592af 369
370 if(!aodEvent && AODEvent() && IsStandardAOD()) {
f2703bd2 371 // In case there is an AOD handler writing a standard AOD, use the AOD
372 // event in memory rather than the input (ESD) event.
373 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
374 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
375 // have to taken from the AOD event hold by the AliAODExtension
376 AliAODHandler* aodHandler = (AliAODHandler*)
377 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
378 if(aodHandler->GetExtensions()) {
379 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
380 AliAODEvent *aodFromExt = ext->GetAOD();
381
382 switch (fDecayChannel){
383 case 2:{
384 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
385 break;
386 }
387 case 21:{
388 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
389 break;
390 }
391 case 31:
392 case 32:
393 case 33:{
394 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
395 break;
396 }
397 case 4:{
398 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
399 break;
400 }
401 default:
402 break;
403 }
404 }
405 }
406 else {
407 switch (fDecayChannel){
408 case 2:{
409 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
410 break;
411 }
412 case 21:{
413 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
414 break;
415 }
416 case 31:
417 case 32:
418 case 33:{
419 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
420 break;
421 }
422 case 4:{
423 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm4Prong");
424 break;
425 }
426 default:
427 break;
428 }
379592af 429 }
ec5288c3 430
379592af 431 AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
432 if (!aodVtx) return;
f2703bd2 433
434 if (!arrayBranch) {
435 AliError("Could not find array of HF vertices");
436 return;
379592af 437 }
438
439 fEvents++;
6e2e4f50 440
379592af 441 fCFManager->SetRecEventInfo(aodEvent);
442 fCFManager->SetMCEventInfo(aodEvent);
443
379592af 444 //******** DEFINE number of variables of the container***** for now set at 13, in the future in the config macro.
379592af 445
379592af 446 //loop on the MC event
447
448 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
449 if (!mcArray) {
450 AliError("Could not find Monte-Carlo in AOD");
451 return;
452 }
453 Int_t icountMC = 0;
454 Int_t icountAcc = 0;
455 Int_t icountReco = 0;
456 Int_t icountVertex = 0;
457 Int_t icountRefit = 0;
458 Int_t icountRecoAcc = 0;
459 Int_t icountRecoITSClusters = 0;
460 Int_t icountRecoPPR = 0;
461 Int_t icountRecoPID = 0;
379592af 462 Int_t cquarks = 0;
f2703bd2 463
379592af 464 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
465 if (!mcHeader) {
466 AliError("Could not find MC Header in AOD");
467 return;
468 }
e11ae259 469
470 Double_t* containerInput = new Double_t[fNvar];
471 Double_t* containerInputMC = new Double_t[fNvar];
472
f2703bd2 473
474 AliCFVertexingHF* cfVtxHF=0x0;
475 switch (fDecayChannel){
476 case 2:{
367e9aa3 477 cfVtxHF = new AliCFVertexingHF2Prong(mcArray, fOriginDselection);
478 break;
f2703bd2 479 }
480 case 21:{
367e9aa3 481 cfVtxHF = new AliCFVertexingHFCascade(mcArray, fOriginDselection);
482 break;
f2703bd2 483 }
484 case 31:
485 case 32:
486 case 33:{
4e600f58 487 cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel);
f2703bd2 488 break;
489 }
490 case 4:{
491 //cfVtxHF = new AliCFVertexingHF4Prong(mcArray, originDselection); // not there yet
492 break;
493 }
494 default:
495 break;
496 }
8a983038 497 if (!cfVtxHF){
498 AliError("No AliCFVertexingHF initialized");
d49092e3 499 delete[] containerInput;
500 delete[] containerInputMC;
8a983038 501 return;
502 }
f2703bd2 503
379592af 504 Double_t zPrimVertex = aodVtx ->GetZ();
505 Double_t zMCVertex = mcHeader->GetVtxZ();
751b8274 506 if (TMath::Abs(zMCVertex) > fCuts->GetMaxVtxZ()){
507 AliDebug(3,Form("z coordinate of MC vertex = %f, it was required to be within [-%f, +%f], skipping event", zMCVertex, fCuts->GetMaxVtxZ(), fCuts->GetMaxVtxZ()));
3655ecf2 508 delete[] containerInput;
509 delete[] containerInputMC;
751b8274 510 return;
511 }
512
2bf2e62b 513 AliESDtrackCuts** trackCuts = new AliESDtrackCuts*[cfVtxHF->GetNProngs()];
514 if (fDecayChannel == 21){
515 // for the D*, setting the third element of the array of the track cuts to those for the soft pion
516 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs()-1; iProng++){
517 trackCuts[iProng]=fCuts->GetTrackCuts();
518 }
367e9aa3 519 trackCuts[2] = fCuts->GetTrackCutsSoftPi();
2bf2e62b 520 }
521 else {
522 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs(); iProng++){
523 trackCuts[iProng]=fCuts->GetTrackCuts();
524 }
525 }
526
379592af 527 //General settings: vertex, feed down and fill reco container with generated values.
528 cfVtxHF->SetRecoPrimVertex(zPrimVertex);
529 cfVtxHF->SetMCPrimaryVertex(zMCVertex);
530 cfVtxHF->SetFillFromGenerated(fFillFromGenerated);
b7af2639 531 cfVtxHF->SetNVar(fNvar);
fbec9fa9 532 cfVtxHF->SetFakeSelection(fFakeSelection);
0ada222f 533 cfVtxHF->SetRejectCandidateIfNotFromQuark(fRejectIfNoQuark);
ec5288c3 534 cfVtxHF->SetConfiguration(fConfiguration);
b7af2639 535
af548a24 536 // switch-off the trigger class selection (doesn't work for MC)
537 fCuts->SetTriggerClass("");
538
0ada222f 539 // MC vertex, to be used, in case, for pp
540 if (fUseMCVertex) fCuts->SetUseMCVertex();
af548a24 541
72156e2e 542 if (fCentralitySelection){ // keep only the requested centrality
661d6c29 543 if(fCuts->IsEventSelectedInCentrality(aodEvent)!=0) {
544 delete[] containerInput;
545 delete[] containerInputMC;
546 delete [] trackCuts;
547 return;
548 }
72156e2e 549 } else { // keep all centralities
550 fCuts->SetMinCentrality(0.);
551 fCuts->SetMaxCentrality(100.);
8b1b055d 552 }
553
b7af2639 554 Float_t centValue = fCuts->GetCentrality(aodEvent);
555 cfVtxHF->SetCentralityValue(centValue);
f2703bd2 556
ec5288c3 557 // number of tracklets - multiplicity estimator
558 Double_t multiplicity = (Double_t)(aodEvent->GetTracklets()->GetNumberOfTracklets()); // casted to double because the CF is filled with doubles
559 cfVtxHF->SetMultiplicity(multiplicity);
560
379592af 561 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
b7af2639 562 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
563 if (!mcPart){
564 AliError("Failed casting particle from MC array!, Skipping particle");
565 continue;
566 }
567 // check the MC-level cuts, must be the desidered particle
568 if (!fCFManager->CheckParticleCuts(0, mcPart)) {
569 continue; // 0 stands for MC level
570 }
571 cfVtxHF->SetMCCandidateParam(iPart);
572
573 //counting c quarks
574 cquarks += cfVtxHF->MCcquarkCounting(mcPart);
575
576 if (!(cfVtxHF->SetLabelArray())){
577 AliDebug(2,Form("Impossible to set the label array (decaychannel = %d)",fDecayChannel));
578 continue;
579 }
2bf2e62b 580
b7af2639 581 //check the candiate family at MC level
582 if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) {
583 AliDebug(2,Form("Check on the family wrong!!! (decaychannel = %d)",fDecayChannel));
584 continue;
585 }
f2703bd2 586 else{
6e2e4f50 587 AliDebug(2,Form("Check on the family OK!!! (decaychannel = %d)",fDecayChannel));
f2703bd2 588 }
589
590 //Fill the MC container
591 Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
f2703bd2 592 if (mcContainerFilled) {
593 if (fUseWeight)fWeight = GetWeight(containerInputMC[0]);
594 if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) continue;
595 //MC Limited Acceptance
596 if (TMath::Abs(containerInputMC[1]) < 0.5) {
597 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc, fWeight);
598 AliDebug(3,"MC Lim Acc container filled\n");
599 }
600
601 //MC
602 fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight);
603 icountMC++;
604 AliDebug(3,"MC cointainer filled \n");
605
606 // MC in acceptance
607 // check the MC-Acceptance level cuts
608 // since standard CF functions are not applicable, using Kine Cuts on daughters
609 Bool_t mcAccepStep = cfVtxHF-> MCAcceptanceStep();
610 if (mcAccepStep){
611 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance, fWeight);
612 AliDebug(3,"MC acceptance cut passed\n");
613 icountAcc++;
614
615 //MC Vertex step
616 if (fCuts->IsEventSelected(aodEvent)){
617 // filling the container if the vertex is ok
618 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex, fWeight) ;
619 AliDebug(3,"Vertex cut passed and container filled\n");
620 icountVertex++;
621
622 //mc Refit requirement
2bf2e62b 623 Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, &trackCuts[0]);
f2703bd2 624 if (mcRefitStep){
625 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight);
626 AliDebug(3,"MC Refit cut passed and container filled\n");
627 icountRefit++;
628 }
629 else{
630 AliDebug(3,"MC Refit cut not passed\n");
631 continue;
632 }
633 }
634 else{
379592af 635 AliDebug (3, "MC vertex step not passed\n");
636 continue;
f2703bd2 637 }
638 }
639 else{
640 AliDebug (3, "MC in acceptance step not passed\n");
641 continue;
642 }
643 }
644 else {
645 AliDebug (3, "MC container not filled\n");
379592af 646 }
647 }
f2703bd2 648
649 if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks));
31da6b05 650 AliDebug(2,Form("Found %i MC particles that are %s!!",icountMC,fPartName.Data()));
651 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Acc cuts!!",icountAcc,fPartName.Data()));
652 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Vertex cuts!!",icountVertex,fPartName.Data()));
653 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Refit cuts!!",icountRefit,fPartName.Data()));
f2703bd2 654
379592af 655 // Now go to rec level
656 fCountMC += icountMC;
657 fCountAcc += icountAcc;
658 fCountVertex+= icountVertex;
659 fCountRefit+= icountRefit;
f2703bd2 660
6e2e4f50 661 AliDebug(2,Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel));
379592af 662
f2703bd2 663 for(Int_t iCandid = 0; iCandid<arrayBranch->GetEntriesFast();iCandid++){
f2703bd2 664 AliAODRecoDecayHF* charmCandidate=0x0;
665 switch (fDecayChannel){
666 case 2:{
667 charmCandidate = (AliAODRecoDecayHF2Prong*)arrayBranch->At(iCandid);
668 break;
669 }
670 case 21:{
671 charmCandidate = (AliAODRecoCascadeHF*)arrayBranch->At(iCandid);
672 break;
673 }
674 case 31:
675 case 32:
676 case 33:{
677 charmCandidate = (AliAODRecoDecayHF3Prong*)arrayBranch->At(iCandid);
678 break;
679 }
680 case 4:{
681 charmCandidate = (AliAODRecoDecayHF4Prong*)arrayBranch->At(iCandid);
682 break;
683 }
684 default:
685 break;
686 }
687
379592af 688 Bool_t unsetvtx=kFALSE;
f2703bd2 689 if(!charmCandidate->GetOwnPrimaryVtx()) {
690 charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables
379592af 691 unsetvtx=kTRUE;
692 }
f2703bd2 693
694 Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate);
379592af 695 if (!signAssociation){
f2703bd2 696 charmCandidate = 0x0;
379592af 697 continue;
698 }
3ee5eb83 699
700 Int_t isPartOrAntipart = cfVtxHF->CheckReflexion(fSign);
701 if (isPartOrAntipart == 0){
702 AliDebug(2, Form("The candidate pdg code doesn't match the requirement set in the task (fSign = %d)",fSign));
703 continue;
704 }
705
6694a2a8 706
379592af 707 Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
708 if (recoContFilled){
6694a2a8 709
710 // weight according to pt
711 if (fUseWeight)fWeight = GetWeight(containerInput[0]);
712
379592af 713 if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue;
f2703bd2 714
379592af 715 //Reco Step
716 Bool_t recoStep = cfVtxHF->RecoStep();
717 Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
f2703bd2 718
379592af 719 if (recoStep && recoContFilled && vtxCheck){
f2703bd2 720 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
379592af 721 icountReco++;
f2703bd2 722 AliDebug(3,"Reco step passed and container filled\n");
723
379592af 724 //Reco in the acceptance -- take care of UNFOLDING!!!!
2bf2e62b 725 Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(&trackCuts[0]);
379592af 726 if (recoAcceptanceStep) {
f2703bd2 727 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
379592af 728 icountRecoAcc++;
f2703bd2 729 AliDebug(3,"Reco acceptance cut passed and container filled\n");
730
379592af 731 if(fAcceptanceUnf){
732 Double_t fill[4]; //fill response matrix
733 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
734 if (bUnfolding) fCorrelation->Fill(fill);
f2703bd2 735 }
736
379592af 737 //Number of ITS cluster requirements
f2703bd2 738 Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
379592af 739 if (recoITSnCluster){
f2703bd2 740 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
379592af 741 icountRecoITSClusters++;
f2703bd2 742 AliDebug(3,"Reco n ITS cluster cut passed and container filled\n");
743
744 Bool_t iscutsusingpid = fCuts->GetIsUsePID();
745 Int_t recoAnalysisCuts = -1, recoPidSelection = -1;
746 fCuts->SetUsePID(kFALSE);
747 recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
748
6c62cb59 749 if (fDecayChannel==33){ // Ds case, where more possibilities are considered
750 Bool_t keepDs=ProcessDs(recoAnalysisCuts);
751 if(keepDs) recoAnalysisCuts=3;
f2703bd2 752 }
6c62cb59 753
379592af 754
f2703bd2 755 fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object
756 if (recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart){
757 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight);
758 icountRecoPPR++;
759 AliDebug(3,"Reco Analysis cuts passed and container filled \n");
760 //pid selection
761 //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID);
762 //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){
763 recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
6c62cb59 764
765 if (fDecayChannel==33){ // Ds case, where more possibilities are considered
766 Bool_t keepDs=ProcessDs(recoPidSelection);
767 if(keepDs) recoPidSelection=3;
768 }
769
f2703bd2 770 if (recoPidSelection == 3 || recoPidSelection == isPartOrAntipart){
771 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight);
772 icountRecoPID++;
773 AliDebug(3,"Reco PID cuts passed and container filled \n");
774 if(!fAcceptanceUnf){
775 Double_t fill[4]; //fill response matrix
776 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
777 if (bUnfolding) fCorrelation->Fill(fill);
778 }
779 }
780 else {
781 AliDebug(3, "Analysis Cuts step not passed \n");
782 continue;
783 }
379592af 784 }
785 else {
f2703bd2 786 AliDebug(3, "PID selection not passed \n");
787 continue;
379592af 788 }
789 }
790 else{
f2703bd2 791 AliDebug(3, "Number of ITS cluster step not passed\n");
792 continue;
379592af 793 }
794 }
795 else{
f2703bd2 796 AliDebug(3, "Reco acceptance step not passed\n");
797 continue;
379592af 798 }
799 }
800 else {
f2703bd2 801 AliDebug(3, "Reco step not passed\n");
802 continue;
379592af 803 }
804 }
805
f2703bd2 806 if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
807 } // end loop on candidate
808
379592af 809 fCountReco+= icountReco;
810 fCountRecoAcc+= icountRecoAcc;
811 fCountRecoITSClusters+= icountRecoITSClusters;
812 fCountRecoPPR+= icountRecoPPR;
813 fCountRecoPID+= icountRecoPID;
814
815 fHistEventsProcessed->Fill(0);
f2703bd2 816
817 delete[] containerInput;
f2703bd2 818 delete[] containerInputMC;
f2703bd2 819 delete cfVtxHF;
2bf2e62b 820 if (trackCuts){
821 // for (Int_t i=0; i<cfVtxHF->GetNProngs(); i++){
822 // delete [] trackCuts[i];
823 // }
824 delete [] trackCuts;
825 }
ec5288c3 826
827
379592af 828}
829
830//___________________________________________________________________________
831void AliCFTaskVertexingHF::Terminate(Option_t*)
832{
833 // The Terminate() function is the last function to be called during
834 // a query. It always runs on the client, it can be used to present
835 // the results graphically or save the results to file.
836
837 AliAnalysisTaskSE::Terminate();
31da6b05 838
839 AliInfo(Form("Found %i MC particles that are %s in MC, in %d events",fCountMC,fPartName.Data(),fEvents));
840 AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, in %d events",fCountAcc,fPartName.Data(),fEvents));
841 AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, and satisfy Vertex requirement in %d events",fCountVertex,fPartName.Data(),fEvents));
842 AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, and satisfy ITS+TPC refit requirementin %d events",fCountRefit,fPartName.Data(),fEvents));
843 AliInfo(Form("Found %i reco %s that are decaying in %s, in %d events",fCountReco,fPartName.Data(),fDauNames.Data(),fEvents));
844 AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and are in the requested acceptance, in %d events",fCountRecoAcc,fPartName.Data(),fDauNames.Data(),fEvents));
845 AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and have at least 5 clusters in ITS, in %d events",fCountRecoITSClusters,fPartName.Data(),fDauNames.Data(),fEvents));
846 AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and satisfy PPR cuts, in %d events",fCountRecoPPR,fPartName.Data(),fDauNames.Data(),fEvents));
847 AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and satisfy PPR+PID cuts, in %d events",fCountRecoPID,fPartName.Data(),fDauNames.Data(),fEvents));
379592af 848
849 // draw some example plots....
ec5288c3 850 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
379592af 851 if(!cont) {
852 printf("CONTAINER NOT FOUND\n");
853 return;
854 }
855 // projecting the containers to obtain histograms
856 // first argument = variable, second argument = step
ec5288c3 857
858 TH1D** h = new TH1D*[3];
859 if (fConfiguration == kSnail){
860 //h = new TH1D[3][12];
861 for (Int_t ih = 0; ih<3; ih++){
862 h[ih] = new TH1D[12];
863 }
864 for(Int_t iC=1;iC<12; iC++){
865 // MC-level
866 h[0][iC] = *(cont->ShowProjection(iC,0));
867 // MC-Acceptance level
868 h[1][iC] = *(cont->ShowProjection(iC,1));
869 // Reco-level
870 h[2][iC] = *(cont->ShowProjection(iC,4));
871 }
31da6b05 872 }
ec5288c3 873 else{
874 //h = new TH1D[3][12];
875 for (Int_t ih = 0; ih<3; ih++){
876 h[ih] = new TH1D[8];
877 }
878 for(Int_t iC=0;iC<8; iC++){
879 // MC-level
880 h[0][iC] = *(cont->ShowProjection(iC,0));
881 // MC-Acceptance level
882 h[1][iC] = *(cont->ShowProjection(iC,1));
883 // Reco-level
884 h[2][iC] = *(cont->ShowProjection(iC,4));
885 }
886 }
887
888 TString* titles;
889 Int_t nvarToPlot = 0;
890 if (fConfiguration == kSnail){
891 nvarToPlot = 12;
892 titles = new TString[nvarToPlot];
893 if(fDecayChannel==31){
894 titles[0]="pT_Dplus (GeV/c)";
895 titles[1]="rapidity";
896 titles[2]="phi (rad)";
897 titles[3]="cT (#mum)";
898 titles[4]="cosPointingAngle";
899 titles[5]="pT_1 (GeV/c)";
900 titles[6]="pT_2 (GeV/c)";
901 titles[7]="pT_3 (GeV/c)";
902 titles[8]="d0_1 (#mum)";
903 titles[9]="d0_2 (#mum)";
904 titles[10]="d0_3 (#mum)";
905 titles[11]="zVertex (cm)";
906 }else{
907 titles[0]="pT_D0 (GeV/c)";
908 titles[1]="rapidity";
909 titles[2]="cosThetaStar";
910 titles[3]="pT_pi (GeV/c)";
911 titles[4]="pT_K (Gev/c)";
912 titles[5]="cT (#mum)";
913 titles[6]="dca (#mum)";
914 titles[7]="d0_pi (#mum)";
915 titles[8]="d0_K (#mum)";
916 titles[9]="d0xd0 (#mum^2)";
917 titles[10]="cosPointingAngle";
918 titles[11]="phi (rad)";
919
920 }
921 }
922 else{
923 nvarToPlot = 8;
924 titles = new TString[nvarToPlot];
925 titles[0]="pT_candidate (GeV/c)";
926 titles[1]="rapidity";
927 titles[2]="cT (#mum)";
928 titles[3]="phi";
929 titles[4]="z_{vtx}";
930 titles[5]="centrality";
931 titles[6]="fake";
932 titles[7]="multiplicity";
31da6b05 933 }
ec5288c3 934
31da6b05 935 Int_t markers[12]={20,24,21,25,27,28,
936 20,24,21,25,27,28};
937 Int_t colors[3]={2,8,4};
ec5288c3 938 for(Int_t iC=0;iC<nvarToPlot; iC++){
939 for(Int_t iStep=0;iStep<3;iStep++){
940 h[iStep][iC].SetTitle(titles[iC].Data());
941 h[iStep][iC].GetXaxis()->SetTitle(titles[iC].Data());
942 Double_t maxh=h[iStep][iC].GetMaximum();
943 h[iStep][iC].GetYaxis()->SetRangeUser(0,maxh*1.2);
944 h[iStep][iC].SetMarkerStyle(markers[iC]);
945 h[iStep][iC].SetMarkerColor(colors[iStep]);
946 }
31da6b05 947 }
379592af 948
949 gStyle->SetCanvasColor(0);
950 gStyle->SetFrameFillColor(0);
951 gStyle->SetTitleFillColor(0);
952 gStyle->SetStatColor(0);
953
954 // drawing in 2 separate canvas for a matter of clearity
ec5288c3 955 TCanvas * c1 =new TCanvas(Form("c1New_%d",fDecayChannel),"Vars 0, 1, 2, 3",1100,1200);
956 c1->Divide(3,4);
31da6b05 957 Int_t iPad=1;
ec5288c3 958 for(Int_t iVar=0; iVar<4; iVar++){
959 c1->cd(iPad++);
960 h[0][iVar].Draw("p");
961 c1->cd(iPad++);
962 h[1][iVar].Draw("p");
963 c1->cd(iPad++);
964 h[2][iVar].Draw("p");
31da6b05 965 }
379592af 966
ec5288c3 967 TCanvas * c2 =new TCanvas(Form("c2New_%d",fDecayChannel),"Vars 4, 5, 6, 7",1100,1200);
968 c2->Divide(3,4);
31da6b05 969 iPad=1;
ec5288c3 970 for(Int_t iVar=4; iVar<8; iVar++){
971 c2->cd(iPad++);
972 h[0][iVar].Draw("p");
973 c2->cd(iPad++);
974 h[1][iVar].Draw("p");
975 c2->cd(iPad++);
976 h[2][iVar].Draw("p");
31da6b05 977 }
978
ec5288c3 979 if (fConfiguration == kSnail){
980 TCanvas * c3 =new TCanvas(Form("c3New_%d",fDecayChannel),"Vars 8, 9, 10, 11",1100,1200);
981 c3->Divide(3,4);
982 iPad=1;
983 for(Int_t iVar=8; iVar<12; iVar++){
984 c3->cd(iPad++);
985 h[0][iVar].Draw("p");
986 c3->cd(iPad++);
987 h[1][iVar].Draw("p");
988 c3->cd(iPad++);
989 h[2][iVar].Draw("p");
990 }
31da6b05 991 }
992
379592af 993
994 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
995
996 TH2D* corr1 =hcorr->Projection(0,2);
997 TH2D* corr2 = hcorr->Projection(1,3);
998
ec5288c3 999 TCanvas * c7 =new TCanvas(Form("c7New_%d",fDecayChannel),"",800,400);
379592af 1000 c7->Divide(2,1);
1001 c7->cd(1);
1002 corr1->Draw("text");
1003 c7->cd(2);
1004 corr2->Draw("text");
1005
379592af 1006 TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
1007
1008 corr1->Write();
1009 corr2->Write();
ec5288c3 1010 for(Int_t iC=0;iC<nvarToPlot; iC++){
1011 for(Int_t iStep=0;iStep<3;iStep++){
1012 h[iStep][iC].Write(Form("Step%d_%s",iStep,titles[iC].Data()));
1013 }
31da6b05 1014 }
379592af 1015 file_projection->Close();
6299d56c 1016
1017 for (Int_t ih = 0; ih<3; ih++) delete h[ih];
dc818294 1018 delete [] h;
6299d56c 1019
379592af 1020
379592af 1021}
1022
1023//___________________________________________________________________________
f2703bd2 1024void AliCFTaskVertexingHF::UserCreateOutputObjects()
1025{
379592af 1026 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1027 //TO BE SET BEFORE THE EXECUTION OF THE TASK
1028 //
1029 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1030
1031 //slot #1
1032 OpenFile(1);
1033 fHistEventsProcessed = new TH1I("CFHFchist0","",1,0,1) ;
0e355b7f 1034
1035 PostData(1,fHistEventsProcessed) ;
1036 PostData(2,fCFManager->GetParticleContainer()) ;
1037 PostData(3,fCorrelation) ;
1038
379592af 1039}
1040
f2703bd2 1041
1042//_________________________________________________________________________
1043Double_t AliCFTaskVertexingHF::GetWeight(Float_t pt)
1044{
1045 //
1046 // calculating the weight to fill the container
1047 //
1048
1049 // FNOLL central:
1050 // p0 = 1.63297e-01 --> 0.322643
1051 // p1 = 2.96275e+00
1052 // p2 = 2.30301e+00
1053 // p3 = 2.50000e+00
1054
1055 // PYTHIA
1056 // p0 = 1.85906e-01 --> 0.36609
1057 // p1 = 1.94635e+00
1058 // p2 = 1.40463e+00
1059 // p3 = 2.50000e+00
1060
1061 Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
1062 Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
1063
1064 Double_t dndpt_func1 = dNdptFit(pt,func1);
1065 Double_t dndpt_func2 = dNdptFit(pt,func2);
1066 AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
1067 return dndpt_func1/dndpt_func2;
1068}
1069
1070//__________________________________________________________________________________________________
1071Double_t AliCFTaskVertexingHF::dNdptFit(Float_t pt, Double_t* par)
1072{
1073 //
1074 // calculating dNdpt
1075 //
1076
1077 Double_t denom = TMath::Power((pt/par[1]), par[3] );
1078 Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
1079
1080 return dNdpt;
1081}
6c62cb59 1082
1083
1084//__________________________________________________________________________________________________
1085Bool_t AliCFTaskVertexingHF::ProcessDs(Int_t recoAnalysisCuts) const{
1086 // processes output of Ds is selected
1087 Bool_t keep=kFALSE;
1088 if(recoAnalysisCuts > 0){
1089 Int_t isKKpi=recoAnalysisCuts&1;
1090 Int_t ispiKK=recoAnalysisCuts&2;
1091 Int_t isPhiKKpi=recoAnalysisCuts&4;
1092 Int_t isPhipiKK=recoAnalysisCuts&8;
1093 Int_t isK0starKKpi=recoAnalysisCuts&16;
1094 Int_t isK0starpiKK=recoAnalysisCuts&32;
1095 if(fDsOption==1){
1096 if(isKKpi && isPhiKKpi) keep=kTRUE;
1097 if(ispiKK && isPhipiKK) keep=kTRUE;
1098 }
1099 else if(fDsOption==2){
1100 if(isKKpi && isK0starKKpi) keep=kTRUE;
1101 if(ispiKK && isK0starpiKK) keep=kTRUE;
1102 }
1103 else if(fDsOption==3)keep=kTRUE;
1104 }
1105 return keep;
1106}