ana_dsdsj2.C File Reference

Go to the source code of this file.

Functions
int	ana_dsdsj2 (TString fname="dsds_10k.evt.root", int num=0)

Function Documentation

int ana_dsdsj2	(	TString	fname = "dsds_10k.evt.root",
		int	num = 0
	)

Definition at line 1 of file ana_dsdsj2.C.

References basiclibs(), c1, ctime, Double_t, f, fMcCands, idx, num, printf(), rtime, t, and timer.

{
    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("pndsim") ;

    // **** for every event, a TCLonesArray with the candidates is stored in pndsim
    //
    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);

  return 0;
}