Removing warnings (Andrea)
[u/mrichter/AliRoot.git] / PWG3 / vertexingOld / AliBtoJPSItoEleAnalysis.cxx
CommitLineData
27e190fb 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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// Implementation of the BtoJPSItoEle reconstruction and analysis class
18// Note: the two decay tracks are labelled: 0 (positive track)
19// 1 (negative 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//----------------------------------------------------------------------------
24#include <TKey.h>
25#include <TFile.h>
26#include <TTree.h>
27#include <TString.h>
28#include <TSystem.h>
29#include <TParticle.h>
5e6f195c 30#include "AliESDEvent.h"
27e190fb 31#include "AliMC.h"
32#include "AliRun.h"
33#include "AliRunLoader.h"
34#include "AliVertexerTracks.h"
35#include "AliESDVertex.h"
36#include "AliESDv0.h"
37#include "AliBtoJPSItoEle.h"
38#include "AliBtoJPSItoEleAnalysis.h"
39#include "AliLog.h"
0ac6b6e8 40#include "AliKFParticleBase.h"
41#include "AliKFParticle.h"
42#include "AliKFVertex.h"
27e190fb 43
44typedef struct {
45 Int_t lab;
46 Int_t pdg;
47 Int_t mumlab;
48 Int_t mumpdg;
49 Int_t gmumlab;
50 Int_t gmumpdg;
51 Int_t mumprongs;
52 Float_t Vx,Vy,Vz;
53 Float_t Px,Py,Pz;
54} REFTRACK;
55
56ClassImp(AliBtoJPSItoEleAnalysis)
57
58//----------------------------------------------------------------------------
dd1ac8a2 59AliBtoJPSItoEleAnalysis::AliBtoJPSItoEleAnalysis():
60fVertexOnTheFly(kFALSE),
61fSim(kFALSE),
62fOnlySignal(kFALSE),
63fOnlyPrimaryJpsi(kFALSE),
64fPID("TRDTPCparam"),
65fPtCut(0.),
66fd0Cut(0.),
67fMassCut(1000.),
0ac6b6e8 68fPidCut(0.),
69//fKFPrimVertex(kFALSE),
70//fKFTopConstr(kFALSE),
71fKFSecondVertex(kFALSE)
dd1ac8a2 72{
27e190fb 73 // Default constructor
74
27e190fb 75 SetBCuts();
76 SetVertex1();
27e190fb 77}
78//----------------------------------------------------------------------------
79AliBtoJPSItoEleAnalysis::~AliBtoJPSItoEleAnalysis() {}
80//----------------------------------------------------------------------------
81void AliBtoJPSItoEleAnalysis::ApplySelection(const Char_t *inName,const Char_t *outName) const {
82 // select candidates that pass fBCuts and write them to a new file
83
84 TFile *inFile = TFile::Open(inName);
85
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();
90
91
92 AliBtoJPSItoEle *d = 0;
93 treeBin->SetBranchAddress("BtoJPSItoEle",&d);
94 Int_t entries = (Int_t)treeBin->GetEntries();
95
96 printf("+++\n+++ Number of B candidates in input tree: %d\n+++\n",entries);
97
98 TTree *treeBout = new TTree("TreeB","Tree with selected B candidates");
99 treeBout->Branch("BtoJPSItoEle","AliBtoJPSItoEle",&d,200000,0);
100
101
102 Int_t okB=0;
103 Int_t nSel = 0;
104
105 for(Int_t i=0; i<entries; i++) {
106 // get event from tree
107 treeBin->GetEvent(i);
108
109 if(fSim && fOnlySignal && !d->IsSignal()) continue;
110
111 // check if candidate passes selection (as B or Bbar)
112 if(d->Select(fBCuts,okB)) {
113 nSel++;
114 treeBout->Fill();
115 }
116
117 }
118
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();
124
125 printf("+++\n+++ Number of B mesons in output tree: %d\n+++\n",nSel);
126
127 TFile* outFile = new TFile(outName,"recreate");
128 treeBout->Write();
129 outAnalysis->Write();
130 outFile->Close();
131
132 return;
133}
134//----------------------------------------------------------------------------
135Double_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;
141 if(a > 1.) {
142 a = TMath::Sqrt(a-1.);
143 } else {
144 a = -TMath::Sqrt(1.-a);
145 }
146
147 return mom*a;
148}
149//----------------------------------------------------------------------------
150void AliBtoJPSItoEleAnalysis::FindCandidates(Int_t evFirst,Int_t evLast,
151 const Char_t *outName) {
152 // Find candidates and calculate parameters
153
154
155 TString esdName="AliESDs.root";
156 if(gSystem->AccessPathName(esdName.Data(),kFileExists)) {
157 printf("AliBtoJPSItoEleAnalysis::FindCandidatesESD(): No ESDs file found!\n");
158 return;
159 }
160
161 TString outName1=outName;
162 TString outName2="nTotEvents.dat";
163
164 Int_t nTotEv=0,nBrec=0,nBrec1ev=0;
165 Double_t dca;
166 Double_t v2[3],mom[6],d0[2];
167 Int_t iTrkP,iTrkN,trkEntries;
168 Int_t nTrksP=0,nTrksN=0;
169 Int_t trkNum[2];
170 Double_t tofmass[2];
171 Int_t okB=0;
172 AliESDtrack *postrack = 0;
173 AliESDtrack *negtrack = 0;
174
175 // create the AliVertexerTracks object
176 // (it will be used only if fVertexOnTheFly=kTrue)
dd1ac8a2 177 AliVertexerTracks *vertexer1 = new AliVertexerTracks();
27e190fb 178 if(fVertexOnTheFly) {
179 // open the mean vertex
180 TFile *invtx = new TFile("AliESDVertexMean.root");
181 AliESDVertex *initVertex = (AliESDVertex*)invtx->Get("vtxmean");
182 invtx->Close();
183 vertexer1->SetVtxStart(initVertex);
184 delete invtx;
185 }
186 Int_t skipped[2];
187 Bool_t goodVtx1;
188
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);
193
194 // open file with tracks
195 TFile *esdFile = TFile::Open(esdName.Data());
5e6f195c 196 AliESDEvent* event = new AliESDEvent();
27e190fb 197 TTree* tree = (TTree*) esdFile->Get("esdTree");
198 if(!tree) {
199 Error("FindCandidatesESD", "no ESD tree found");
200 return;
201 }
5e6f195c 202 event->ReadFromTree(tree);
27e190fb 203
0ac6b6e8 204/* if (fKFPrimVertex)
205 AliRunLoader* runLoader = 0;
206 {
207 if (gAlice) {
208 delete gAlice->GetRunLoader();
209 delete gAlice;
210 gAlice=0;
211 }
212 runLoader = AliRunLoader::Open(galName.Data());
213 if (runLoader == 0x0) {
214 cerr<<"Can not open session"<<endl;
215 return;
216 }
217 cout << "Ok open galice.root" << endl;
218 runLoader->LoadgAlice();
219
220 gAlice = runLoader->GetAliRun();
221 runLoader->LoadKinematics();
222 runLoader->LoadHeader();
223 } */
224
27e190fb 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
232 nTotEv++;
233
234 trkEntries = (Int_t)event->GetNumberOfTracks();
235 printf(" Number of tracks: %d\n",trkEntries);
236 if(trkEntries<2) continue;
237
238 // retrieve primary vertex from file
239 if(!event->GetPrimaryVertex()) {
240 printf(" No vertex\n");
241 continue;
242 }
243 event->GetPrimaryVertex()->GetXYZ(fV1);
244
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);
254
255 printf(" pos. tracks: %d neg .tracks: %d\n",nTrksP,nTrksN);
256
257
258 nBrec1ev = 0;
259
0ac6b6e8 260/*
261//===================== PRIMARY VERTEX USING KF METHODS ==========================//
262
263if (fKFPrimVertex)
264{
265 AliStack* stack = runLoader->Stack();
266
267 class TESDTrackInfo
268 {
269 public:
270 TESDTrackInfo(){}
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
276 };
277
278 // TESDTrackInfo ESDTrackInfo[trkEntries];
279 TESDTrackInfo ESDTrackInfo[1000];
280 for (Int_t iTr=0; iTr<trkEntries; iTr++)
281 {
282 TESDTrackInfo &info = ESDTrackInfo[iTr];
283 info.fOK = 0;
284 //info.fPrimUsedFlag = 0;
285 info.mcPDG = -1;
286 info.mcMotherID = -1;
287
288 // track quality check
289
290 AliESDtrack *pTrack = event->GetTrack(iTr);
291 if( !pTrack ) continue;
292 if (pTrack->GetKinkIndex(0)>0) continue;
293 if ( !( pTrack->GetStatus()&AliESDtrack::kITSrefit ) ) continue;
294 //Int_t indi[12];
295 //if( pTrack->GetITSclusters(indi) <5 ) continue;
296 //Int_t PDG = ( pTrack->GetSigned1Pt() <0 ) ?321 :211;
297
298 // take MC PDG
299
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();
305
306
307 // Construct KFParticle for the track
308
309 info.fParticle = AliKFParticle( *pTrack, PDG );
310 info.fOK = 1;
311 }
312
313 // Find event primary vertex with KF methods
314
315 AliKFVertex primVtx;
316 {
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
323
324 Int_t nSelected = 0;
325 for( Int_t i = 0; i<trkEntries; i++){
326 if(ESDTrackInfo[i].fOK ){
327 vSelected[nSelected] = &(ESDTrackInfo[i].fParticle);
328 vIndex[nSelected] = i;
329 nSelected++;
330 }
331 }
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;
335 //}
336 if( primVtx.GetNDF() <1 ) return; // Less then two tracks in primary vertex
337 }
338
339}*/
340
27e190fb 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);
344
345 // get track from track array
346 postrack = (AliESDtrack*)trksP.UncheckedAt(iTrkP);
347 trkNum[0] = trkEntryP[iTrkP];
348
349 // LOOP ON NEGATIVE TRACKS
350 for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
351
352 // get track from tracks array
353 negtrack = (AliESDtrack*)trksN.UncheckedAt(iTrkN);
354 trkNum[1] = trkEntryN[iTrkN];
355
356 //
357 // ----------- DCA MINIMIZATION ------------------
358 //
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);
363
364 // define the AliESDv0 object
365 AliESDv0 vertex2(nt,trkNum[0],pt,trkNum[1]);
366
367 // get position of the secondary vertex
0ac6b6e8 368
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
376 v2[0] = V0.GetX();
377 v2[1] = V0.GetY();
378 v2[2] = V0.GetZ();
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();
381 }else{
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]);
386 }
27e190fb 387 goodVtx1 = kTRUE;
388
389 // no vertexing if DeltaMass > fMassCut
390 if(fVertexOnTheFly) {
391 goodVtx1 = kFALSE;
392 if(SelectInvMass(mom)) {
393 // primary vertex from *other* tracks in the event
dd1ac8a2 394 vertexer1->SetFieldkG(event->GetMagneticField());
27e190fb 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();
403 delete vertex1onfly;
404 }
405 }
406
0ac6b6e8 407/* if (fKFSecondVertex&&fKFTopConstr&&fKFPrimVertex){
408//====================== TOPOLOGICAL CONSTRAINT !!=====================================
409 //Primary vertex constructed from ESD using KF methods!!!
410 AliKFVertex primVtxCopy(*(event->GetPrimaryVertex()));
411
412 //Subtract Daughters from primary vertex
413 primVtxCopy -= trackP;
414 primVtxCopy -= trackN;
415
416 //Add V0 to the vertex in order to improve primary vertex resolution
417 primVtxCopy += V0;
418
419 //Set production vertex for V0
420 V0.SetProductionVertex(primVtxCopy);
421
422 //Recalculate primary vertex
423 fV1[0] = primVtxCopy.GetX();
424 fV1[1] = primVtxCopy.GetY();
425 fV1[2] = primVtxCopy.GetZ();
426//=====================================================================================
427 }*/
428
27e190fb 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);
432
433 // create the object AliBtoJPSItoEle
434 AliBtoJPSItoEle theB(ev,trkNum,fV1,v2,dca,mom,d0);
435 // select B's
436 if(goodVtx1 && theB.Select(fBCuts,okB)) {
437 // get PID info from ESD
438 AliESDtrack *t0 = (AliESDtrack*)event->GetTrack(trkNum[0]);
439 Double_t esdpid0[5];
440 t0->GetESDpid(esdpid0);
441 if(t0->GetStatus()&AliESDtrack::kTOFpid) {
442 tofmass[0] = CalculateTOFmass(t0->GetP(),
443 t0->GetIntegratedLength(),
444 t0->GetTOFsignal());
445 } else {
446 tofmass[0] = -1000.;
447 }
448 AliESDtrack *t1 = (AliESDtrack*)event->GetTrack(trkNum[1]);
449 Double_t esdpid1[5];
450 t1->GetESDpid(esdpid1);
451 if(t1->GetStatus()&AliESDtrack::kTOFpid) {
452 tofmass[1] = CalculateTOFmass(t1->GetP(),
453 t1->GetIntegratedLength(),
454 t1->GetTOFsignal());
455 } else {
456 tofmass[1] = -1000.;
457 }
458
459 theB.SetPIDresponse(esdpid0,esdpid1);
460 theB.SetTOFmasses(tofmass);
461
462 // fill the tree
463 ioBtoJPSItoEle=&theB;
464 treeB->Fill();
465
466 nBrec++; nBrec1ev++;
467 ioBtoJPSItoEle=0;
468 }
469
470 negtrack = 0;
471 } // loop on negative tracks
472 postrack = 0;
473 } // loop on positive tracks
474
475 delete [] trkEntryP;
476 delete [] trkEntryN;
477 delete trkTree;
478
479 printf(" Number of B candidates: %d\n",nBrec1ev);
480 } // loop on events in file
481
482
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);
485
486 delete vertexer1;
487
488 esdFile->Close();
489
490 // create a copy of this class to be written to output file
491 AliBtoJPSItoEleAnalysis *copy = (AliBtoJPSItoEleAnalysis*)this->Clone("BtoJPSItoEleAnalysis");
492
493 // add PDG codes to decay tracks in found candidates (in simulation mode)
494 // and store tree in the output file
495 if(!fSim) {
496 TFile *outroot = new TFile(outName1.Data(),"recreate");
497 treeB->Write();
498 copy->Write();
499 outroot->Close();
500 delete outroot;
501 } else {
502 printf(" Now adding information from simulation (PDG codes) ...\n");
503 TTree *treeBsim = new TTree("TreeB","Tree with B candidates");
504 SimulationInfo(treeB,treeBsim);
505 delete treeB;
506 TFile *outroot = new TFile(outName1.Data(),"recreate");
507 treeBsim->Write();
508 copy->Write();
509 outroot->Close();
510 delete outroot;
511 }
512
513 // write to a file the total number of events
514 FILE *outdat = fopen(outName2.Data(),"w");
515 fprintf(outdat,"%d\n",nTotEv);
516 fclose(outdat);
517
518 return;
519}
520//-----------------------------------------------------------------------------
521void AliBtoJPSItoEleAnalysis::PrintStatus() const {
522 // Print parameters being used
523
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");
530 if(fSim) {
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");
535 }
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]);
546
547 return;
548}
549//-----------------------------------------------------------------------------
550Bool_t AliBtoJPSItoEleAnalysis::SelectInvMass(const Double_t p[6]) const {
551 // Apply preselection in the invariant mass of the pair
552
553 Double_t mJPsi = 3.096916;
554 Double_t mel = 0.00510998902;
555
556 Double_t energy[2];
557 Double_t mom2[2],momTot2;
558
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];
561
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]);
565
566 // J/Psi -> e+ e-
567 energy[1] = TMath::Sqrt(mel*mel+mom2[1]);
568 energy[0] = TMath::Sqrt(mel*mel+mom2[0]);
569
570 Double_t minvJPsi = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
571
572 if(TMath::Abs(minvJPsi-mJPsi) < fMassCut) return kTRUE;
573 return kFALSE;
574}
575//-----------------------------------------------------------------------------
5e6f195c 576void AliBtoJPSItoEleAnalysis::SelectTracks(AliESDEvent *event,
27e190fb 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
581
582 nTrksP=0,nTrksN=0;
583
584 Int_t entr = event->GetNumberOfTracks();
585
586 // transfer ITS tracks from ESD to arrays and to a tree
587 for(Int_t i=0; i<entr; i++) {
588
589 AliESDtrack *esdtrack = event->GetTrack(i);
590 UInt_t status = esdtrack->GetStatus();
591
592 if(!(status&AliESDtrack::kITSin)) continue;
593
594 // single track selection
595 if(!SingleTrkCuts(*esdtrack,event->GetMagneticField())) continue;
596
597 if(esdtrack->GetSign()<0) { // negative track
598 trksN.AddLast(esdtrack);
599 trkEntryN[nTrksN] = i;
600 nTrksN++;
601 } else { // positive track
602 trksP.AddLast(esdtrack);
603 trkEntryP[nTrksP] = i;
604 nTrksP++;
605 }
606
607 } // loop on ESD tracks
608
609 return;
610}
611//-----------------------------------------------------------------------------
612void 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
617 fBCuts[0] = cut0;
618 fBCuts[1] = cut1;
619 fBCuts[2] = cut2;
620 fBCuts[3] = cut3;
621 fBCuts[4] = cut4;
622 fBCuts[5] = cut5;
623 fBCuts[6] = cut6;
624 fBCuts[7] = cut7;
625 fBCuts[8] = cut8;
626
627 return;
628}
629//-----------------------------------------------------------------------------
630void AliBtoJPSItoEleAnalysis::SetBCuts(const Double_t cuts[9]) {
631 // Set the cuts for B selection
632
633 for(Int_t i=0; i<9; i++) fBCuts[i] = cuts[i];
634
635 return;
636}
637//-----------------------------------------------------------------------------
638Bool_t
639AliBtoJPSItoEleAnalysis::SingleTrkCuts(const AliESDtrack& trk, Double_t b) const {
640 // Check if track passes some kinematical cuts
641 // Magnetic field "b" (kG)
642
643 if(TMath::Abs(1./trk.GetParameter()[4]) < fPtCut)
644 return kFALSE;
645 if(TMath::Abs(10000.*trk.GetD(fV1[0],fV1[1],b)) < fd0Cut)
646 return kFALSE;
647 //select only tracks with the "combined PID"
648 UInt_t status = trk.GetStatus();
649 if ((status&AliESDtrack::kESDpid)==0) return kTRUE;
650 Double_t r[5];
651 trk.GetESDpid(r);
652 if(r[0] < fPidCut) return kFALSE;
653
654 return kTRUE;
655}
656//----------------------------------------------------------------------------
71d22fe8 657void AliBtoJPSItoEleAnalysis::MakeTracksRefFile(AliRun *mygAlice,
27e190fb 658 Int_t evFirst,Int_t evLast) const {
659 // Create a file with simulation info for the reconstructed tracks
660
661 TFile *outFile = TFile::Open("BTracksRefFile.root","recreate");
662 TFile *esdFile = TFile::Open("AliESDs.root");
663
71d22fe8 664 AliMC *mc = mygAlice->GetMCApp();
27e190fb 665
666 Int_t label;
667 TParticle *part;
668 TParticle *mumpart;
669 TParticle *gmumpart;
670 REFTRACK reftrk;
671
5e6f195c 672 AliESDEvent* event = new AliESDEvent();
27e190fb 673 TTree* tree = (TTree*) esdFile->Get("esdTree");
5e6f195c 674 event->ReadFromTree(tree);
27e190fb 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();
680
71d22fe8 681 mygAlice->GetEvent(ev);
27e190fb 682
683 Int_t nentr=(Int_t)event->GetNumberOfTracks();
684
685 // Tree for true track parameters
686 char ttname[100];
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");
691
692 for(Int_t i=0; i<nentr; i++) {
693 AliESDtrack *esdtrack = (AliESDtrack*)event->GetTrack(i);
694 label = TMath::Abs(esdtrack->GetLabel());
695
696 part = (TParticle*)mc->Particle(label);
697 reftrk.lab = label;
698 reftrk.pdg = part->GetPdgCode();
699 reftrk.mumlab = part->GetFirstMother();
700 if(part->GetFirstMother()>=0) {
71d22fe8 701 mumpart = (TParticle*)mygAlice->GetMCApp()->Particle(part->GetFirstMother());
27e190fb 702 reftrk.mumpdg = mumpart->GetPdgCode();
703 reftrk.mumprongs = mumpart->GetNDaughters();
704 reftrk.gmumlab = mumpart->GetFirstMother();
705 if(mumpart->GetFirstMother()>=0) {
71d22fe8 706 gmumpart = (TParticle*)mygAlice->GetMCApp()->Particle(mumpart->GetFirstMother());
27e190fb 707 reftrk.gmumpdg = gmumpart->GetPdgCode();
708 }
709 } else {
710 reftrk.mumpdg=-1;
711 reftrk.mumprongs=-1;
712 reftrk.gmumpdg=-1;
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!
715 }
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();
722
723 reftree->Fill();
724
725 } // loop on tracks
726
727 outFile->cd();
728 reftree->Write();
729
730 delete reftree;
731 } // loop on events
732
733 esdFile->Close();
734 outFile->Close();
735
736 return;
737}
738//-----------------------------------------------------------------------------
739void AliBtoJPSItoEleAnalysis::SimulationInfo(TTree *treeBin,TTree *treeBout) const {
740 // add pdg codes to candidate decay tracks (for sim)
741
742 TString refFileName("BTracksRefFile.root");
743 if(fSim && gSystem->AccessPathName(refFileName.Data(),kFileExists)) {
744 printf("AliBtoJPSItoEleAnalysis::SimulationInfo: no reference file found!\n");
745 return;
746 }
747 TFile *refFile = TFile::Open(refFileName.Data());
748
749 Char_t refTreeName[100];
750 Int_t event;
751 Int_t pdg[2],mumpdg[2],mumlab[2],gmumpdg[2],gmumlab[2];
752 REFTRACK reftrk;
753
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);
758
759 AliBtoJPSItoEle *theB = 0;
760 treeBin->SetBranchAddress("BtoJPSItoEle",&theB);
761 treeBout->Branch("BtoJPSItoEle","AliBtoJPSItoEle",&theB,200000,0);
762
763 Int_t entries = (Int_t)treeBin->GetEntries();
764
765 for(Int_t i=0; i<entries; i++) {
766 if(i%100==0) printf(" done %d candidates of %d\n",i,entries);
767
768 treeBin->GetEvent(i);
769 event = theB->EventNo();
770
771 if(event==0) { // always true for Pb-Pb (avoid to read-in tree every time)
772 refTree0->GetEvent(theB->GetTrkNum(0));
773 pdg[0] = reftrk.pdg;
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));
779 pdg[1] = reftrk.pdg;
780 mumpdg[1] = reftrk.mumpdg;
781 mumlab[1] = reftrk.mumlab;
782 gmumpdg[1] = reftrk.gmumpdg;
783 gmumlab[1] = reftrk.gmumlab;
784 } else {
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));
789 pdg[0] = reftrk.pdg;
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));
795 pdg[1] = reftrk.pdg;
796 mumpdg[1] = reftrk.mumpdg;
797 mumlab[1] = reftrk.mumlab;
798 gmumpdg[1] = reftrk.gmumpdg;
799 gmumlab[1] = reftrk.gmumlab;
800 delete refTree;
801 }
802
803 theB->SetPdgCodes(pdg);
804 theB->SetMumPdgCodes(mumpdg);
805 theB->SetGMumPdgCodes(gmumpdg);
806
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
831 ) &&
832 mumpdg[0]==443 && mumpdg[1]== 443 &&
833 mumlab[0]==mumlab[1] &&
834 reftrk.mumprongs==2 &&
835 pdg[0]==-11 && pdg[1]==11
836 ) theB->SetSignal();
837
838 else if ( // here consider the case of primary J/psi
27e190fb 839 mumpdg[0]==443 && mumpdg[1]== 443 &&
0ac6b6e8 840 pdg[0]==-11 && pdg[1]==11 &&
27e190fb 841 mumlab[0]==mumlab[1] &&
842 reftrk.mumprongs==2 &&
0ac6b6e8 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)
850 )
27e190fb 851 ) theB->SetJpsiPrimary();
852
853 // if(!fOnlySignal || theB->IsSignal()) treeBout->Fill();
854
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();
860 }
861
862 delete refTree0;
863
864 refFile->Close();
865
866 return;
867}
868
869
870
871
872