1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //----------------------------------------------------------------------------
17 // Implementation of the BtoJPSItoEle reconstruction and analysis class
18 // Note: the two decay tracks are labelled: 0 (positive track)
20 // An example of usage can be found in the macro AliBtoJPSItoEleTest.C
21 // Origin: G.E. Bruno giuseppe.bruno@ba.infn.it
22 // based on Class for charm golden channel (D0->Kpi)
23 //----------------------------------------------------------------------------
29 #include <TParticle.h>
30 #include "AliESDEvent.h"
33 #include "AliRunLoader.h"
34 #include "AliVertexerTracks.h"
35 #include "AliESDVertex.h"
37 #include "AliBtoJPSItoEle.h"
38 #include "AliBtoJPSItoEleAnalysis.h"
40 #include "AliKFParticleBase.h"
41 #include "AliKFParticle.h"
42 #include "AliKFVertex.h"
56 ClassImp(AliBtoJPSItoEleAnalysis)
58 //----------------------------------------------------------------------------
59 AliBtoJPSItoEleAnalysis::AliBtoJPSItoEleAnalysis():
60 fVertexOnTheFly(kFALSE),
63 fOnlyPrimaryJpsi(kFALSE),
69 //fKFPrimVertex(kFALSE),
70 //fKFTopConstr(kFALSE),
71 fKFSecondVertex(kFALSE)
73 // Default constructor
78 //----------------------------------------------------------------------------
79 AliBtoJPSItoEleAnalysis::~AliBtoJPSItoEleAnalysis() {}
80 //----------------------------------------------------------------------------
81 void AliBtoJPSItoEleAnalysis::ApplySelection(const Char_t *inName,const Char_t *outName) const {
82 // select candidates that pass fBCuts and write them to a new file
84 TFile *inFile = TFile::Open(inName);
86 TTree *treeBin=(TTree*)inFile->Get("TreeB");
87 AliBtoJPSItoEleAnalysis *inAnalysis = (AliBtoJPSItoEleAnalysis*)inFile->Get("BtoJPSItoEleAnalysis");
88 printf("+++\n+++ I N P U T S T A T U S:\n+++\n");
89 inAnalysis->PrintStatus();
92 AliBtoJPSItoEle *d = 0;
93 treeBin->SetBranchAddress("BtoJPSItoEle",&d);
94 Int_t entries = (Int_t)treeBin->GetEntries();
96 printf("+++\n+++ Number of B candidates in input tree: %d\n+++\n",entries);
98 TTree *treeBout = new TTree("TreeB","Tree with selected B candidates");
99 treeBout->Branch("BtoJPSItoEle","AliBtoJPSItoEle",&d,200000,0);
105 for(Int_t i=0; i<entries; i++) {
106 // get event from tree
107 treeBin->GetEvent(i);
109 if(fSim && fOnlySignal && !d->IsSignal()) continue;
111 // check if candidate passes selection (as B or Bbar)
112 if(d->Select(fBCuts,okB)) {
119 AliBtoJPSItoEleAnalysis *outAnalysis = (AliBtoJPSItoEleAnalysis*)inAnalysis->Clone("BtoJPSItoEleAnalysis");
120 outAnalysis->SetBCuts(fBCuts);
121 printf("------------------------------------------\n");
122 printf("+++\n+++ O U T P U T S T A T U S:\n+++\n");
123 outAnalysis->PrintStatus();
125 printf("+++\n+++ Number of B mesons in output tree: %d\n+++\n",nSel);
127 TFile* outFile = new TFile(outName,"recreate");
129 outAnalysis->Write();
134 //----------------------------------------------------------------------------
135 Double_t AliBtoJPSItoEleAnalysis::CalculateTOFmass(Double_t mom,Double_t length,
136 Double_t time) const {
137 // calculated the mass from momentum, track length from vertex to TOF
138 // and time measured by the TOF
139 if(length==0.) return -1000.;
140 Double_t a = time*time/length/length;
142 a = TMath::Sqrt(a-1.);
144 a = -TMath::Sqrt(1.-a);
149 //----------------------------------------------------------------------------
150 void AliBtoJPSItoEleAnalysis::FindCandidates(Int_t evFirst,Int_t evLast,
151 const Char_t *outName) {
152 // Find candidates and calculate parameters
155 TString esdName="AliESDs.root";
156 if(gSystem->AccessPathName(esdName.Data(),kFileExists)) {
157 printf("AliBtoJPSItoEleAnalysis::FindCandidatesESD(): No ESDs file found!\n");
161 TString outName1=outName;
162 TString outName2="nTotEvents.dat";
164 Int_t nTotEv=0,nBrec=0,nBrec1ev=0;
166 Double_t v2[3],mom[6],d0[2];
167 Int_t iTrkP,iTrkN,trkEntries;
168 Int_t nTrksP=0,nTrksN=0;
172 AliESDtrack *postrack = 0;
173 AliESDtrack *negtrack = 0;
175 // create the AliVertexerTracks object
176 // (it will be used only if fVertexOnTheFly=kTrue)
177 AliVertexerTracks *vertexer1 = new AliVertexerTracks();
178 if(fVertexOnTheFly) {
179 // open the mean vertex
180 TFile *invtx = new TFile("AliESDVertexMean.root");
181 AliESDVertex *initVertex = (AliESDVertex*)invtx->Get("vtxmean");
183 vertexer1->SetVtxStart(initVertex);
189 // create tree for reconstructed decayes
190 AliBtoJPSItoEle *ioBtoJPSItoEle=0;
191 TTree *treeB = new TTree("TreeB","Tree with candidates");
192 treeB->Branch("BtoJPSItoEle","AliBtoJPSItoEle",&ioBtoJPSItoEle,200000,0);
194 // open file with tracks
195 TFile *esdFile = TFile::Open(esdName.Data());
196 AliESDEvent* event = new AliESDEvent();
197 TTree* tree = (TTree*) esdFile->Get("esdTree");
199 Error("FindCandidatesESD", "no ESD tree found");
202 event->ReadFromTree(tree);
204 /* if (fKFPrimVertex)
205 AliRunLoader* runLoader = 0;
208 delete gAlice->GetRunLoader();
212 runLoader = AliRunLoader::Open(galName.Data());
213 if (runLoader == 0x0) {
214 cerr<<"Can not open session"<<endl;
217 cout << "Ok open galice.root" << endl;
218 runLoader->LoadgAlice();
220 gAlice = runLoader->GetAliRun();
221 runLoader->LoadKinematics();
222 runLoader->LoadHeader();
225 // loop on events in file
226 for(Int_t iEvent = evFirst; iEvent < tree->GetEntries(); iEvent++) {
227 if(iEvent > evLast) break;
228 tree->GetEvent(iEvent);
229 Int_t ev = (Int_t)event->GetEventNumberInFile();
230 printf("--- Finding B -> JPSI -> e+ e- in event %d\n",ev);
231 // count the total number of events
234 trkEntries = (Int_t)event->GetNumberOfTracks();
235 printf(" Number of tracks: %d\n",trkEntries);
236 if(trkEntries<2) continue;
238 // retrieve primary vertex from file
239 if(!event->GetPrimaryVertex()) {
240 printf(" No vertex\n");
243 event->GetPrimaryVertex()->GetXYZ(fV1);
245 // call function which applies sigle-track selection and
246 // separetes positives and negatives
247 TObjArray trksP(trkEntries/2);
248 Int_t *trkEntryP = new Int_t[trkEntries];
249 TObjArray trksN(trkEntries/2);
250 Int_t *trkEntryN = new Int_t[trkEntries];
251 TTree *trkTree = new TTree();
252 SelectTracks(event,trksP,trkEntryP,nTrksP,
253 trksN,trkEntryN,nTrksN);
255 printf(" pos. tracks: %d neg .tracks: %d\n",nTrksP,nTrksN);
261 //===================== PRIMARY VERTEX USING KF METHODS ==========================//
265 AliStack* stack = runLoader->Stack();
271 AliKFParticle fParticle; // KFParticle constructed from ESD track
272 //Bool_t fPrimUsedFlag; // flag says that the particle was used for primary vertex fit
273 Bool_t fOK; // is the track good enough
274 Int_t mcPDG; // Monte Carlo PDG code of the particle
275 Int_t mcMotherID; // Monte Carlo ID of its mother
278 // TESDTrackInfo ESDTrackInfo[trkEntries];
279 TESDTrackInfo ESDTrackInfo[1000];
280 for (Int_t iTr=0; iTr<trkEntries; iTr++)
282 TESDTrackInfo &info = ESDTrackInfo[iTr];
284 //info.fPrimUsedFlag = 0;
286 info.mcMotherID = -1;
288 // track quality check
290 AliESDtrack *pTrack = event->GetTrack(iTr);
291 if( !pTrack ) continue;
292 if (pTrack->GetKinkIndex(0)>0) continue;
293 if ( !( pTrack->GetStatus()&AliESDtrack::kITSrefit ) ) continue;
295 //if( pTrack->GetITSclusters(indi) <5 ) continue;
296 //Int_t PDG = ( pTrack->GetSigned1Pt() <0 ) ?321 :211;
300 Int_t mcID = TMath::Abs(pTrack->GetLabel());
301 TParticle * part = stack->Particle(TMath::Abs(mcID));
302 info.mcPDG = part->GetPdgCode();
303 Int_t PDG = info.mcPDG;
304 if( mcID>=0 ) info.mcMotherID = part->GetFirstMother();
307 // Construct KFParticle for the track
309 info.fParticle = AliKFParticle( *pTrack, PDG );
313 // Find event primary vertex with KF methods
317 // const AliKFParticle * vSelected[trkEntries]; // Selected particles for vertex fit
318 // Int_t vIndex[trkEntries]; // Indices of selected particles
319 // Bool_t vFlag[trkEntries]; // Flags returned by the vertex finder
320 const AliKFParticle * vSelected[1000]; // Selected particles for vertex fit
321 Int_t vIndex[1000]; // Indices of selected particles
322 Bool_t vFlag[1000]; // Flags returned by the vertex finder
325 for( Int_t i = 0; i<trkEntries; i++){
326 if(ESDTrackInfo[i].fOK ){
327 vSelected[nSelected] = &(ESDTrackInfo[i].fParticle);
328 vIndex[nSelected] = i;
332 primVtx.ConstructPrimaryVertex( vSelected, nSelected, vFlag, 3. );
333 //for( Int_t i = 0; i<nSelected; i++){
334 // if( vFlag[i] ) ESDTrackInfo[vIndex[i]].fPrimUsedFlag = 1;
336 if( primVtx.GetNDF() <1 ) return; // Less then two tracks in primary vertex
341 // LOOP ON POSITIVE TRACKS
342 for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
343 if(iTrkP%1==0) printf(" Processing positive track number %d of %d\n",iTrkP,nTrksP);
345 // get track from track array
346 postrack = (AliESDtrack*)trksP.UncheckedAt(iTrkP);
347 trkNum[0] = trkEntryP[iTrkP];
349 // LOOP ON NEGATIVE TRACKS
350 for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
352 // get track from tracks array
353 negtrack = (AliESDtrack*)trksN.UncheckedAt(iTrkN);
354 trkNum[1] = trkEntryN[iTrkN];
357 // ----------- DCA MINIMIZATION ------------------
359 // find the DCA and propagate the tracks to the DCA
360 Double_t b=event->GetMagneticField();
361 AliESDtrack nt(*negtrack), pt(*postrack);
362 dca = nt.PropagateToDCA(&pt,b);
364 // define the AliESDv0 object
365 AliESDv0 vertex2(nt,trkNum[0],pt,trkNum[1]);
367 // get position of the secondary vertex
369 if (fKFSecondVertex){
370 //Define the AliKFParticle Objects
371 AliKFParticle trackP = AliKFParticle(pt,-11);
372 AliKFParticle trackN = AliKFParticle(nt,11);
373 //Construct the V0like mother
374 AliKFParticle V0(trackP,trackN);
375 //Get global position of the secondary vertex using KF methods
379 mom[0] = trackP.GetPx(); mom[1] = trackP.GetPy(); mom[2] = trackP.GetPz();
380 mom[3] = trackN.GetPx(); mom[4] = trackN.GetPy(); mom[5] = trackN.GetPz();
382 //Get position of the secondary vertex
383 vertex2.GetXYZ(v2[0],v2[1],v2[2]);
384 vertex2.GetPPxPyPz(mom[0],mom[1],mom[2]);
385 vertex2.GetNPxPyPz(mom[3],mom[4],mom[5]);
389 // no vertexing if DeltaMass > fMassCut
390 if(fVertexOnTheFly) {
392 if(SelectInvMass(mom)) {
393 // primary vertex from *other* tracks in the event
394 vertexer1->SetFieldkG(event->GetMagneticField());
395 skipped[0] = trkEntryP[iTrkP];
396 skipped[1] = trkEntryN[iTrkN];
397 vertexer1->SetSkipTracks(2,skipped);
398 AliESDVertex *vertex1onfly =
399 (AliESDVertex*)vertexer1->FindPrimaryVertex(event);
400 if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE;
401 vertex1onfly->GetXYZ(fV1);
402 //vertex1onfly->PrintStatus();
407 /* if (fKFSecondVertex&&fKFTopConstr&&fKFPrimVertex){
408 //====================== TOPOLOGICAL CONSTRAINT !!=====================================
409 //Primary vertex constructed from ESD using KF methods!!!
410 AliKFVertex primVtxCopy(*(event->GetPrimaryVertex()));
412 //Subtract Daughters from primary vertex
413 primVtxCopy -= trackP;
414 primVtxCopy -= trackN;
416 //Add V0 to the vertex in order to improve primary vertex resolution
419 //Set production vertex for V0
420 V0.SetProductionVertex(primVtxCopy);
422 //Recalculate primary vertex
423 fV1[0] = primVtxCopy.GetX();
424 fV1[1] = primVtxCopy.GetY();
425 fV1[2] = primVtxCopy.GetZ();
426 //=====================================================================================
429 // impact parameters of the tracks w.r.t. the primary vertex
430 d0[0] = 10000.*pt.GetD(fV1[0],fV1[1],b);
431 d0[1] = -10000.*nt.GetD(fV1[0],fV1[1],b);
433 // create the object AliBtoJPSItoEle
434 AliBtoJPSItoEle theB(ev,trkNum,fV1,v2,dca,mom,d0);
436 if(goodVtx1 && theB.Select(fBCuts,okB)) {
437 // get PID info from ESD
438 AliESDtrack *t0 = (AliESDtrack*)event->GetTrack(trkNum[0]);
440 t0->GetESDpid(esdpid0);
441 if(t0->GetStatus()&AliESDtrack::kTOFpid) {
442 tofmass[0] = CalculateTOFmass(t0->GetP(),
443 t0->GetIntegratedLength(),
448 AliESDtrack *t1 = (AliESDtrack*)event->GetTrack(trkNum[1]);
450 t1->GetESDpid(esdpid1);
451 if(t1->GetStatus()&AliESDtrack::kTOFpid) {
452 tofmass[1] = CalculateTOFmass(t1->GetP(),
453 t1->GetIntegratedLength(),
459 theB.SetPIDresponse(esdpid0,esdpid1);
460 theB.SetTOFmasses(tofmass);
463 ioBtoJPSItoEle=&theB;
471 } // loop on negative tracks
473 } // loop on positive tracks
479 printf(" Number of B candidates: %d\n",nBrec1ev);
480 } // loop on events in file
483 printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
484 printf("\n+++\n+++ Total number of B candidates: %d\n+++\n",nBrec);
490 // create a copy of this class to be written to output file
491 AliBtoJPSItoEleAnalysis *copy = (AliBtoJPSItoEleAnalysis*)this->Clone("BtoJPSItoEleAnalysis");
493 // add PDG codes to decay tracks in found candidates (in simulation mode)
494 // and store tree in the output file
496 TFile *outroot = new TFile(outName1.Data(),"recreate");
502 printf(" Now adding information from simulation (PDG codes) ...\n");
503 TTree *treeBsim = new TTree("TreeB","Tree with B candidates");
504 SimulationInfo(treeB,treeBsim);
506 TFile *outroot = new TFile(outName1.Data(),"recreate");
513 // write to a file the total number of events
514 FILE *outdat = fopen(outName2.Data(),"w");
515 fprintf(outdat,"%d\n",nTotEv);
520 //-----------------------------------------------------------------------------
521 void AliBtoJPSItoEleAnalysis::PrintStatus() const {
522 // Print parameters being used
524 printf("Preselections:\n");
525 printf(" fPtCut = %f GeV\n",fPtCut);
526 printf(" fd0Cut = %f micron\n",fd0Cut);
527 printf(" fMassCut = %f GeV\n",fMassCut);
528 printf(" fPidCut > %f \n",fPidCut);
529 if(fVertexOnTheFly) printf("Primary vertex on the fly\n");
531 printf("Simulation mode\n");
532 if(fOnlySignal && !(fOnlyPrimaryJpsi)) printf(" Only signal goes to file\n");
533 if(fOnlyPrimaryJpsi && !(fOnlySignal)) printf(" Only primary Jpsi go to file\n");
534 if(fOnlyPrimaryJpsi && fOnlySignal) printf(" Both signal and primary Jpsi go to file\n");
536 printf("Cuts on candidates:\n");
537 printf(" |M-MJPsi| [GeV] < %f\n",fBCuts[0]);
538 printf(" dca [micron] < %f\n",fBCuts[1]);
539 printf(" cosThetaStar < %f\n",fBCuts[2]);
540 printf(" pTP [GeV] > %f\n",fBCuts[3]);
541 printf(" pTN [GeV] > %f\n",fBCuts[4]);
542 printf(" |d0P| [micron] < %f\n",fBCuts[5]);
543 printf(" |d0N| [micron] < %f\n",fBCuts[6]);
544 printf(" d0d0 [micron^2] < %f\n",fBCuts[7]);
545 printf(" cosThetaPoint > %f\n",fBCuts[8]);
549 //-----------------------------------------------------------------------------
550 Bool_t AliBtoJPSItoEleAnalysis::SelectInvMass(const Double_t p[6]) const {
551 // Apply preselection in the invariant mass of the pair
553 Double_t mJPsi = 3.096916;
554 Double_t mel = 0.00510998902;
557 Double_t mom2[2],momTot2;
559 mom2[0] = p[0]*p[0] + p[1]*p[1] + p[2]*p[2];
560 mom2[1] = p[3]*p[3] + p[4]*p[4] + p[5]*p[5];
562 momTot2 = (p[0]+p[3])*(p[0]+p[3])+
563 (p[1]+p[4])*(p[1]+p[4])+
564 (p[2]+p[5])*(p[2]+p[5]);
567 energy[1] = TMath::Sqrt(mel*mel+mom2[1]);
568 energy[0] = TMath::Sqrt(mel*mel+mom2[0]);
570 Double_t minvJPsi = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
572 if(TMath::Abs(minvJPsi-mJPsi) < fMassCut) return kTRUE;
575 //-----------------------------------------------------------------------------
576 void AliBtoJPSItoEleAnalysis::SelectTracks(AliESDEvent *event,
577 TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
578 TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
579 // Create two TObjArrays with positive and negative tracks and
580 // apply single-track preselection
584 Int_t entr = event->GetNumberOfTracks();
586 // transfer ITS tracks from ESD to arrays and to a tree
587 for(Int_t i=0; i<entr; i++) {
589 AliESDtrack *esdtrack = event->GetTrack(i);
590 UInt_t status = esdtrack->GetStatus();
592 if(!(status&AliESDtrack::kITSin)) continue;
594 // single track selection
595 if(!SingleTrkCuts(*esdtrack,event->GetMagneticField())) continue;
597 if(esdtrack->GetSign()<0) { // negative track
598 trksN.AddLast(esdtrack);
599 trkEntryN[nTrksN] = i;
601 } else { // positive track
602 trksP.AddLast(esdtrack);
603 trkEntryP[nTrksP] = i;
607 } // loop on ESD tracks
611 //-----------------------------------------------------------------------------
612 void AliBtoJPSItoEleAnalysis::SetBCuts(Double_t cut0,Double_t cut1,
613 Double_t cut2,Double_t cut3,Double_t cut4,
614 Double_t cut5,Double_t cut6,
615 Double_t cut7,Double_t cut8) {
616 // Set the cuts for B selection
629 //-----------------------------------------------------------------------------
630 void AliBtoJPSItoEleAnalysis::SetBCuts(const Double_t cuts[9]) {
631 // Set the cuts for B selection
633 for(Int_t i=0; i<9; i++) fBCuts[i] = cuts[i];
637 //-----------------------------------------------------------------------------
639 AliBtoJPSItoEleAnalysis::SingleTrkCuts(const AliESDtrack& trk, Double_t b) const {
640 // Check if track passes some kinematical cuts
641 // Magnetic field "b" (kG)
643 if(TMath::Abs(1./trk.GetParameter()[4]) < fPtCut)
645 if(TMath::Abs(10000.*trk.GetD(fV1[0],fV1[1],b)) < fd0Cut)
647 //select only tracks with the "combined PID"
648 UInt_t status = trk.GetStatus();
649 if ((status&AliESDtrack::kESDpid)==0) return kTRUE;
652 if(r[0] < fPidCut) return kFALSE;
656 //----------------------------------------------------------------------------
657 void AliBtoJPSItoEleAnalysis::MakeTracksRefFile(AliRun *gAlice,
658 Int_t evFirst,Int_t evLast) const {
659 // Create a file with simulation info for the reconstructed tracks
661 TFile *outFile = TFile::Open("BTracksRefFile.root","recreate");
662 TFile *esdFile = TFile::Open("AliESDs.root");
664 AliMC *mc = gAlice->GetMCApp();
672 AliESDEvent* event = new AliESDEvent();
673 TTree* tree = (TTree*) esdFile->Get("esdTree");
674 event->ReadFromTree(tree);
675 // loop on events in file
676 for(Int_t iEvent=evFirst; iEvent<tree->GetEntries(); iEvent++) {
677 if(iEvent>evLast) break;
678 tree->GetEvent(iEvent);
679 Int_t ev = (Int_t)event->GetEventNumberInFile();
681 gAlice->GetEvent(ev);
683 Int_t nentr=(Int_t)event->GetNumberOfTracks();
685 // Tree for true track parameters
687 sprintf(ttname,"Tree_Ref_%d",ev);
688 TTree *reftree = new TTree(ttname,"Tree with true track params");
689 reftree->Branch("rectracks",&reftrk,"lab/I:pdg:mumlab:mumpdg:gmumlab:gmumpdg:mumprongs:Vx/F:Vy:Vz:Px:Py:Pz");
690 // reftree->Branch("rectracks",&reftrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
692 for(Int_t i=0; i<nentr; i++) {
693 AliESDtrack *esdtrack = (AliESDtrack*)event->GetTrack(i);
694 label = TMath::Abs(esdtrack->GetLabel());
696 part = (TParticle*)mc->Particle(label);
698 reftrk.pdg = part->GetPdgCode();
699 reftrk.mumlab = part->GetFirstMother();
700 if(part->GetFirstMother()>=0) {
701 mumpart = (TParticle*)gAlice->GetMCApp()->Particle(part->GetFirstMother());
702 reftrk.mumpdg = mumpart->GetPdgCode();
703 reftrk.mumprongs = mumpart->GetNDaughters();
704 reftrk.gmumlab = mumpart->GetFirstMother();
705 if(mumpart->GetFirstMother()>=0) {
706 gmumpart = (TParticle*)gAlice->GetMCApp()->Particle(mumpart->GetFirstMother());
707 reftrk.gmumpdg = gmumpart->GetPdgCode();
713 reftrk.gmumlab=part->GetFirstMother(); // If it hasn't any mother, then it has neither Gmother!
714 // reftrk.gmumlab=-1; // If it hasn't any mother, then it has neither Gmother!
716 reftrk.Vx = part->Vx();
717 reftrk.Vy = part->Vy();
718 reftrk.Vz = part->Vz();
719 reftrk.Px = part->Px();
720 reftrk.Py = part->Py();
721 reftrk.Pz = part->Pz();
738 //-----------------------------------------------------------------------------
739 void AliBtoJPSItoEleAnalysis::SimulationInfo(TTree *treeBin,TTree *treeBout) const {
740 // add pdg codes to candidate decay tracks (for sim)
742 TString refFileName("BTracksRefFile.root");
743 if(fSim && gSystem->AccessPathName(refFileName.Data(),kFileExists)) {
744 printf("AliBtoJPSItoEleAnalysis::SimulationInfo: no reference file found!\n");
747 TFile *refFile = TFile::Open(refFileName.Data());
749 Char_t refTreeName[100];
751 Int_t pdg[2],mumpdg[2],mumlab[2],gmumpdg[2],gmumlab[2];
754 // read-in reference tree for event 0 (the only event for Pb-Pb)
755 sprintf(refTreeName,"Tree_Ref_%d",0);
756 TTree *refTree0 = (TTree*)refFile->Get(refTreeName);
757 refTree0->SetBranchAddress("rectracks",&reftrk);
759 AliBtoJPSItoEle *theB = 0;
760 treeBin->SetBranchAddress("BtoJPSItoEle",&theB);
761 treeBout->Branch("BtoJPSItoEle","AliBtoJPSItoEle",&theB,200000,0);
763 Int_t entries = (Int_t)treeBin->GetEntries();
765 for(Int_t i=0; i<entries; i++) {
766 if(i%100==0) printf(" done %d candidates of %d\n",i,entries);
768 treeBin->GetEvent(i);
769 event = theB->EventNo();
771 if(event==0) { // always true for Pb-Pb (avoid to read-in tree every time)
772 refTree0->GetEvent(theB->GetTrkNum(0));
774 mumpdg[0] = reftrk.mumpdg;
775 mumlab[0] = reftrk.mumlab;
776 gmumpdg[0] = reftrk.gmumpdg;
777 gmumlab[0] = reftrk.gmumlab;
778 refTree0->GetEvent(theB->GetTrkNum(1));
780 mumpdg[1] = reftrk.mumpdg;
781 mumlab[1] = reftrk.mumlab;
782 gmumpdg[1] = reftrk.gmumpdg;
783 gmumlab[1] = reftrk.gmumlab;
785 sprintf(refTreeName,"Tree_Ref_%d",event);
786 TTree *refTree = (TTree*)refFile->Get(refTreeName);
787 refTree->SetBranchAddress("rectracks",&reftrk);
788 refTree->GetEvent(theB->GetTrkNum(0));
790 mumpdg[0] = reftrk.mumpdg;
791 mumlab[0] = reftrk.mumlab;
792 gmumpdg[0] = reftrk.gmumpdg;
793 gmumlab[0] = reftrk.gmumlab;
794 refTree->GetEvent(theB->GetTrkNum(1));
796 mumpdg[1] = reftrk.mumpdg;
797 mumlab[1] = reftrk.mumlab;
798 gmumpdg[1] = reftrk.gmumpdg;
799 gmumlab[1] = reftrk.gmumlab;
803 theB->SetPdgCodes(pdg);
804 theB->SetMumPdgCodes(mumpdg);
805 theB->SetGMumPdgCodes(gmumpdg);
807 if(gmumpdg[0]==gmumpdg[1] && // Both GrandMothers are of the same sign
808 (TMath::Abs(gmumpdg[0])==521 || TMath::Abs(gmumpdg[0])==511 || // GrandMother Bplus/Bminus or B0/B0bar
809 TMath::Abs(gmumpdg[0])==523 || TMath::Abs(gmumpdg[0])==513 || // B0s/B0sbar
810 TMath::Abs(gmumpdg[0])==515 || TMath::Abs(gmumpdg[0])==525 || //
811 TMath::Abs(gmumpdg[0])==531 || TMath::Abs(gmumpdg[0])==533 || //
812 TMath::Abs(gmumpdg[0])==535 || TMath::Abs(gmumpdg[0])==541 || //
813 TMath::Abs(gmumpdg[0])==543 || TMath::Abs(gmumpdg[0])==545 || //
814 TMath::Abs(gmumpdg[0])==10521 || TMath::Abs(gmumpdg[0])==10511 || // all possible
815 TMath::Abs(gmumpdg[0])==10523 || TMath::Abs(gmumpdg[0])==10513 || // B mesons
816 TMath::Abs(gmumpdg[0])==20523 || TMath::Abs(gmumpdg[0])==20513 || //
817 TMath::Abs(gmumpdg[0])==10531 || TMath::Abs(gmumpdg[0])==10533 || //
818 TMath::Abs(gmumpdg[0])==20533 || TMath::Abs(gmumpdg[0])==10541 || //
819 TMath::Abs(gmumpdg[0])==20543 || TMath::Abs(gmumpdg[0])==10543 || //
820 TMath::Abs(gmumpdg[0])==4122 || TMath::Abs(gmumpdg[0])==4222 || // All possible B baryons
821 TMath::Abs(gmumpdg[0])==4212 || TMath::Abs(gmumpdg[0])==4112 || // All possible B baryons
822 TMath::Abs(gmumpdg[0])==4224 || TMath::Abs(gmumpdg[0])==4214 || // All possible B baryons
823 TMath::Abs(gmumpdg[0])==4114 || TMath::Abs(gmumpdg[0])==4232 || // All possible B baryons
824 TMath::Abs(gmumpdg[0])==4132 || TMath::Abs(gmumpdg[0])==4322 || // All possible B baryons
825 TMath::Abs(gmumpdg[0])==4312 || TMath::Abs(gmumpdg[0])==4324 || // All possible B baryons
826 TMath::Abs(gmumpdg[0])==4314 || TMath::Abs(gmumpdg[0])==4332 || // All possible B baryons
827 TMath::Abs(gmumpdg[0])==4334 || TMath::Abs(gmumpdg[0])==4412 || // All possible B baryons
828 TMath::Abs(gmumpdg[0])==4422 || TMath::Abs(gmumpdg[0])==4414 || // All possible B baryons
829 TMath::Abs(gmumpdg[0])==4424 || TMath::Abs(gmumpdg[0])==4432 || // All possible B baryons
830 TMath::Abs(gmumpdg[0])==4434 || TMath::Abs(gmumpdg[0])==4444 // All possible B baryons
832 mumpdg[0]==443 && mumpdg[1]== 443 &&
833 mumlab[0]==mumlab[1] &&
834 reftrk.mumprongs==2 &&
835 pdg[0]==-11 && pdg[1]==11
838 else if ( // here consider the case of primary J/psi
839 mumpdg[0]==443 && mumpdg[1]== 443 &&
840 pdg[0]==-11 && pdg[1]==11 &&
841 mumlab[0]==mumlab[1] &&
842 reftrk.mumprongs==2 &&
843 ( gmumlab[0]<0 || // really primary J/psi (without family. e.g. from Cocktail)
844 TMath::Abs(gmumpdg[0])==100443 || // from Psi(2S)
845 TMath::Abs(gmumpdg[0])==10441 || // from Csi_c0(1P)
846 TMath::Abs(gmumpdg[0])==20443 || // from Csi_c1(1P)
847 TMath::Abs(gmumpdg[0])==10443 || // from h_c(1P)
848 TMath::Abs(gmumpdg[0])==445 || // from Csi_c2(1P)
849 (gmumpdg[0]>=81 && gmumpdg[0]<=100) // this is for MC internal use (e.g. J/psi from string)
851 ) theB->SetJpsiPrimary();
853 // if(!fOnlySignal || theB->IsSignal()) treeBout->Fill();
855 // write it out 1) always if you have not asked for onlySignal or OnlyPrimaryJpsi (or both)
856 // or 2) if you have asked for Signal and it is Signal
857 // or 3) if you have asked for Primary Jpsi and it is a Primary Jpsi
858 if ( (!fOnlySignal && !fOnlyPrimaryJpsi) || (fOnlySignal && theB->IsSignal())
859 || (fOnlyPrimaryJpsi && theB->IsJpsiPrimary()) ) treeBout->Fill();