void ana_dsdsj2(TString fname="dsds_10k.evt.root",int num=0)
{
    TStopwatch timer;
    timer.Start();
    
    gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
    basiclibs();
    
    // **** the library is needed in order to use the TCandidate and TCandList etc
    gSystem->Load("libRho");
    
    TCanvas *c1=new TCanvas("c1","c1",900,600);
    c1->Divide(3,2);
    
    // **** access the datastructure holding the particle lists
    //
    TFile* f = new TFile(fname.Data());
    TTree *t=f->Get("cbmsim") ;
    
    // **** for every event, a TCLonesArray with the candidates is stored in cbmsim
    //
    TClonesArray *fChrgCands=new TClonesArray("TCandidate");
    TClonesArray *fNeutCands=new TClonesArray("TCandidate");
    TClonesArray *fMcCands=new TClonesArray("TCandidate");
    
    t->SetBranchAddress("PndChargedCandidates",&fChrgCands);
    t->SetBranchAddress("PndNeutralCandidates",&fNeutCands);
    t->SetBranchAddress("PndMcTracks",&fMcCands);
    
    TCandidate *tc;
    
    // **** create and setup some histos for QA plots
    //
    TH1F *phimass = new TH1F("phimass","phi cands",100,0.95,1.1);
    TH1F *pi0mass = new TH1F("pi0mass","pi0 cands",100,0.135-0.03,0.135+0.03);
    TH1F *dsmass  = new TH1F("dsmass","Ds cands",100,1.968-0.03,1.968+0.03); 
    TH1F *ds0mass = new TH1F("ds0mass","Ds0 cands",100,2.317-0.05,2.317+0.05);
    TH1F *ppmass = new TH1F("ppmass","pbarp cands",100,4.306-0.1,4.306+0.1);
    
    TH1F *nmult=new TH1F("nmult","# neutrals",15,0,15);
    
    TH1F *pdiff=new TH1F("pdiff","momentum difference",100,-0.05,0.05);
    
    phimass->SetMinimum(0);
    pi0mass->SetMinimum(0);
    dsmass->SetMinimum(0);
    ds0mass->SetMinimum(0);
    ppmass->SetMinimum(0);
    
    if (num==0) num= t->GetEntriesFast();
    cout <<"\n####### Processing "<<num <<" events...\n"<<endl;
    
    // **** create all the particle lists we'll need for rebuilding the decay tree
    //
    TCandList mcCands;
    
    TCandList allCands,neutralCands,chargedCands, plusCands,minusCands;
    
    TCandList kpCands,kmCands,piCands;
    
    TCandList phiCands,pi0Cands,dsCands,ds0Cands,ppCands;
    
    // **** create and configure the selectors/filters we'd like to use later
    //
    TPidChargedSelector *chargedSel = new TPidChargedSelector;
    TPidNeutralSelector *neutralSel = new TPidNeutralSelector;
    TPidPlusSelector    *plusSel    = new TPidPlusSelector;
    TPidMinusSelector   *minusSel   = new TPidMinusSelector;
    
    // **** mass selectors for the resonances/composites
    //
    TPidMassSelector *phiMSel       = new TPidMassSelector("phiSelector" , 1.0195 , 0.01);
    TPidMassSelector *pi0MSel       = new TPidMassSelector("pi0Selector" , 0.135  , 0.005);
    TPidMassSelector *dsMSel        = new TPidMassSelector("dsSelector"  , 1.9685 , 0.01);
    
    TPidSimpleKaonSelector *kSel    = new TPidSimpleKaonSelector();
    kSel->SetCriterion("loose");
    TPidSimplePionSelector *piSel   = new TPidSimplePionSelector();
    piSel->SetCriterion("veryLoose");
    
    int i1,i2;
    
    // **** loop over all _events_
    //
    for (Int_t j=0; j< num;j++)
    {
       if ((j%100)==0) cout <<"evt "<<j<<endl;

	t->GetEntry(j);
	TFactory::Instance()->Reset();
	
	mcCands.Cleanup();
	chargedCands.Cleanup();
	neutralCands.Cleanup();
	
	// **** loop over all Candidates and add them to the list allCands
	//    
	
	for (i1=0; i1<fChrgCands->GetEntriesFast(); i1++){
	    tc = (TCandidate *)fChrgCands->At(i1);
	    chargedCands.Add(*tc);
	}
	
	for (i1=0; i1<fNeutCands->GetEntriesFast(); i1++){
	    tc = (TCandidate *)fNeutCands->At(i1);
	    neutralCands.Add(*tc);
	}
	
	for (i1=0; i1<fMcCands->GetEntriesFast(); i1++){
	    tc = (TCandidate *)fMcCands->At(i1);
	    mcCands.Add(*tc);
	}
	
	
	// **** select all the basic list
	//
	
	//cout <<"c:"<<chargedCands.GetLength()<<" n:"<<neutralCands.GetLength()<<endl;
	
	nmult->Fill(neutralCands.GetLength());
		
	plusCands.Select(chargedCands  ,plusSel);
	minusCands.Select(chargedCands ,minusSel);
	
	// **** pid selection
	//
	kpCands.Select(plusCands    ,kSel);
	kmCands.Select(minusCands   ,kSel);
	piCands.Select(chargedCands ,piSel);
	
	// **** now start combining all composits; inbetween plot masses
	//      before using the mass selectors
	//
	phiCands.Combine(kpCands,kmCands);
	
	TCandListIterator iterPhi(phiCands);
	while (tc=iterPhi.Next()) phimass->Fill(tc->M());
	phiCands.Select(phiMSel);
	
	dsCands.Combine(phiCands,piCands);
	
	TCandListIterator iterDs(dsCands);
	while (tc=iterDs.Next()) dsmass->Fill(tc->M()); 
	dsCands.Select(dsMSel);
	
	pi0Cands.Combine(neutralCands,neutralCands);
	
	TCandListIterator iterPi0(pi0Cands);
	while (tc=iterPi0.Next()) pi0mass->Fill(tc->M()); 
	pi0Cands.Select(pi0MSel);
	
	ds0Cands.Combine(dsCands,pi0Cands);
	
	TCandListIterator iterDs0(ds0Cands);
	while (tc=iterDs0.Next()) ds0mass->Fill(tc->M());  
	
	ppCands.Combine(ds0Cands,dsCands);
	//cout << ds0Cands.GetLength()<<" "<<dsCands.GetLength()<<" "<<ppCands.GetLength()<<endl;
	ppCands.Select(neutralSel);
	
	// **** since we didn't care about the D_s charge, combinations might appear in two
	//      different ways; RemoveClones removes double candidates based on the same final states
	ppCands.RemoveClones();
	
	TCandListIterator iterPp(ppCands);
	
	while (tc=iterPp.Next()) 
	{
	    ppmass->Fill(tc->M());
	    //printTree(tc);
	}
	
	TCandListIterator iterChrg(chargedCands);
	while(tc=iterChrg.Next()) 
	{
	    int idx=tc->GetMcIdx();
	    if (idx!=-1) pdiff->Fill(tc->P()-mcCands[idx]->P());
	}
	
	
    }
    
    // **** plot all that stuff
    //
    c1->cd(1);
    phimass->Draw();
    c1->cd(2);
    pi0mass->Draw();
    c1->cd(3);
    dsmass->Draw();
    c1->cd(4);
    ds0mass->Draw();
    c1->cd(5);
    ppmass->Draw();
    c1->cd(6);
    pdiff->Draw();
    
    timer.Stop();
    Double_t rtime = timer.RealTime();
    Double_t ctime = timer.CpuTime();
    
    printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
    
}

/*void printTree(TCandidate *tc, int level=0)
{
  int i=0;
  int nd=tc->NDaughters();
  if (nd==0) return;
  cout <<tc->Uid()<<"("<<level<<") -> ";
  for (i=0;i<nd;i++) cout<<tc->Daughter(i)->Uid()<<"  ";
  cout <<endl;
  for (i=0;i<nd;i++) printTree(tc->Daughter(i),level+1);
  if (level==0) cout <<endl;
}
*/