memory consumption measures only memory used by preprocessor
[u/mrichter/AliRoot.git] / PWG3 / vertexingOld / AliBtoJPSItoEleAnalysis.cxx
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>
30 #include "AliESDEvent.h"
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"
40 #include "AliKFParticleBase.h"
41 #include "AliKFParticle.h"
42 #include "AliKFVertex.h"
43
44 typedef 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
56 ClassImp(AliBtoJPSItoEleAnalysis)
57
58 //----------------------------------------------------------------------------
59 AliBtoJPSItoEleAnalysis::AliBtoJPSItoEleAnalysis():
60 fVertexOnTheFly(kFALSE),
61 fSim(kFALSE),
62 fOnlySignal(kFALSE),
63 fOnlyPrimaryJpsi(kFALSE),
64 fPID("TRDTPCparam"),
65 fPtCut(0.),
66 fd0Cut(0.), 
67 fMassCut(1000.),
68 fPidCut(0.),
69 //fKFPrimVertex(kFALSE),
70 //fKFTopConstr(kFALSE),
71 fKFSecondVertex(kFALSE)
72 {
73   // Default constructor
74
75   SetBCuts();
76   SetVertex1();
77 }
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
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 //----------------------------------------------------------------------------
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;
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 //----------------------------------------------------------------------------
150 void 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)
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");
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());
196   AliESDEvent* event = new AliESDEvent();
197   TTree* tree = (TTree*) esdFile->Get("esdTree");
198   if(!tree) {
199     Error("FindCandidatesESD", "no ESD tree found");
200     return;
201   }
202   event->ReadFromTree(tree);
203
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
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
260 /*
261 //===================== PRIMARY VERTEX USING KF METHODS ==========================//
262
263 if (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
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
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         }
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
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();
403             delete vertex1onfly;        
404           }
405         }
406        
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
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 //-----------------------------------------------------------------------------
521 void 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 //-----------------------------------------------------------------------------
550 Bool_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 //-----------------------------------------------------------------------------
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
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 //-----------------------------------------------------------------------------
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
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 //-----------------------------------------------------------------------------
630 void 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 //-----------------------------------------------------------------------------
638 Bool_t 
639 AliBtoJPSItoEleAnalysis::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 //----------------------------------------------------------------------------
657 void AliBtoJPSItoEleAnalysis::MakeTracksRefFile(AliRun *mygAlice,
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
664   AliMC *mc = mygAlice->GetMCApp();
665   
666   Int_t      label;
667   TParticle *part;  
668   TParticle *mumpart;
669   TParticle *gmumpart;
670   REFTRACK   reftrk;
671   
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();
680
681     mygAlice->GetEvent(ev);
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) {
701         mumpart = (TParticle*)mygAlice->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*)mygAlice->GetMCApp()->Particle(mumpart->GetFirstMother());
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 //-----------------------------------------------------------------------------
739 void 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 
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)
850        )
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